GitHub Issue NOAA-EMC/GSI#468 Enhancements to SDL and VDL for simultaneous multiscale EnVar and parallel ensemble IO for EnVar for FV3-LAM

[Wangy1111]

The following capabilities developed by OU MAP lab are included

(1) Further development for simultaneous multiscale EnVar for both global and regional DA
(1a) spatial scale-dependent localization (SDL; contributed by Ting Lei and Daryl Kleist/EMC) is implemented in EnVar as described in Huang et al 2021, MWR for the global NWP application.
(1b) variable-dependent localization (VDL) method by Wang and Wang 2022, JAMES is implemented in EnVar.

(2)Development of parallel ensemble IO for EnVar for FV3-LAM 
Implement an approach to simultaneously read in all ensemble members for EnVar. Specifically, parallel ensemble IO for both conventional and radar EnVar for FV3-LAM is implemented by reading in all ensemble members simultaneously.

(3) Direct assimilation of radar reflectivity for EnVar for RRFS
The direct radar reflectivity assimilation approach by Wang and Wang 2017, MWR is implemented and tested for FV3-LAM.

Fixes #468

[hu5970]

I ran regression test suite on Hera. 18 out of 19 cases passed.

The one crashed is fv3lam netcdf case, which report the error information on missing physics file.
We need to update the regression test case for fv3lam case. I will work on this.

[TingLei-NOAA]

@hu5970 Thanks for this finding. However, should we let the program don't read physics files when they are not needed?
I think the physics files are only needed when dbz is assimilated , right?

[TingLei-NOAA]

@hu5970 I will open an issue on regression tests soon or you can go ahead too:).

[Wangy1111]

Thank you for your comments. @hu5970 @TingLei-NOAA We should have the codes to let this branch read/write physics files only when if_model_dbz is set to .true. (direct assimilation of dbz). Please let us know if you need anything from us to proceed with the regression tests.

[TingLei-NOAA]

@Wangy1111 Thanks for your update. Yes, it will be better to use if_model_dbz to control if physics files are needed.
About the regression test issue, sorry, my comment should be replied in another PR when Ming and myself were talking about opening an issue to describe/update on how to run regression tests. I am not sure if @hu5970 would hope you and your colleagues would attach regression tests for this PR, while they would be surely helpful for developers/users in the future.

[daviddowellNOAA]

@Wangy1111 I'm testing the GSI-EnKF code represented by this pull request for the experimental RRFS_B. Currently, RRFS_B assimilates conventional (non-zero static B) and then radar-reflectivity observations (zero static B) in separate executions of GSI. I would like to propose changing line 2014 of gsimod.F90 from

      if ( index(dtype(i), 'dbz') /= 0 )then

to

      if ( if_model_dbz .and. ( index(dtype(i), 'dbz') /= 0 ) ) then

This change would allow the conventional and radar DA in RRFS_B to use the same convinfo file, and avoid a failure from the non-zero static B in the conventional DA. What are your thoughts on this proposed change?

[wangyongm]

@Wangy1111 I'm testing the GSI-EnKF code represented by this pull request for the experimental RRFS_B. Currently, RRFS_B assimilates conventional (non-zero static B) and then radar-reflectivity observations (zero static B) in separate executions of GSI. I would like to propose changing line 2014 of gsimod.F90 from

      if ( index(dtype(i), 'dbz') /= 0 )then

to

      if ( if_model_dbz .and. ( index(dtype(i), 'dbz') /= 0 ) ) then

This change would allow the conventional and radar DA in RRFS_B to use the same convinfo file, and avoid a failure from the non-zero static B in the conventional DA. What are your thoughts on this proposed change?

@daviddowellNOAA Thanks for your comments. We have changed it as you suggested.

[hu5970]

Yes, it will be better to use if_model_dbz to control if physics files are needed.

@Wangy1111 Yongming, are you OK to make this change to the code also? Thanks, Ming

[Wangy1111]

Yes, it will be better to use if_model_dbz to control if physics files are needed.

@Wangy1111 Yongming, are you OK to make this change to the code also? Thanks, Ming

@hu5970 I forgot to mention that using if_model_dbz has been added to the code. Sorry about that.
In the original code, the default if_model_dbz is .true. Such a setting may lead to the failure of the regression tests when if_model_dbz is not specified to be .false. In the latest code, the default if_model_dbz is changed to .false. Therefore, the physics file won't be required when if_model_dbz is not specified to be .true. or not used in the namelist. Thanks for your comment.

