Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding I/O interfaces in GSI for Analysis of Significant Wave Height (HOWV) and near-surface Wind Gust (GUST) for WRF-ARW based 3DRTMA #835

Open
wants to merge 4 commits into
base: develop
Choose a base branch
from

Conversation

GangZhao-NOAA
Copy link
Contributor

@GangZhao-NOAA GangZhao-NOAA commented Feb 7, 2025

Description

Motivation:

  • The analysis of significant wave height (HOWV) and 10-m wind gust (GUST) are listed as standard product in the operational (2D)RTMA/URMA product suites. To keep providing these two products in next generation RTMA/URMA (ie, 3DRTMA) system, the capability to directly analyze wave height and wind gust must be added into GSI 3DVar analysis.

Method:

  • The kernel part of the analysis code (such as stphowv.f90, setuphowv.f90, inthowv.f90, gsi_howvOper.f90 and m_howvNode.f90 for HOWV, and similar code for GUST) for wave height and wind gust had already been introduced into GSI through the implementation of 2DRTMA. So most of the code work for analysis of HOWV and GUST in 3DVar, are mainly focused on the I/O interfaces for the firstguess provided by various models (i.e. reading background of HOWV/GUST from firstguess file and writing the updated HOWV/GUST back to the final analysis file).
  • The I/O interfaces for FV3-LAM model based firstguess had been done in previous work (see issue Adding I/O for Analysis of Significant Wave Height for 3DRTMA #601 and issue Adding I/O for direct analysis of near-surface wind gust for RRFS-based 3DRTMA #726). Since RRFS-3DRTMAv1 has been switched to be based on HRRR model, so in this work, some efforts are made to add the I/O interfaces in GSI for WRF-ARW model based firstguess.

Related Issues:

Resolves #834

Type of change

Please delete options that are not relevant.

  • New feature (non-breaking change which adds functionality)

How Has This Been Tested?

Checklist

  • My code follows the style guidelines of this project
  • I have performed a self-review of my own code
  • I have commented my code, particularly in hard-to-understand areas
  • New and existing tests pass with my changes
  • Any dependent changes have been merged and published
  • Regression tests on WCOSS2 (Cactus) and Hera (Rocky-8) : my updated GSI commit #47b3115) vs control/original GSI code (commit #92165a4)
    Here is the reports of the regression tests on WCOSS2 (Dogwood):
[gang.zhao@dlogin09:build] (feature/rtma3d_wavhgtwndgst_wrfarw)$ ctest -N
Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/ProdGSI_dev/gsi_dev/build
  Test #1: global_4denvar
  Test #2: rtma
  Test #3: rrfs_3denvar_rdasens
  Test #4: hafs_4denvar_glbens
  Test #5: hafs_3denvar_hybens
  Test #6: global_enkf

Total Tests: 6
[gang.zhao@dlogin09:build] (feature/rtma3d_wavhgtwndgst_wrfarw)$ ctest -j 6
Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/ProdGSI_dev/gsi_dev/build
    Start 1: global_4denvar
    Start 2: rtma
    Start 3: rrfs_3denvar_rdasens
    Start 4: hafs_4denvar_glbens
    Start 5: hafs_3denvar_hybens
    Start 6: global_enkf
1/6 Test #3: rrfs_3denvar_rdasens .............   Passed  1155.27 sec
2/6 Test #5: hafs_3denvar_hybens ..............   Passed  1586.49 sec
3/6 Test #4: hafs_4denvar_glbens ..............   Passed  1650.16 sec
4/6 Test #6: global_enkf ......................   Passed  1957.18 sec
5/6 Test #2: rtma .............................   Passed  3319.15 sec
6/6 Test #1: global_4denvar ...................   Passed  3546.50 sec

100% tests passed, 0 tests failed out of 6

Total Test time (real) = 3546.52 sec

Here is the reports of the regression tests on Hera:

(work-env2) [Gang.Zhao@hfe01:build] (feature/rtma3d_wavhgtwndgst_wrfarw)$ ctest -N
Test project /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build
  Test #1: global_4denvar
  Test #2: rtma
  Test #3: rrfs_3denvar_rdasens
  Test #4: hafs_4denvar_glbens
  Test #5: hafs_3denvar_hybens
  Test #6: global_enkf

Total Tests: 6
(work-env2) [Gang.Zhao@hfe01:build] (feature/rtma3d_wavhgtwndgst_wrfarw)$ ctest -j 6
Test project /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build
    Start 4: hafs_4denvar_glbens
    Start 1: global_4denvar
    Start 6: global_enkf
    Start 5: hafs_3denvar_hybens
    Start 2: rtma
    Start 3: rrfs_3denvar_rdasens
1/6 Test #3: rrfs_3denvar_rdasens .............***Failed  2218.19 sec
2/6 Test #5: hafs_3denvar_hybens ..............   Passed  3687.92 sec
3/6 Test #2: rtma .............................   Passed  3828.48 sec
4/6 Test #4: hafs_4denvar_glbens ..............   Passed  4112.07 sec
5/6 Test #1: global_4denvar ...................   Passed  4118.91 sec
6/6 Test #6: global_enkf ......................   Passed  4119.37 sec

83% tests passed, 1 tests failed out of 6

Total Test time (real) = 4119.38 sec

The following tests FAILED:
          3 - rrfs_3denvar_rdasens (Failed)
Errors while running CTest
Output from these tests are in: /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build/Testing/Temporary/LastTest.log
Use "--rerun-failed --output-on-failure" to re-run the failed cases verbosely.
(work-env2) [Gang.Zhao@hfe01:build] (feature/rtma3d_wavhgtwndgst_wrfarw)$ ctest -R rrfs_3denvar_rdasens
Test project /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build
    Start 3: rrfs_3denvar_rdasens
1/1 Test #3: rrfs_3denvar_rdasens .............   Passed  619.88 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) = 619.90 sec

! contribution to the total BE of howv, so the total BE of howv is actually
! just the reduced static BE of howv. If to make the analysis of howv
! in hyrbid run is as similar as the analysis of howv in pure 3dvar run,
! the static BE of howv used in hybrid run needs to be tuned (inflated actually).
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@GangZhao-NOAA Then, how can we do the inflation for static B for howv? Did I miss it?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@TingLei-NOAA
Hi Ting,

Sorry for my confusing comments in the source code.

What I mean is to give a larger background error for HOWV when running in hybrid EnVar analysis.

For example, in the pure 3dvar run of GSI the background error of HOWV is 0.30 meters. If we would like to get the similar analysis of HOWV in a hybrid EnVar run of GSI with the weight for ensemble covariances = 50%, because there is no ensemble for HOWV yet, and if we want to get the analysis of HOWV as similar as from pure 3DVar, we could NOT use 0.3 meters as the static background error, because the static background error variances is reduced by the weight (100%-50%). To compensate for this reduction in background error, we need to "inflate" 0.3 to some larger values.

My explanation comment in the source code is really misleading, I will clean off these comments before merging the code.

Thank you very much for pointing this out!

-Gang

@TingLei-NOAA
Copy link
Contributor

@GangZhao-NOAA Thanks. Your inline comments are good and I just wondered if there are already existing options for user to do the inflation. If there is no, that is ok. That comment will help users/developer understand the situation.

Copy link
Contributor

@TingLei-NOAA TingLei-NOAA left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for continual enhancements of GSI

@GangZhao-NOAA
Copy link
Contributor Author

@GangZhao-NOAA Thanks. Your inline comments are good and I just wondered if there are already existing options for user to do the inflation. If there is no, that is ok. That comment will help users/developer understand the situation.

@TingLei-NOAA Hi Ting, I realized that I should not use the word "inflat" here, since "covariance inflation" is somehow a specific term used in ensemble DA. I will change it to "increase" in my next PR. Since the process to increase the static error for HOWV and GUST are done outside of GSI code, not in GSI, so I'd better remove these comments in my next PR. Thank you again!

@@ -2085,25 +2119,25 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype)
end if
if ( laeroana_gocart ) then
do iv = 1, n_gocart_var
i_chem(iv)=i_tt+(iv-1)*lm+1
i_chem(iv)=i_tt+(iv-1)*lm+1 ! ???? (tt is 3D field, should i_chem be i_tt + iv*lm ?)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do you still need the comment here?

