Bugfix/pseudo rh #602

ClaraDraper-NOAA · 2023-07-28T20:20:07Z

Description

The EnKF uses the pseudo_RH=True option, which switches the humidity analysis variable to q/qsat (pseudo_RH=False would use q)
The default (and in operations) settings also includes use use_qsatensmean=True, which uses q/ensmean(qsat) for each ensemble member.

If use_qsatensmean is true, in controlvec.f90 read_control (then write_control) the q for each ensemble member are divided (then multiplied) by the same number (ensemble mean qsat) at each location before (then after) the state update. This multipilcation/division of each ensemble member by the same number has no net effect, and use_qsatensmean = true reduces to being the same as use pseudo_rh = false (using q as the control variable). Experiments show identical output from both options. See example in slides:

https://docs.google.com/presentation/d/1mLenEtuoN_9PI0--BEyIHN-E3rUW2WwcsHgzgqi2r0M/edit?usp=sharing

This PR then removes the use_qsatensmean option (including from the namelists). This required adding code to write out the ensemble mean for the pseudo_RH option without relying on the ensemble mean qsat.

In the PR I have the changed the pseudo_RH option to False in the regression test, which amounts to no science change. However, True would be better option and is what we thought we were doing up to now (more consistent with what is done in the Var solver, less heterogeneity in ensemble covariances with temperature). I suggest that we switch this going forward.

I ran both cases (but also decreased the obs error) in an EnKF experiment. There was no significant impact on skill. See slides linked above.

Note: @jswhit introduced the use_qsatensmean option to avoid problems with very small qsat at some locations. Could this become a problem again if we revert back to pseudo_RH=True.

Fixes #561

Type of change

Please delete options that are not relevant.

[ x] Bug fix (non-breaking change which fixes an issue)
New feature (non-breaking change which adds functionality)
Breaking change (fix or feature that would cause existing functionality to not work as expected)
This change requires a documentation update

How Has This Been Tested?

See experiments in slides above.
All regerssion tests have passed on hera.

Checklist

[x ] My code follows the style guidelines of this project
[ x] I have performed a self-review of my own code
[x ] I have commented my code, particularly in hard-to-understand areas
[x ] New and existing tests pass with my changes
[x ] Any dependent changes have been merged and published

DUE DATE for this PR is 9/8/2023. If this PR is not merged into develop by this date, the PR will be closed and returned to the developer.

CoryMartin-NOAA

Code changes look fine, but I will need to test in the global-workflow since the workflow scripts uses use_qsatensmean=.true. in a script that will need to be removed.

jswhit · 2023-07-31T15:32:16Z

I do recall there being crashes when pseudo_RH was True (before use_qsatmean was added). I wonder if we could/should multiply by the qsat computed from the updated p,T after the EnKF update when pseudo_RH=T? In that case use_qsatmean=T would not be the same as pseudo_RH=F.

ClaraDraper-NOAA · 2023-07-31T17:31:17Z

@jswhit . I've been thinking about your suggestion, since I do plan to introduce the use of qsat_anal when we scale the RH updated (from pseudo_RH=True) back to q. What your suggestion would effectively do is perform the update on q, then multip the update by qsat_back / qsat_anal. By "perform the update on q", I mean the ensemble correlations would effectively be with q not RH (since every ensemble member is divided by the same qsat). For the assimilation of screen level obs the correlations between q and T are very heterogenous (to the extend that they the sign varies), so using q will likely not give good q and T analyses (for other observations too - but the problem is more acute near the surface). We do need to address the crashing though - I think you said this was happening when q (or spread in q) was too small? Can we do something else to protect against that?

Depending on your answer to the above, I can re-introduce the use_qsatensmean code so that we have it in case we need it. This PR would then reduce to adding the code to write out the ensemble mean when use_qsatensmean is false.

jswhit · 2023-07-31T17:38:20Z

@ClaraDraper-NOAA since you have verified that use_qsatmean=T/pseudo_RH=T is the same as pseudo_RH=F I'm OK with removing that option (with the caveat that global_workflow will need to be updated). Sounds like the issue of how to renormalize q might be rather involved, and is best done separately.

ClaraDraper-NOAA · 2023-07-31T17:44:30Z

@ClaraDraper-NOAA since you have verified that use_qsatmean=T/pseudo_RH=T is the same as pseudo_RH=F I'm OK with removing that option (with the caveat that global_workflow will need to be updated). Sounds like the issue of how to renormalize q might be rather involved, and is best done separately.

Yes, good point.Are you OK with merging this as it is then?

jswhit · 2023-07-31T17:51:31Z