[shoyokota]

@Wangy1111 I'm still on the way of the review, but I let you know what I noticed at this moment.

• Could you explain how to set the options you added (l_sum_spc_weights, vdl_scale, vloc_varlist, smooth_scales_num, global_spectral_filter_sd, assign_vdl_nml)? It is better to add comments in gsimod.F90 and hybrid_ensemble_parameters.f90. In particular, "vdl_scale" and "vloc_varlist" are difficult to understand.
• "l_sum_spc_weights" is not probably needed because it is the same as "i_ensloccov4scl".
• "smooth_scales_num" is not probably needed because it can be obtained as "smooth_scales_num=naensloc-naensgrp".
• I think "ngvarloc=1" should be forced in the case of "assign_vdl_nml=.true."
• Refactoring you did last week looks reverted by the update you did yesterday.

@hu5970 I confirmed that this branch brings different results from the current RRFS branch in dbz assimilation. It is probably caused by a bug of the current RRFS branch ( https://github.com/NOAA-GSL/GSI/blob/feature/rrfs_dev/src/gsi/gsi_rfv3io_mod.f90#L2207 ). I guess it should be fixed from uu2d(4:nxcase-3,4:nycase-3) to uu2d(1:nxcase,1:nycase).

@TingLei-NOAA This branch does not include the modification to output netCDF diag files. ( NOAA-GSL/GSI#20 ) Is it better to add the modification also in this branch?

[Wangy1111]

@Wangy1111 I'm still on the way of the review, but I let you know what I noticed at this moment.

@hu5970 I confirmed that this branch brings different results from the current RRFS branch in dbz assimilation. It is probably caused by a bug of the current RRFS branch ( https://github.com/NOAA-GSL/GSI/blob/feature/rrfs_dev/src/gsi/gsi_rfv3io_mod.f90#L2207 ). I guess it should be fixed from uu2d(4:nxcase-3,4:nycase-3) to uu2d(1:nxcase,1:nycase).

@shoyokota Thanks for your comments and review.
We will first explain the above question. To achieve the LBC update, the grid size of fv3_phyvar is smaller than that of fv3_dynvar and fv3_tracer. Therefore, uu2d(4:nxcase-3,4:nycase-3) is used. We also notice that someone may use the same grid size for all files, so uu2d(1:nxcase,1:nycase) will be needed. Thus, this version code allows GSI to use uu2d(1:nxcase,1:nycase) or uu2d(4:nxcase-3,4:nycase-3) automatically depending on the grid size of fv3_phyvar.

We will address the remaining comments soon.

[shoyokota]

Thank you for the clarification. Yes, coding in the branch of this PR looks correct. I just wanted to point out Ming that coding in the current RRFS branch ( https://github.com/NOAA-GSL/GSI/tree/feature/rrfs_dev/src/gsi/gsi_rfv3io_mod.f90 ) is not correct. In my recognision, the current RRFS uses the same grid size for all files but uses uu2d(4:nxcase-3,4:nycase-3) probably by mistake.

[daviddowellNOAA]

@shoyokota Good catch regarding the dbz assimilation. My tests agree with yours. Specifically,

(1) For conventional GSI EnVar data assimilation with RRFS background ensemble, I get identical results when I use the RRFS GSI and the new GSI represented by this pull request (PR 504).

(2) For radar-reflectivity DA, I get different results with the two versions of GSI. Here are some analysis difference plots (analysis with PR 504 GSI minus analysis with RRFS GSI) for the November 4 case. The plots are for T in the southern plains (model level 64) and U in the northwest (model level 54).

As you suggested, we'll need the capability to output netcdf diag files in order to test EnKF.

[Wangy1111]

@Wangy1111 I'm still on the way of the review, but I let you know what I noticed at this moment.

* Could you explain how to set the options you added (l_sum_spc_weights, vdl_scale, vloc_varlist, smooth_scales_num, global_spectral_filter_sd, assign_vdl_nml)? It is better to add comments in gsimod.F90 and hybrid_ensemble_parameters.f90. In particular, "vdl_scale" and "vloc_varlist" are difficult to understand.

* "l_sum_spc_weights" is not probably needed because it is the same as "i_ensloccov4scl".

* "smooth_scales_num" is not probably needed because it can be obtained as "smooth_scales_num=naensloc-naensgrp".

* I think "ngvarloc=1" should be forced in the case of "assign_vdl_nml=.true."

* Refactoring you did last week looks reverted by the update you did yesterday.

@shoyokota We have addressed the above comments in the code.
As you suggested, l_sum_spc_weights and smooth_scales_num have been removed from the namelist part. The option ngvarloc is forced to 1 when assign_vdl_nml=.true.
We also added explanations for these new options in gsimod.F90.

Thank you.

[shoyokota]

I confirmed that the branch of this PR and the current RRFS branch bring the same results in dbz assimilation if using same "uu2d". I also confirmed that "assign_vdl_nml=.true., nsclgrp=2, ngvarloc=1" and "assign_vdl_nml=.false., nsclgrp=1, ngvarloc=2" can bring the same results. It is great!

@Wangy1111 Could you check the additional comments as follows?

• The dimension of "vdl_scale" is redundant. For example, vdl_scale(2:3) and vdl_scale(5:6) are not used in the case of "vdl_scale=3,3,3,3,3,3". In this case, the dimension can be reduced like "vdl_scale=3,3" (the number of dimension is the number of scales of SDL and each component is the number of variable groups of each scale of SDL). Could you modify to use such smaller "vdl_scale"?
• gsimod.F90#L1872 should be modified from "if(naensloc<naensgrp+nsclgrp-1)" to "if(.not.assign_vdl_nml .and. naensloc<naensgrp+nsclgrp-1)".
• I guess this VDL does not work correctly in the case of "ntlevs_ens>1" because "en_pertstmp(n)%valuesr4=en_pertstmp(n)%valuesr4-en_pertstmp1(n)%valuesr4" is in the loop of "do m=1,ntlevs_ens". I think "en_pertstmp(n)" and "en_pertstmp1(n)" should be "en_pertstmp(n,m)" and "en_pertstmp1(n,m)", respectively.
• Could you modify gsimod.F90#L1383 from "recursive filter" to "Gaussian filter"? The horizontal recursive filter is not used for the global application.
• Is "spcwgt_params( : , 2 )" required? It is not used in apply_scaledepwgts.f90#L65.
• "work1", "outwork", "iflg", "igrp", and "nensid" are not needed in subroutine apply_scaledepwgts.
• In general_specmod.f90, "JMAX" and "L" are still uppercase and "L=0" is not nesessary.

[RussTreadon-NOAA]

@Wangy1111 , it is not necessary to force push changes to your forked develop. We can squash merge this PR into the authoritative repository when it is finally approved. We loose history information with force pushes. This is contrary to a core tenant of version control. We want a history of changes to a file. Despite the many changes to apply_scaledepwgts.f90, a check of the file history only shows one revision, d7be093.

[Wangy1111]

@Wangy1111 , it is not necessary to force push changes to your forked develop. We can squash merge this PR into the authoritative repository when it is finally approved. We loose history information with force pushes. This is contrary to a core tenant of version control. We want a history of changes to a file. Despite the many changes to apply_scaledepwgts.f90, a check of the file history only shows one revision, d7be093.

@RussTreadon-NOAA Thanks for your information.
I did not know the issue you mentioned. Is there anything I can do to remedy it?

[RussTreadon-NOAA]

@Wangy1111 , I'm not a git guru but my understanding is that git --force push forcibly overwrites the commit history. Previous commit history is lost.

[Wangy1111]

@Wangy1111 , I'm not a git guru but my understanding is that git --force push forcibly overwrites the commit history. Previous commit history is lost.

I use git commit --amend and git push origin develop --force to commit the changes.

[Wangy1111]

@Wangy1111 , it is not necessary to force push changes to your forked develop. We can squash merge this PR into the authoritative repository when it is finally approved. We loose history information with force pushes. This is contrary to a core tenant of version control. We want a history of changes to a file. Despite the many changes to apply_scaledepwgts.f90, a check of the file history only shows one revision, d7be093.

@RussTreadon-NOAA Could you please let us know if "the file history" here means the history from the other PR's previous commits OR the history only from this PR's commit?
If you mean that from the other PR's previous commits, the history still exists. However, the history from this PR's commit was indeed overwritten.

For apply_scaledepwgts.f90, this is a new file, so only one revision is shown.

Thank you for your time.

[RussTreadon-NOAA]

I am referring to the entire development history of this file. Yes, this is a new file, but it has gone through several revisions in your development. As you state the development history of this file in your PR has been reduce to only the most recent change. The full development history of this file has been lost.

[Wangy1111]

I am referring to the entire development history of this file. Yes, this is a new file, but it has gone through several revisions in your development. As you state the development history of this file in your PR has been reduce to only the most recent change. The full development history of this file has been lost.

@RussTreadon-NOAA Thank you for the clarification.

[TingLei-NOAA]

@Wangy1111 when the capabilities in this PR were all under the name of "The following capabilities developed by OU MAP lab ", I hope EMC's major contribution to the GSI coding for the global SDL application (1a in the summary) could be stated.

[hu5970]

@Wangy1111 Could you sync with EMC develop branch again? There is a PR merged during weekend for fixing the regression tests. This will not impact your PR. But need to be synced to conduct regression tests.

[Wangy1111]

@Wangy1111 Could you sync with EMC develop branch again? There is a PR merged during weekend for fixing the regression tests. This will not impact your PR. But need to be synced to conduct regression tests.

@hu5970 Done.

[Wangy1111]

@Wangy1111 when the capabilities in this PR were all under the name of "The following capabilities developed by OU MAP lab ", I hope EMC's major contribution to the GSI coding for the global SDL application (1a in the summary) could be stated.

Added, see above.

[hu5970]

@TingLei-NOAA @CatherineThomas-NOAA Ting and Cathy, could you review the PR and approve it again if the PR is OK.
I will run regression tests soon. Thanks, Ming

[hu5970]

The regression tests finished on Orion:

1/9 Test #9: [=[global_enkf]=] ................   Passed  486.80 sec
2/9 Test #8: [=[netcdf_fv3_regional]=] ........   Passed  603.20 sec
3/9 Test #7: [=[rrfs_3denvar_glbens]=] ........***Failed  605.81 sec
4/9 Test #4: [=[hwrf_nmm_d2]=] ................   Passed  607.14 sec
5/9 Test #5: [=[hwrf_nmm_d3]=] ................***Failed  734.72 sec
6/9 Test #6: [=[rtma]=] .......................***Failed  1208.81 sec
7/9 Test #3: [=[global_4denvar]=] .............   Passed  1561.95 sec
8/9 Test #2: [=[global_4dvar]=] ...............   Passed  1801.50 sec
9/9 Test #1: [=[global_3dvar]=] ...............   Passed  1861.77 sec

67% tests passed, 3 tests failed out of 9

Total Test time (real) = 1861.79 sec

The following tests FAILED:
	  5 - [=[hwrf_nmm_d3]=] (Failed)
	  6 - [=[rtma]=] (Failed)
	  7 - [=[rrfs_3denvar_glbens]=] (Failed)

For three failed cases, they all failed because of
"The case has Failed the scalability test"

They are not a fatal failure.

[RussTreadon-NOAA]

Hera ctests
NOAA-EMC/GSI:develop at e55a937 and Wangy1111/GSI:develop at d511d44 were installed on Hera and the suite of 9 regression (ctests) run. 8 of the 9 tests passed

Hera(hfe10):/scratch1/NCEPDEV/da/Russ.Treadon/git/gsi/pr504/build$ ctest -j 9
Test project /scratch1/NCEPDEV/da/Russ.Treadon/git/gsi/pr504/build
    Start 1: global_3dvar
    Start 2: global_4dvar
    Start 3: global_4denvar
    Start 4: hwrf_nmm_d2
    Start 5: hwrf_nmm_d3
    Start 6: rtma
    Start 7: rrfs_3denvar_glbens
    Start 8: netcdf_fv3_regional
    Start 9: global_enkf
1/9 Test #7: rrfs_3denvar_glbens ..............   Passed  849.26 sec
2/9 Test #8: netcdf_fv3_regional ..............   Passed  1505.97 sec
3/9 Test #4: hwrf_nmm_d2 ......................   Passed  1510.61 sec
4/9 Test #5: hwrf_nmm_d3 ......................***Failed  1520.41 sec
5/9 Test #9: global_enkf ......................   Passed  1602.62 sec
6/9 Test #2: global_4dvar .....................   Passed  2765.06 sec
7/9 Test #3: global_4denvar ...................   Passed  2867.59 sec
8/9 Test #6: rtma .............................   Passed  3493.32 sec
9/9 Test #1: global_3dvar .....................   Passed  3587.52 sec

89% tests passed, 1 tests failed out of 9

Total Test time (real) = 3587.53 sec

The following tests FAILED:
          5 - hwrf_nmm_d3 (Failed)
Errors while running CTest
Output from these tests are in: /scratch1/NCEPDEV/da/Russ.Treadon/git/gsi/pr504/build/Testing/Temporary/LastTest.log
Use "--rerun-failed --output-on-failure" to re-run the failed cases verbosely.

Examination of hwrf_nmm_d3 showed that it failed due to a wall time scalability test

The case has Failed the scalability test.
The slope for the update (10.731268 seconds per node) is less than that for the control (12.317250 seconds per node).

A check of the gsi.x wall times for the updat and contrl do not show any anomalies

Hera(hfe07):/scratch1/NCEPDEV/stmp2/Russ.Treadon/pr504/tmpreg_hwrf_nmm_d3$ grep "wall time" hwrf*/stdout
hwrf_nmm_d3_hiproc_contrl/stdout:The total amount of wall time                        = 62.380179
hwrf_nmm_d3_hiproc_updat/stdout:The total amount of wall time                        = 62.599070
hwrf_nmm_d3_loproc_contrl/stdout:The total amount of wall time                        = 74.697429
hwrf_nmm_d3_loproc_updat/stdout:The total amount of wall time                        = 69.753249

A rerun of the hwrf_nmm_d3 test on Hera passed

Hera(hfe10):/scratch1/NCEPDEV/da/Russ.Treadon/git/gsi/pr504/build$ ctest -R hwrf_nmm_d3
Test project /scratch1/NCEPDEV/da/Russ.Treadon/git/gsi/pr504/build
    Start 5: hwrf_nmm_d3
1/1 Test #5: hwrf_nmm_d3 ......................   Passed  623.37 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) = 623.39 sec

[hu5970]

The regression tests on WCOSS2 finishedL

1/9 Test #8: netcdf_fv3_regional ..............***Failed  1503.88 sec
2/9 Test #9: global_enkf ......................***Failed  1512.73 sec
3/9 Test #7: rrfs_3denvar_glbens ..............   Passed  1568.96 sec
4/9 Test #4: hwrf_nmm_d2 ......................   Passed  1872.95 sec
5/9 Test #5: hwrf_nmm_d3 ......................   Passed  1873.43 sec
6/9 Test #2: global_4dvar .....................   Passed  2162.62 sec
7/9 Test #6: rtma .............................   Passed  2531.94 sec
8/9 Test #3: global_4denvar ...................   Passed  2703.34 sec
9/9 Test #1: global_3dvar .....................   Passed  2707.39 sec

78% tests passed, 2 tests failed out of 9

Total Test time (real) = 2707.39 sec

The following tests FAILED:
	  8 - netcdf_fv3_regional (Failed)
	  9 - global_enkf (Failed)

The "netcdf_fv3_regional" failure is because of "scalability"
The "global_enkf" failure is because of "timethresh2"

Both are not fatal failures.

[TingLei-NOAA]

@Wangy1111 Thanks for the prompt response, Would you change from Ting Lei to Ting Lei and Daryl Kleist? Thanks.

[Wangy1111]

@Wangy1111 Thanks for the prompt response, Would you change from Ting Lei to Ting Lei and Daryl Kleist? Thanks.

Changed.

[TingLei-NOAA]

@Wangy1111 Please allow me make another suggestion. Would it be more accurate to say " as described by Huang et al 2021, MW..." rather than "following Huang et al 2021, MW..". EMC 's work on this was mainly done in 2018 and 2019.

[Wangy1111]

@Wangy1111 Please allow me make another suggestion. Would it be more accurate to say " as described by Huang et al 2021, MW..." rather than "following Huang et al 2021, MW..". EMC 's work on this was mainly done in 2018 and 2019.

Changed as suggested.

[hu5970]

Thank Ting, Sho, Cathy for reviewing this big PR thoroughly and Yongming for quick response to all comments.
All regression cases passed without fatal failure.