end do
endif

if ( wrf_pm2_5 ) then
iv=1
i_chem(iv)=i_tt+(iv-1)*lm+1
i_chem(iv)=i_tt+(iv-1)*lm+1 ! ???? (tt is 3D field, should i_chem be i_tt + iv*lm ?)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here

Copy link
Contributor

@RussTreadon-NOAA RussTreadon-NOAA left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did not test code. Most comments are minor.

Question: Do any of the existing GSI ctests exercise the functionality added by this PR? If not, (a) how do we know the changes work as intended? and (b) what's the strategy to ensure future PRs do not break the functionality added by this PR?

Comment on lines +1361 to +1362
! 2025-01-28 zhao - added code to read wave height (howv) & wind gust (gust)
! from intermediate binary file for the analysis in WRF-ARW based 3DRTMA
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Adding changelog entries to source code is no longer needed. Version control along with well documented issues record who changed the code and what the change is for.

It's fine to keep the new changelog entries added by this PR. I suggest that future PRs do not add changelogs to the source code.

@@ -62,6 +62,8 @@ subroutine convert_regional_guess_wrf(this,mype,ctph0,stph0,tlm0)
! 2005-05-24 pondeca - add 2dvar only surface analysis option
! 2005-07-06 parrish - add variable update_pint
! 2012-10-11 parrish - add byte_swap, which is set only on pe 0 and must be broadcast to all pes.
! 2025-01-29 zhao - add i_howv_3dda/i_gust_3dda/i_howv_mask, which are set only on pe 0 and
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment as before.