@ClaraDraper-NOAA since you have verified that use_qsatmean=T/pseudo_RH=T is the same as pseudo_RH=F I'm OK with removing that option (with the caveat that global_workflow will need to be updated). Sounds like the issue of how to renormalize q might be rather involved, and is best done separately.

Yes, good point.Are you OK with merging this as it is then?

started a review with a couple of comments

jswhit · 2023-07-31T17:52:36Z

May also have to revise regression tests if they have use_qsatmean=T

ClaraDraper-NOAA · 2023-07-31T18:00:04Z

I've removed use_qsatensmean from the regression tests, and switched pseudo_RH from T to F - so all tests pass. Switching to the recommended pseudo_RH= True would fail the current tests, and require a new baseline be generated.

ClaraDraper-NOAA · 2023-08-02T22:57:21Z

@CoryMartin-NOAA @jswhit Are there any more issues that I should address for this? If not, can you both please approve it?

jswhit · 2023-08-02T23:41:10Z

@ClaraDraper-NOAA did you see my review comments above?

CoryMartin-NOAA · 2023-08-03T12:37:50Z

@ClaraDraper-NOAA I am running a few days of cycles in our workflow just to confirm things are working as expected. I should be able to see results today or tomorrow and will approve it then. Thanks for making these changes.

CoryMartin-NOAA · 2023-08-03T15:21:54Z

@ClaraDraper-NOAA just to clarify, if I run with PSEUDO_RH=.false. in this branch, and .true. in develop, with USE_QSATENSMEAN=.true., the results should be identical, or are they expected to be different? I ask because I am seeing changes in the q ensemble mean increment.

ClaraDraper-NOAA · 2023-08-03T15:29:54Z

@ClaraDraper-NOAA just to clarify, if I run with PSEUDO_RH=.false. in this branch, and .true. in develop, with USE_QSATENSMEAN=.true., the results should be identical, or are they expected to be different? I ask because I am seeing changes in the q ensemble mean increment.

Yes, I'd expect those to the same. Are you using the write ensemble mean option within the GSI EnKF, or using a different program to calculate the mean of the ensemble members after they're written out?

CoryMartin-NOAA · 2023-08-03T15:30:37Z

@ClaraDraper-NOAA disregard, it is because the workflow checks out a GSI version by default that does not include #568 yet. This tracer clipping difference is likely the cause, but let me verify.

ClaraDraper-NOAA · 2023-08-03T19:07:01Z

src/enkf/params.f90

@@ -681,10 +675,6 @@ subroutine read_namelist()
       letkf_flag) then
     print *,'warning: no time localization in LETKF!'
   endif
-   if ((write_ensmean .and. pseudo_rh) .and. .not. use_qsatensmean) then


@jswhit Responding to "Do we need to raise on error and stop if write_ensmean=T and pseudo_RH=T?".

No, I moved the re-normalization by qsat to the top of the writecontrol routine, so it's now OK to use write_ensmean=T and psuedo_RH=T. Note: I got slightly different results when doing the normalization here then writing out the mean, than when I calculated the mean from the ensemble q/qsat output.

ClaraDraper-NOAA · 2023-08-03T19:13:50Z

Pushed a change with the qsatmean and qsat_tmp declarations removed at Jeff's request.

jswhit · 2023-07-31T17:39:49Z

src/enkf/params.f90

-   if ((write_ensmean .and. pseudo_rh) .and. .not. use_qsatensmean) then
-      print *,'write_ensmean=T requires use_qsatensmean=T when pseudo_rh=T'
-      call stop2(19)
-   endif


Do we need to raise an error and stop if write_ensmean=T and pseudo_RH=T (just remove the .and. not. use_qsatensmean part)?

jswhit · 2023-07-31T17:42:35Z

src/enkf/controlvec.f90

@@ -229,56 +229,18 @@ subroutine read_control()
   end if
   !print *,'min/max qsat',nanal,'=',minval(qsat),maxval(qsat)
   q_ind = getindex(cvars3d, 'q')
-   if (use_qsatensmean .and. q_ind>0 ) then
-       allocate(qsatmean(npts,nlevs,nbackgrounds))


Can also remove the variable declaration for qsatmean and qsat_tmp

jswhit · 2023-08-03T19:29:45Z

src/enkf/params.f90

@@ -681,10 +675,6 @@ subroutine read_namelist()
       letkf_flag) then
     print *,'warning: no time localization in LETKF!'
   endif
-   if ((write_ensmean .and. pseudo_rh) .and. .not. use_qsatensmean) then


CoryMartin-NOAA

I ran some low resolution cycling test to make sure everything is working as intended, looks good to me!

NOAA-EMC/GSI#602 was merged in without the subsequent script changes needed in the workflow. This PR fixes that.

ClaraDraper-NOAA added 2 commits July 27, 2023 22:47

