From 728d006c36ae6b17f08b380ff23b3422b97624fe Mon Sep 17 00:00:00 2001 From: Gang <53267411+GangZhao-NOAA@users.noreply.github.com> Date: Fri, 29 Sep 2023 09:11:39 -0400 Subject: [PATCH] Adding code to analyze the siginificant wave heigh in GSI 3D Analysis (see issue #601) (#614) Adding code to analyze the siginificant wave heigh in GSI 3D Analysis, esp. for FV3-LAM model based DA, eg. RRFS-DA, RRFS-3DRTMA. (Also see the issue in EMC GSI github repository: #601 Adding I/O for Analysis of Significant Wave Height for 3DRTMA) **Description** Significant Wave Height (hereafter as SWH) is one of the standard products provided by the operational (2D)RTMA. To continuously provide the same products in 3DRTMA, the next-generation RTMA, some efforts in GSI code need to be made in order to analyze the SWH in 3D analysis of GSI. The kernel subroutines to assimilate SWH in GSI (such as stphowv.f90, setuphowv.f90, inthowv.f90, gsi_howvOper.f90 and m_howvNode.f90) already had been added for (2D)RTMA years ago by Manuel Pondeca, so for this issue, the code work mainly focus on adding the I/O of SWH in background and analysis fields for 3DRTMA (esp. RRFS-based 3DRTMA), and some necessary modifications in background error, options, variables related to analysis of SWH, etc. Modified code in GSI: 1. rapidrefresh_cldsurf_mod.f90: adding a few variables related to the analysis of howv in 3D analysis 2. gsimod.F90: adding namelist options used for analysis of howv in 3D analysis 3. m_berror_stats_reg.f90: added some code for the special treatment to the static background error (BE) of howv 4. read_prepbufr.f90: adding code to decode the observation of howv in prepbufr file when howv is available in firstguess 5. setuphowv.f90: adding code to use obs of howv when howv is available in firstguess 6. gsi_rfv3io_mod.f90: adding I/O code to read in howv from firstguess and write out howv into analysis. No dependencies are required for this change. This PR is addressing the issue [#601](https://github.com/NOAA-EMC/GSI/issues/601): Adding code to analyze the siginificant wave heigh in GSI 3D Analysis". Fixes #601 **Type of change** Please delete options that are not relevant. - [*] New feature (non-breaking change which adds functionality) **How Has This Been Tested?** - Brief results from ctest (regression test) with the modified code (on WCOSS2 - Cactus): [gang.zhao@clogin07:build] (feature/3drtma_howv)$ ctest -N Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Test #1: global_3dvar Test #2: global_4dvar Test #3: global_4denvar Test #4: hwrf_nmm_d2 Test #5: hwrf_nmm_d3 Test #6: rtma Test #7: rrfs_3denvar_glbens Test #8: netcdf_fv3_regional Test #9: global_enkf Total Tests: 9 Test #1: global_3dvar [gang.zhao@clogin04:build] (feature/3drtma_howv)$ ctest -R global_3dvar Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 1: global_3dvar 1/1 Test #1: global_3dvar ..................... Passed 1631.12 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 1631.14 sec Test #2: global_4dvar [gang.zhao@clogin09:build] (feature/3drtma_howv)$ ctest -R global_4dvar Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 2: global_4dvar 1/1 Test #2: global_4dvar ..................... Passed 2462.19 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 2462.23 sec Test #3: global_4denvar [gang.zhao@clogin04:build] (feature/3drtma_howv)$ ctest -R global_4denvar Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 3: global_4denvar 1/1 Test #3: global_4denvar ................... Passed 1922.43 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 1922.46 sec Test #4: hwrf_nmm_d2 [gang.zhao@clogin09:build] (feature/3drtma_howv)$ ctest -R hwrf_nmm_d2 Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 4: hwrf_nmm_d2 1/1 Test #4: hwrf_nmm_d2 ...................... Passed 1214.10 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 1214.20 sec Test #5: hwrf_nmm_d3 [gang.zhao@clogin09:build] (feature/3drtma_howv)$ ctest -R hwrf_nmm_d3 Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 5: hwrf_nmm_d3 1/1 Test #5: hwrf_nmm_d3 ...................... Passed 736.38 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 736.50 sec Test #6: rtma [gang.zhao@clogin05:build] (feature/3drtma_howv)$ ctest -R rtma Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 6: rtma 1/1 Test #6: rtma ............................. Passed 1027.01 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 1027.01 sec Test #7: rrfs_3denvar_glbens [gang.zhao@clogin06:build] (feature/3drtma_howv)$ ctest -R rrfs_3denvar_glbens Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 7: rrfs_3denvar_glbens 1/1 Test #7: rrfs_3denvar_glbens .............. Passed 484.69 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 484.70 sec Test #8: netcdf_fv3_regional [gang.zhao@clogin03:build] (feature/3drtma_howv)$ ctest -R netcdf_fv3_regional Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 8: netcdf_fv3_regional 1/1 Test #8: netcdf_fv3_regional .............. Passed 483.08 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 483.11 sec Test #9: global_enkf [gang.zhao@clogin03:build] (feature/3drtma_howv)$ ctest -R global_enkf Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/WaveHgt/develop/build Start 9: global_enkf 1/1 Test #9: global_enkf ...................... Passed 488.50 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 488.57 sec - The modified GSI code passed the regression tests (all 9 tasks) on Hera and WCOSS2 (Cactus). - adding the analysis of howv only has very trial influences on the analyses of other variables. Here is the statistics of the differences of other variables (u/v/t/ps/q/t2m/q2m) from the runs of GSI without howv vs. with howv (from a testing case 2023-07-12_14:00:00 UTC in 3km North-American domain): comparing two netcdf files: fcst_fv3lam_hyb_betas/INPUT/fv_core.res.tile1.nc fcst_fv3lam_nodata_noinfo/INPUT/fv_core.res.tile1.nc ... Variable Group Count Sum AbsSum Min Max Range Mean StdDev u / 602135550 3926.84 25760.8 -0.1026 0.485788 0.588388 6.52152e-06 0.00115817 v / 620166777 -4891.34 32582.5 -0.835774 0.268402 1.10418 -7.88714e-06 0.00197793 T / 155987083 178.048 6497.51 -0.0246582 0.0384064 0.0630646 1.14143e-06 0.000218737 delp / 19559676 -281.532 3008.29 -0.00292969 0.00219727 0.00512695 -1.43935e-05 0.000183727 comparing two netcdf files: fcst_fv3lam_hyb_betas/INPUT/fv_tracer.res.tile1.nc fcst_fv3lam_nodata_noinfo/INPUT/fv_tracer.res.tile1.nc ... Variable Group Count Sum AbsSum Min Max Range Mean StdDev sphum / 430707614 0.594287 2.77816 -2.6139e-05 3.1759e-05 5.7898e-05 1.37979e-09 8.03072e-08 comparing two netcdf files: fcst_fv3lam_hyb_betas/INPUT/sfc_data.nc fcst_fv3lam_nodata_noinfo/INPUT/sfc_data.nc ... Variable Group Count Sum AbsSum Min Max Range Mean StdDev t2m / 10665000 43.3899 135.095 -0.00152825 0.00686629 0.00839454 4.06844e-06 5.02866e-05 q2m / 10665000 0.0192553 0.124707 -3.1476e-06 1.77554e-05 2.0903e-05 1.80547e-09 5.89657e-08 It could be seen that the differences are trivial and ignorable. The regression tests were done by following the instructions of "[GSI Ctests (regression tests)](https://github.com/NOAA-EMC/GSI/wiki/GSI-Ctests-(regression-tests))" in [GSI Wiki](https://github.com/NOAA-EMC/GSI/wiki) The modified code had also been tested with a testing case 2023-07-12_14:00:00 UTC for 3km North-American domain Here is a brief summary of the test results: 1. Here is the analysis increment of Significant Wave Height (aka howv hereafter): pure 3dvar, static background error of howv is 0.42 meters, and the de-correlation length scale is 170km. ![HOWV_var_inc_maprll_datll_reg_ncf](https://github.com/NOAA-EMC/GSI/assets/53267411/4fdeeb82-7258-4344-be69-cce747474312) 2. The following figure shows the distribution of howv in the analysis (used obs is in green, rejected in red). Obviously the location of used obs of howv match the area of non-zero analysis increments of howv. ![var_obs_2023071214_howv_maprll_datll_reg_ncf](https://github.com/NOAA-EMC/GSI/assets/53267411/d4ed6013-cfc8-486e-8f47-db07ec0e4e53) 3. The following figure is the analysis increment of howv with hybrid envar analysis (using gdas ensemble 80 members and the ensemble weight is 84%), and the static BE of howv is tuned/inflated. The analysis increments are very similar to the results from pure 3dvar run (see the first figure) ![HOWV_hyb_betas016_inc_maprll_datll_reg_ncf](https://github.com/NOAA-EMC/GSI/assets/53267411/e6e696e8-932b-42ab-9001-3472e970b21c) 4. The last figure shows the analysis increments of howv with hybrid envar analysis (using gdas ensemble 80 members and the ensemble weight is 84%), but the static BE of howv is NOT tuned. It can be observed that the analysis increments is less than the results from the hybrid run with tuning the static BE of howv. That is because the weight of static BE (16%) reduced the background error of howv (ensemble of howv is not available yet), so the impact of obs is decreased. ![HOWV_hyb_betas016_noTune_inc_maprll_datll_reg_ncf](https://github.com/NOAA-EMC/GSI/assets/53267411/ca25d068-fc86-4d47-a9d2-46e02ac22dac) **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 **DUE DATE for this PR is 10/5/2023.** If this PR is not merged into `develop` by this date, the PR will be closed and returned to the developer. --- src/gsi/gsi_rfv3io_mod.f90 | 78 +++++++++++++++++++++++++--- src/gsi/gsimod.F90 | 13 ++++- src/gsi/m_berror_stats_reg.f90 | 41 +++++++++++---- src/gsi/rapidrefresh_cldsurf_mod.f90 | 41 ++++++++++++++- src/gsi/read_prepbufr.f90 | 7 +++ 5 files changed, 160 insertions(+), 20 deletions(-) diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 62b23ee713..6d16be7c13 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -22,6 +22,8 @@ module gsi_rfv3io_mod ! used as background for surface observation operator ! 2022-04-15 Wang - add IO for regional FV3-CMAQ (RRFS-CMAQ) model ! 2022-08-10 Wang - add IO for regional FV3-SMOKE (RRFS-SMOKE) model +! 2023-07-30 Zhao - add IO for the analysis of the significant wave height +! (SWH, aka howv in GSI) in fv3-lam based DA (eg., RRFS-3DRTMA) ! ! subroutines included: ! sub gsi_rfv3io_get_grid_specs @@ -56,6 +58,7 @@ module gsi_rfv3io_mod use rapidrefresh_cldsurf_mod, only: i_use_2mq4b,i_use_2mt4b use chemmod, only: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3,laeroana_fv3cmaq use chemmod, only: naero_smoke_fv3,aeronames_smoke_fv3,laeroana_fv3smoke + use rapidrefresh_cldsurf_mod, only: i_howv_3dda implicit none public type_fv3regfilenameg @@ -133,7 +136,7 @@ module gsi_rfv3io_mod public :: mype_u,mype_v,mype_t,mype_q,mype_p,mype_oz,mype_ql public :: mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w public :: k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc - public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m + public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv public :: ijns,ijns2d,displss,displss2d,ijnz,displsz_g public :: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv public :: fv3lam_io_tracerchemvars3d_nouv,fv3lam_io_tracersmokevars3d_nouv @@ -144,7 +147,7 @@ module gsi_rfv3io_mod integer(i_kind) mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w integer(i_kind) k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc - integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m + integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv parameter( & k_f10m =1, & !fact10 k_stype=2, & !soil_type @@ -159,7 +162,8 @@ module gsi_rfv3io_mod k_t2m =11, & ! 2 m T k_q2m =12, & ! 2 m Q k_orog =13, & !terrain - n2d=13 ) + k_howv =14, & ! significant wave height (aka howv in GSI) + n2d=14 ) logical :: grid_reverse_flag character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_dynmetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields @@ -767,6 +771,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) ! 2022-04-01 Y. Wang and X. Wang - add capability to read reflectivity ! for direct radar EnVar DA using reflectivity as state ! variable, poc: xuguang.wang@ou.edu +! 2023-07-30 Zhao - added code to read significant wave height (howv) field +! from the 2D fv3-lam firstguess file (fv3_sfcdata). ! attributes: ! language: f90 ! machine: ibm RS/6000 SP @@ -816,6 +822,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) real(r_kind),pointer,dimension(:,:,:):: ges_delp =>NULL() real(r_kind),dimension(:,:),pointer::ges_t2m=>NULL() real(r_kind),dimension(:,:),pointer::ges_q2m=>NULL() + real(r_kind),dimension(:,:),pointer::ges_howv=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_ql=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qi=>NULL() @@ -1093,6 +1100,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if(mype == 0) write(6,*)'the metvarname ',trim(vartem),' will be dealt separately' else if(trim(vartem)=='t2m') then else if(trim(vartem)=='q2m') then + else if(trim(vartem)=='howv') then else write(6,*)'the metvarname2 ',trim(vartem),' has not been considered yet, stop' call stop2(333) @@ -1110,7 +1118,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) do i=1,size(name_metvars2d) vartem=trim(name_metvars2d(i)) if(.not.( (trim(vartem)=='ps'.and.fv3sar_bg_opt==0).or.(trim(vartem)=="z") & - .or.(trim(vartem)=="t2m").or.(trim(vartem)=="q2m"))) then !z is treated separately + .or.(trim(vartem)=="t2m").or.(trim(vartem)=="q2m") & + .or.(trim(vartem)=="howv"))) then ! z is treated separately if (ifindstrloc(vardynvars,trim(vartem)) > 0) then jdynvar=jdynvar+1 fv3lam_io_dynmetvars2d_nouv(jdynvar)=trim(vartem) @@ -1365,6 +1374,13 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it),'t2m',ges_t2m,istatus );ier=ier+istatus if (ier/=0) call die(trim(myname),'cannot get pointers for t2m,ier=',ier) endif + +!--- significant wave height (howv) + if ( i_howv_3dda == 1 ) then + call GSI_BundleGetPointer(GSI_MetGuess_Bundle(it),'howv',ges_howv,istatus ); ier=ier+istatus + if (ier/=0) call die(trim(myname),'cannot get pointers for howv, ier=',ier) + endif + if(mype == 0 ) then call check(nf90_open(fv3filenamegin(it)%dynvars,nf90_nowrite,loc_id)) call check(nf90_inquire(loc_id,formatNum=ncfmt)) @@ -1546,7 +1562,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif - call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m) + call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m,ges_howv) if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then ! Convert 2m guess mixing ratio to specific humidity @@ -1782,7 +1798,7 @@ end subroutine gsi_bundlegetpointer_fv3lam_tracerchem_nouv end subroutine read_fv3_netcdf_guess -subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) +subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf2d_read @@ -1792,6 +1808,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) ! Scatter the field to each PE ! program history log: ! 2023-02-14 Hu - Bug fix for read in subdomain surface restart files +! 2023-07-30 Zhao - added IO to read significant wave height (howv) from 2D FV3-LAM +! firstguess file (fv3_sfcdata) ! ! input argument list: ! it - time index for 2d fields @@ -1805,7 +1823,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) ! !$$$ end documentation block use kinds, only: r_kind,i_kind - use mpimod, only: ierror,mpi_comm_world,npe,mpi_rtype,mype + use mpimod, only: ierror,mpi_comm_world,npe,mpi_rtype,mype,mpi_itype + use mpeu_util, only: die use guess_grids, only: fact10,soil_type,veg_frac,veg_type,sfc_rough, & sfct,sno,soil_temp,soil_moi,isli use gridmod, only: lat2,lon2,itotsub,ijn_s @@ -1813,8 +1832,11 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable + use netcdf, only: nf90_inq_varid + use netcdf, only: nf90_noerr use mod_fv3_lola, only: fv3_h_to_ll,nxa,nya use constants, only: grav + use constants, only: zero implicit none @@ -1822,6 +1844,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) real(r_kind),intent(in),dimension(:,:),pointer::ges_z real(r_kind),intent(in),dimension(:,:),pointer::ges_t2m real(r_kind),intent(in),dimension(:,:),pointer::ges_q2m + real(r_kind),intent(in),dimension(:,:),pointer::ges_howv type (type_fv3regfilenameg),intent(in) :: fv3filenamegin character(len=max_varname_length) :: name integer(i_kind),allocatable,dimension(:):: dim @@ -1835,6 +1858,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) integer(i_kind) kk,n,ns,j,ii,jj,mm1 character(len=:),allocatable :: sfcdata !='fv3_sfcdata' character(len=:),allocatable :: dynvars !='fv3_dynvars' +! for checking the existence of howv in firstguess file + integer(i_kind) id_howv + integer(i_kind) iret_bcast ! for io_layout > 1 real(r_kind),allocatable,dimension(:,:):: sfc_fulldomain @@ -1850,6 +1876,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) allocate(work(itotsub*n2d)) allocate( sfcn2d(lat2,lon2,n2d)) +!-- initialisation of the array for howv + sfcn2d(:,:,k_howv) = zero + if(mype==mype_2d ) then allocate(sfc_fulldomain(nx,ny)) @@ -1877,6 +1906,24 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) iret=nf90_inquire_dimension(gfile_loc,k,name,len) dim(k)=len enddo + +!--- check the existence of significant wave height (howv) in 2D FV3-LAM firstguess file +! if howv is set in anavinfo (as i_howv_3dda=1), then check its existence in firstguess, +! but if it is not found in firstguess, then stop GSI run and set i_howv_3dda = 0. + if ( i_howv_3dda == 1 ) then + iret = nf90_inq_varid(gfile_loc,'howv',id_howv) + if ( iret /= nf90_noerr ) then + iret = nf90_inq_varid(gfile_loc,'HOWV',id_howv) ! double check with name in uppercase + end if + if ( iret /= nf90_noerr ) then + i_howv_3dda = 0 ! howv does not exist in firstguess, then stop GSI run. + call die('gsi_fv3ncdf2d_read','Warning: CANNOT find howv in firstguess, aborting..., iret = ', iret) + else + write(6,'(1x,A,1x,A,1x,A,1x,I4,1x,I4,1x,A,1x,I4.4,A)') 'gsi_fv3ncdf2d_read:: Found howv in firstguess ', & + trim(sfcdata), ', iret, varid = ',iret, id_howv,' (on pe: ', mype,').' + end if + end if + !!!!!!!!!!!! read in 2d variables !!!!!!!!!!!!!!!!!!!!!!!!!! do i=ndimensions+1,nvariables iret=nf90_inquire_variable(gfile_loc,i,name,len) @@ -1904,6 +1951,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) k=k_t2m else if( trim(name)=='Q2M'.or.trim(name)=='q2m' ) then k=k_q2m + else if( trim(name)=='HOWV'.or.trim(name)=='howv' ) then + k=k_howv else cycle endif @@ -2036,6 +2085,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) if(allocated(sfc_fulldomain)) deallocate (sfc_fulldomain) endif ! mype +!-- broadcast the updated i_howv_3dda to all tasks (!!!!) + call mpi_bcast(i_howv_3dda, 1, mpi_itype, mype_2d, mpi_comm_world, iret_bcast) !!!!!!! scatter !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call mpi_scatterv(work,ijns2d,displss2d,mpi_rtype,& @@ -2058,6 +2109,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) ges_t2m(:,:)=sfcn2d(:,:,k_t2m) ges_q2m(:,:)=sfcn2d(:,:,k_q2m) endif + if ( i_howv_3dda == 1 ) then + ges_howv(:,:)=sfcn2d(:,:,k_howv) + endif deallocate (sfcn2d,a) return end subroutine gsi_fv3ncdf2d_read @@ -3192,6 +3246,8 @@ subroutine wrfv3_netcdf(fv3filenamegin) ! 2019-11-22 CAPS(C. Tong) - modify "add_saved" to properly output analyses ! 2021-01-05 x.zhang/lei - add code for updating delz analysis in regional da ! 2022-04-01 Y. Wang and X. Wang - add code for updating reflectivity +! 2023-07-30 Zhao - added code for the output of the analysis of +! significant wave height (howv) ! ! input argument list: ! @@ -3234,6 +3290,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) real(r_kind),pointer,dimension(:,:,:):: ges_q =>NULL() real(r_kind),pointer,dimension(:,: ):: ges_t2m =>NULL() real(r_kind),pointer,dimension(:,: ):: ges_q2m =>NULL() + real(r_kind),pointer,dimension(:,: ):: ges_howv =>NULL() integer(i_kind) i,k @@ -3350,6 +3407,9 @@ subroutine wrfv3_netcdf(fv3filenamegin) call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'q2m',ges_q2m,istatus); ier=ier+istatus call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'t2m',ges_t2m,istatus );ier=ier+istatus endif + if ( i_howv_3dda == 1 ) then + call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'howv',ges_howv,istatus); ier=ier+istatus + endif if (ier/=0) call die('wrfv3_netcdf','cannot get pointers for fv3 met-fields, ier =',ier) if (laeroana_fv3cmaq) then @@ -3559,6 +3619,10 @@ subroutine wrfv3_netcdf(fv3filenamegin) call gsi_fv3ncdf_write_sfc(fv3filenamegin,'t2m',ges_t2m,add_saved) call gsi_fv3ncdf_write_sfc(fv3filenamegin,'q2m',ges_q2m,add_saved) endif +!-- output analysis of howv + if ( i_howv_3dda == 1 ) then + call gsi_fv3ncdf_write_sfc(fv3filenamegin,'howv',ges_howv,add_saved) + endif if(allocated(g_prsi)) deallocate(g_prsi) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 2656a2dce4..70618120d0 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -176,7 +176,8 @@ module gsimod i_coastline,i_gsdqc,qv_max_inc,ioption,l_precip_clear_only,l_fog_off,& cld_bld_coverage,cld_clr_coverage,& i_cloud_q_innovation,i_ens_mean,DTsTmax,& - i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check + i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check, & + corp_howv, hwllp_howv use gsi_metguess_mod, only: gsi_metguess_init,gsi_metguess_final use gsi_chemguess_mod, only: gsi_chemguess_init,gsi_chemguess_final use tcv_mod, only: init_tcps_errvals,tcp_refps,tcp_width,tcp_ermin,tcp_ermax @@ -503,6 +504,9 @@ module gsimod ! 3. fv3_cmaq_regional = .true. ! 4. berror_fv3_cmaq_regional = .true. ! 09-15-2022 yokota - add scale/variable/time-dependent localization +! 2023-07-30 Zhao - added namelist options for analysis of significant wave height +! (aka howv in GSI code): corp_howv, hwllp_howv +! (in namelist session rapidrefresh_cldsurf) ! !EOP !------------------------------------------------------------------------- @@ -1560,6 +1564,10 @@ module gsimod ! = 2(clean Qg as in 1, and adjustment to the retrieved Qr/Qs/Qnr throughout the whole profile) ! = 3(similar to 2, but adjustment to Qr/Qs/Qnr only below maximum reflectivity level ! and where the dbz_obs is missing); +! corp_howv - real, static background error of howv (stddev error) +! = 0.42 meters (default) +! hwllp_howv - real, background error de-correlation length scale of howv +! = 170,000.0 meters (default 170 km) ! namelist/rapidrefresh_cldsurf/dfi_radar_latent_heat_time_period, & metar_impact_radius,metar_impact_radius_lowcloud, & @@ -1580,7 +1588,8 @@ module gsimod i_coastline,i_gsdqc,qv_max_inc,ioption,l_precip_clear_only,l_fog_off,& cld_bld_coverage,cld_clr_coverage,& i_cloud_q_innovation,i_ens_mean,DTsTmax, & - i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check + i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check, & + corp_howv, hwllp_howv ! chem(options for gsi chem analysis) : ! berror_chem - .true. when background for chemical species that require diff --git a/src/gsi/m_berror_stats_reg.f90 b/src/gsi/m_berror_stats_reg.f90 index 2ff8a6aa94..601339e1ac 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -870,16 +870,17 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt hwllp(i,n)=hwllp(i,nrf2_ps) end do else if (n==nrf2_howv) then - call read_howv_stats(mlat,1,2,cov_dum) + call read_howv_stats(mlat,1,2,cov_dum,mype) do i=1,mlat corp(i,n)=cov_dum(i,1,1) !#ww3 hwllp(i,n) = cov_dum(i,1,2) end do hwllp(0,n) = hwllp(1,n) hwllp(mlat+1,n) = hwllp(mlat,n) - - if (mype==0) print*, 'corp(i,n) = ', corp(:,n) - if (mype==0) print*, ' hwllp(i,n) = ', hwllp(:,n) + if (mype==0) then + print*, myname_, ' static BE corp( :,n) (for ', trim(adjustl(cvars2d(n))), ')= ', corp(:,n) + print*, myname_, ' static BE hwllp(:,n) (for ', trim(adjustl(cvars2d(n))), ')= ', hwllp(:,n) + end if ! corp(:,n)=cov_dum(:,1) !do i=1,mlat ! corp(i,n)=0.4_r_kind !#ww3 @@ -1055,7 +1056,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt end subroutine berror_read_wgt_reg !++++ -subroutine read_howv_stats(nlat,nlon,npar,arrout) +subroutine read_howv_stats(nlat,nlon,npar,arrout,mype) !$$$ subprogram documentation block ! . . . . ! subprogram: read_howv_stats @@ -1090,6 +1091,9 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) ! program history log: ! 2016-08-03 stelios ! 2016-08-26 stelios : Compatible with GSI. +! 2023-07-30 Zhao - added code to set the background error +! standard deviation (corp_howv) and de-correlation +! length scale (hwllp_howv) for non-2DRTMA run ! input argument list: ! filename - The name of the file ! output argument list: @@ -1102,10 +1106,14 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) !$$$ end documentation block ! use kinds,only : r_kind, i_kind + use gridmod, only : twodvar_regional + use rapidrefresh_cldsurf_mod, only : corp_howv, hwllp_howv + use gsi_io, only : verbose ! implicit none ! Declare passed variables integer(i_kind), intent(in )::nlat,nlon,npar + integer(i_kind), intent(in ) :: mype ! "my" processor ID real(r_kind), dimension(nlat ,nlon, npar), intent( out)::arrout ! Declare local variables integer(i_kind) :: reclength,i,j,i_npar @@ -1117,12 +1125,18 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) ! filename(1) = 'howv_var_berr.bin' filename(2) = 'howv_lng_berr.bin' -! - arrout(:,:,1)=0.42_r_kind - arrout(:,:,2)=50000.0_r_kind +!-- first, assign the pre-defined values to corp and hwllp + if ( twodvar_regional ) then + arrout(:,:,1)=0.42_r_kind ! values were specified by Manuel and Stelio for 2DRTMA + arrout(:,:,2)=50000.0_r_kind ! values were specified by Manuel and Stelio for 2DRTMA + else + arrout(:,:,1) = corp_howv ! 0.42_r_kind used in 3dvar (default) if not set in namelist + arrout(:,:,2) = hwllp_howv ! 17000.0_r_kind used in 3dvar (default) if not set in namelist + end if reclength=nlat*r_kind -! +!-- secondly, if files for corp and hwllp are available, then read them in for +! corp and hwllp. If the files are not found, then use the pre-defined values. do i_npar = 1,npar inquire(file=trim(filename(i_npar)), exist=file_exists) if (file_exists)then @@ -1132,9 +1146,16 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) read(unit=lun34 ,rec=j) (arrout(i,j,i_npar), i=1,nlat) enddo close(unit=lun34) + if (verbose .and. mype .eq. 0) then + write(6,'(1x,A,1x,A2,1x,A)') trim(adjustl(myname)), '::', & + trim(filename(i_npar))//' is used for background error of howv.' + end if else - print*,myname, trim(filename(i_npar)) // ' does not exist' + if (verbose .and. mype .eq. 0) then + write(6,'(1x,A,1x,A2,1x,A)') trim(adjustl(myname)), '::', & + trim(filename(i_npar))//' does not exist for static BE of howv, using pre-defined values.' + end if end if end do end subroutine read_howv_stats diff --git a/src/gsi/rapidrefresh_cldsurf_mod.f90 b/src/gsi/rapidrefresh_cldsurf_mod.f90 index 1ee35fffba..122d2872d0 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -28,7 +28,11 @@ module rapidrefresh_cldsurf_mod ! option for checking and adjusting the profile of Qr/Qs/Qg/Qnr ! retrieved through cloud analysis to reduce the background ! reflectivity ghost in analysis. (default is 0) -! +! 2023-07-30 Zhao added options for analysis of significant wave height +! (SWH, aka howv in GSI code): +! corp_howv: to set the static background error of howv +! hwllp_howv: to set the de-correlation length scale +! i_howv_3dda: control the analysis of howv in 3D analysis (if howv is in anavinfo) ! ! Subroutines Included: ! sub init_rapidrefresh_cldsurf - initialize RR related variables to default values @@ -181,6 +185,18 @@ module rapidrefresh_cldsurf_mod ! = 2(clean Qg as in 1, and adjustment to the retrieved Qr/Qs/Qnr throughout the whole profile) ! = 3(similar to 2, but adjustment to Qr/Qs/Qnr only below maximum reflectivity level ! and where the dbz_obs is missing); +! corp_howv - namelist real, static BE of howv (standard error deviation) +! hwllp_howv - namelist real, static BE de-correlation length scale of howv +! i_howv_3dda - integer, control the analysis of howv in 3D analysis (either var or hybrid) +! = 0 (howv-off: default) : no analysis of howv in 3D analysis. +! = 1 (howv-on) : if variable name "howv" is found in anavinfo, +! set it to be 1 to turn on analysis of howv; +! note: in hybrid envar run, the static BE is redueced by beta_s (<1.0), +! since there is no ensemble of howv currently yet, then no ensemble +! 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). ! ! attributes: ! language: f90 @@ -252,6 +268,8 @@ module rapidrefresh_cldsurf_mod public :: l_saturate_bkCloud public :: l_rtma3d public :: i_precip_vertical_check + public :: corp_howv, hwllp_howv + public :: i_howv_3dda logical l_hydrometeor_bkio real(r_kind) dfi_radar_latent_heat_time_period @@ -310,6 +328,8 @@ module rapidrefresh_cldsurf_mod logical l_saturate_bkCloud logical l_rtma3d integer(i_kind) i_precip_vertical_check + real(r_kind) :: corp_howv, hwllp_howv + integer(i_kind) :: i_howv_3dda contains @@ -325,6 +345,8 @@ subroutine init_rapidrefresh_cldsurf ! 2008-06-03 Hu initial build for cloud analysis ! 2010-03-29 Hu change names to init_rapidrefresh_cldsurf ! 2011--5-04 Todling inquire MetGuess for presence of hyrometeors & set default +! 2023-07-30 Zhao added code for initialization of some variables used +! in analysis of significant wave height ! ! input argument list: ! @@ -337,8 +359,12 @@ subroutine init_rapidrefresh_cldsurf !$$$ use kinds, only: i_kind use gsi_metguess_mod, only: gsi_metguess_get + use mpimod, only: mype + use state_vectors, only: ns2d,svars2d + implicit none integer(i_kind) ivar,i,ier + integer(i_kind) i2 logical have_hmeteor(5) character(len=2),parameter :: hydrometeors(5) = (/ 'qi', & 'ql', & @@ -418,6 +444,19 @@ subroutine init_rapidrefresh_cldsurf l_saturate_bkCloud= .true. l_rtma3d = .false. ! turn configuration for rtma3d off i_precip_vertical_check = 0 ! No check and adjustment to retrieved Qr/Qs/Qg (default) + corp_howv = 0.42_r_kind ! 0.42 meters (default) + hwllp_howv = 170000.0_r_kind ! 170,000.0 meters (170km as default for 3DRTMA, 50km is used in 2DRTMA) + i_howv_3dda = 0 ! no analysis of significant wave height (howv) in 3D analysis (default) + +!-- searching for specific variable in state variable list (reading from anavinfo) + do i2=1,ns2d + if ( trim(svars2d(i2))=='howv' .or. trim(svars2d(i2))=='HOWV' ) then + i_howv_3dda = 1 + if ( mype == 0 ) then + write(6,'(1x,A,1x,A8,1x,A,1x,I4)')"init_rapidrefresh_cldsurf: anavinfo svars2d (state variable): ",trim(adjustl(svars2d(i2))), " is found in anavinfo, set i_howv_3dda = ", i_howv_3dda + end if + end if + end do ! i2 : looping over 2-D anasv return end subroutine init_rapidrefresh_cldsurf diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index 9efd06418c..304fa62590 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -149,6 +149,8 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! 2020-05-04 wu - no rotate_wind for fv3_regional ! 2020-09-05 CAPS(C. Tong) - add flag for new vadwind obs to assimilate around the analysis time only ! 2023-03-23 draper - add code for processing T2m and q2m for global system +! 2023-07-30 Zhao - added code to extract obs of significant wave height (howvob) from bufr record +! in prepbufr file for 3D analysis ! input argument list: ! infile - unit from which to read BUFR data @@ -1132,6 +1134,11 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if (howvob) call ufbint(lunin,owave,1,255,levs,owavestr) if (cldchob) call ufbint(lunin,cldceilh,1,255,levs,cldceilhstr) endif +! Extract obs of howv in 3D Analysis +! (if-block is to avoid potential issue if decoding the bufr record twice in 2DRTMA run) + if ( .not. twodvar_regional ) then + if (howvob) call ufbint(lunin,owave,1,255,levs,owavestr) + endif if(kx==224 .and. newvad) then call ufbint(lunin,fcstdat,3,255,levs,'UFC VFC TFC ') end if