@@ -47,6 +47,8 @@ subroutine convert_netcdf_mass_wrf(this)
! - add code for hydrometer variables needed for direct reflectivity DA
! - add CV transform option on hydrometer variables for direct reflectivity DA
!
! 2025-01-27 zhao - added code for the input of wave height (howv) and wind gust (gust) fields from
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment as above.

! write(6,'(1x,A,A2,A,1x,I6)')trim(adjustl(myname_)),'::',trim(adjustl(derr_msg)),ierr
! call stop2(ierr)
call die(trim(adjustl(myname_)),trim(adjustl(derr_msg)),ierr)
i_howv_3dda = 0 ! skipping analysis of this variable
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Question: Doesn't die execute mpi_abort and terminate program execution? If so, execution will not reach the i_howv_3dda = 0 line, right?

Comment: The commented out line ! call stop2(ierr) can be removed.

! write(6,'(1x,A,A2,A,1x,I6)')trim(adjustl(myname_)),'::',trim(adjustl(derr_msg)),ierr
! call stop2(ierr)
call die(trim(adjustl(myname_)),trim(adjustl(derr_msg)),ierr)
i_gust_3dda = 0 ! skipping analysis of this variable
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similar to comment above.

@@ -2420,6 +2541,8 @@ subroutine update_netcdf_mass_wrf(this)
! 2020-09-13 CAPS(C. Liu, L. Chen, and H. Li)
! - add code for hydrometer variables needed for direct reflectivity DA
! - add CV transform option on hydrometer variables for direct reflectivity DA
! 2025-01-27 zhao - added code for the output of wave height (howv) and wind gust (gust) fields to
! the WRF-ARW (RAP/HRRR based) analysis file (netcdf fortmat).
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Source code changelog is not necessary.

lakemask,i_howv_mask)
!$$$ subprogram documentation block
! . . . .
! subprogram: unfill_mass_grid2t opposite of fill_mass_grid2
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This subroutine is unfill_mass_grid2t_ldlkmask, not unfill_mass_grid2t

! but screen off analysis increment over land and/or lake .
!
!
! program history log:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If this is a new subroutine added to GSI, why does it have a program history going back to 2004?

! Only add analysis increment over sea water (no land and lake)
select case (i_howv_mask)
case(2)
write(6,*) 'unfill_mass_grid2t_ldlkmask: using landlakemask to screen off HOWV analysis increment over land and lake.'
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do multiple mpi tasks call this routine? If so, write(6,*) could generate a lot of output. Do we need this output from all tasks? Should we place write statements in this routine inside verbose blocks?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
5 participants