Commented out use_qsatensmean code. Generates expected analysis.

a533251

completely removed use_qsatensmean option.

b459667

RussTreadon-NOAA assigned ClaraDraper-NOAA Jul 28, 2023

RussTreadon-NOAA requested review from jswhit2 and CoryMartin-NOAA July 28, 2023 20:30

ClaraDraper-NOAA marked this pull request as ready for review July 28, 2023 21:51

CoryMartin-NOAA reviewed Jul 31, 2023

View reviewed changes

ClaraDraper-NOAA commented Aug 3, 2023

View reviewed changes

removed qsatmean and qsat_tmp declarations.

79e8c67

jswhit approved these changes Aug 3, 2023

View reviewed changes

jswhit2 approved these changes Aug 3, 2023

View reviewed changes

CoryMartin-NOAA approved these changes Aug 4, 2023

View reviewed changes

CoryMartin-NOAA merged commit accb07e into NOAA-EMC:develop Aug 7, 2023
4 checks passed

guoqing-noaa mentioned this pull request Aug 23, 2023

Gsi fed #590

Merged

9 tasks

CoryMartin-NOAA mentioned this pull request Sep 1, 2023

Fix eupd after GSI PR 602 NOAA-EMC/global-workflow#1835

Merged

7 tasks

WalterKolczynski-NOAA pushed a commit to NOAA-EMC/global-workflow that referenced this pull request Sep 8, 2023

Fix eupd after GSI PR 602 (#1835)

af876c2

NOAA-EMC/GSI#602 was merged in without the subsequent script changes needed in the workflow. This PR fixes that.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Bugfix/pseudo rh #602

Bugfix/pseudo rh #602

ClaraDraper-NOAA commented Jul 28, 2023 •

edited

Loading

CoryMartin-NOAA left a comment

jswhit commented Jul 31, 2023 •

edited

Loading

ClaraDraper-NOAA commented Jul 31, 2023 •

edited

Loading

jswhit commented Jul 31, 2023

ClaraDraper-NOAA commented Jul 31, 2023

jswhit commented Jul 31, 2023 •

edited

Loading

jswhit commented Jul 31, 2023

ClaraDraper-NOAA commented Jul 31, 2023

ClaraDraper-NOAA commented Aug 2, 2023

jswhit commented Aug 2, 2023

CoryMartin-NOAA commented Aug 3, 2023

CoryMartin-NOAA commented Aug 3, 2023

ClaraDraper-NOAA commented Aug 3, 2023

CoryMartin-NOAA commented Aug 3, 2023

ClaraDraper-NOAA Aug 3, 2023

jswhit Aug 3, 2023

ClaraDraper-NOAA commented Aug 3, 2023

jswhit Jul 31, 2023

jswhit Aug 3, 2023

jswhit Jul 31, 2023

jswhit Aug 3, 2023

CoryMartin-NOAA left a comment

Bugfix/pseudo rh #602

Bugfix/pseudo rh #602

Conversation

ClaraDraper-NOAA commented Jul 28, 2023 • edited Loading

CoryMartin-NOAA left a comment

Choose a reason for hiding this comment

jswhit commented Jul 31, 2023 • edited Loading

ClaraDraper-NOAA commented Jul 31, 2023 • edited Loading

jswhit commented Jul 31, 2023

ClaraDraper-NOAA commented Jul 31, 2023

jswhit commented Jul 31, 2023 • edited Loading

jswhit commented Jul 31, 2023

ClaraDraper-NOAA commented Jul 31, 2023

ClaraDraper-NOAA commented Aug 2, 2023

jswhit commented Aug 2, 2023

CoryMartin-NOAA commented Aug 3, 2023

CoryMartin-NOAA commented Aug 3, 2023

ClaraDraper-NOAA commented Aug 3, 2023

CoryMartin-NOAA commented Aug 3, 2023

ClaraDraper-NOAA Aug 3, 2023

Choose a reason for hiding this comment

jswhit Aug 3, 2023

Choose a reason for hiding this comment

ClaraDraper-NOAA commented Aug 3, 2023

jswhit Jul 31, 2023

Choose a reason for hiding this comment

jswhit Aug 3, 2023

Choose a reason for hiding this comment

jswhit Jul 31, 2023

Choose a reason for hiding this comment

jswhit Aug 3, 2023

Choose a reason for hiding this comment

CoryMartin-NOAA left a comment

Choose a reason for hiding this comment

ClaraDraper-NOAA commented Jul 28, 2023 •

edited

Loading

jswhit commented Jul 31, 2023 •

edited

Loading

ClaraDraper-NOAA commented Jul 31, 2023 •

edited

Loading

jswhit commented Jul 31, 2023 •

edited

Loading