Skip to content

Commit

Permalink
Adding code to analyze the siginificant wave heigh in GSI 3D Analysis…
Browse files Browse the repository at this point in the history
… (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)

<!-- PLEASE READ -->
<!--
Before opening a PR, please note these guidelines:

- Each PR should only address ONE topic and have an associated issue
- No hardcoded or paths to personal directories should be present
- No temporary or backup files should be committed
- Any code that was disabled by being commented out should be removed
-->

**Description**

<!-- Please include relevant motivation and context. -->
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.
<!-- Please include a summary of the change and which issue is fixed.
-->
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.


<!-- List any dependencies that are required for this change. -->
No dependencies are required for this change.
<!-- Please provide reference to the issue this pull request is
addressing. -->
This PR is addressing the issue
[#601](#601): Adding code to
analyze the siginificant wave heigh in GSI 3D Analysis".
<!-- For e.g. Fixes #IssueNumber -->
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?**

<!-- Please describe the tests that you ran to verify your changes and
on the platforms these tests were conducted. -->
- 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.

<!-- Provide instructions so we can reproduce. -->
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)
<!-- Please also list any relevant details for your test configuration
-->
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.
  • Loading branch information
GangZhao-NOAA committed Sep 29, 2023
1 parent ca19008 commit 728d006
Show file tree
Hide file tree
Showing 5 changed files with 160 additions and 20 deletions.
78 changes: 71 additions & 7 deletions src/gsi/gsi_rfv3io_mod.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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()
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -1805,23 +1823,28 @@ 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
use general_commvars_mod, only: ltosi_s,ltosj_s
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

integer(i_kind),intent(in) :: it
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
Expand All @@ -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
Expand All @@ -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))

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,&
Expand All @@ -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
Expand Down Expand Up @@ -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:
!
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down
13 changes: 11 additions & 2 deletions src/gsi/gsimod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
!-------------------------------------------------------------------------
Expand Down Expand Up @@ -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, &
Expand All @@ -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
Expand Down
41 changes: 31 additions & 10 deletions src/gsi/m_berror_stats_reg.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
Loading

0 comments on commit 728d006

Please sign in to comment.