From 29d9d8fdea7e5139bca6bdd135711a33e7bf25a6 Mon Sep 17 00:00:00 2001 From: Gang Zhao <53267411+GangZhao-NOAA@users.noreply.github.com> Date: Fri, 5 Apr 2024 09:39:54 -0400 Subject: [PATCH] Adding I/O for direct analysis of near-surface wind gust for RRFS-based 3DRTMA (#730) **Description** To improve the analysis of the near-surface wind gust in 3DRTMA, the observations of near-surface wind gust would be analyzed directly in GSI (3DVar and Hybrid 3DEnVar), instead of being a derived product from the near-surface wind analysis. Since the core subroutines for direct variational assimilation of wind gust (e.g., setupgust.f90, intgust.f90, stpgust.f90, etc.) had already been implemented in GSI for 2DRTMA, so in the work the development in GSI mainly focus on adding I/O of 2-D wind gust firstguess and analysis fields for RRFS-based 3DRTMA, and some minor modifications in observation and background error for wind gust, options to control the analysis of wind gust, etc. This PR is to address the issue #726 : Adding I/O for direct analysis of near-surface wind gust for RRFS-based 3DRTMA Fixes #726 **Type of change** Please delete options that are not relevant. - [x] New feature (non-breaking change which adds functionality) **How Has This Been Tested?** **Checklist** - [x] My code follows the style guidelines of this project - [x] I have performed a self-review of my own code - [x] I have commented my code, particularly in hard-to-understand areas - [x] New and existing tests pass with my changes tested with a real case - 2024-02-20_12:00Z, 1. dry-run (using my updated GSI code with wind gust analysis, but actually no wind gust obs is analyzed, so-called dryrun) is compared to control-run (original GSI code running without wind gust obs): the results are identical. This indicates that if without analyzing wind-gust obs, then the updated code generates the analysis identical to the analysis of original/control code. Or say, the added code does not have influence on the other part of code. 2. real case run with updated GSI code to analyze the obs of wind gust: The following figure shows the used observations of near-surface wind gust: ![var_obs_2024022012_gust_used_maprll_datll_reg_ncf](https://github.com/NOAA-EMC/GSI/assets/53267411/ecbe479a-03c6-490f-a179-9e0027291468) the following figure shows the analysis increments: ![GUST_hyb_hwllp90_corptuned_inc_incrintrp_maprll_datrll_reg_grb2](https://github.com/NOAA-EMC/GSI/assets/53267411/a01fca0d-dc1f-438b-b8eb-e624de35a631) - [x] Any dependent changes have been merged and published - [x] Regression tests on WCOSS2 (Cactus) and Hera (Rocky-8) : my updated GSI commit [#f91f247d](https://github.com/GangZhao-NOAA/GSI/commit/f91f247dde65dea0b659c0b78d747433d9ca9559)) vs control/original GSI code (commit [#6d9ebbb7](https://github.com/NOAA-EMC/GSI/commit/6d9ebbb7896b92a93959ce63c7a1ad9e9a0aab4f)) Here is the reports of the regression tests on WCOSS2 (Cactus): ~~~ [gang.zhao@clogin02:build] (feature/windgust_in_3dvar_for_3drtma)$ ctest -j 7 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_glbens Start 4: netcdf_fv3_regional Start 5: hafs_4denvar_glbens Start 6: hafs_3denvar_hybens Start 7: global_enkf 1/7 Test #4: netcdf_fv3_regional .............. Passed 483.15 sec 2/7 Test #3: rrfs_3denvar_glbens .............. Passed 486.74 sec 3/7 Test #7: global_enkf ...................... Passed 850.98 sec 4/7 Test #2: rtma ............................. Passed 970.28 sec 5/7 Test #6: hafs_3denvar_hybens .............. Passed 1152.82 sec 6/7 Test #5: hafs_4denvar_glbens .............. Passed 1213.93 sec 7/7 Test #1: global_4denvar ................... Passed 1683.16 sec 100% tests passed, 0 tests failed out of 7 Total Test time (real) = 1683.19 sec ~~~ Here is the reports of the regression tests on Hera (Rocky8): ~~~ (base) [Gang.Zhao@hfe11:build] (feature/windgust_in_3dvar_for_3drtma)$ ctest -j 7 Test project /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build Start 1: global_4denvar Start 2: rtma Start 3: rrfs_3denvar_glbens Start 4: netcdf_fv3_regional Start 5: hafs_4denvar_glbens Start 6: hafs_3denvar_hybens Start 7: global_enkf 1/7 Test #4: netcdf_fv3_regional .............. Passed 491.53 sec 2/7 Test #3: rrfs_3denvar_glbens ..............***Failed 495.27 sec 3/7 Test #2: rtma ............................. Passed 982.45 sec 4/7 Test #6: hafs_3denvar_hybens .............. Passed 1168.99 sec 5/7 Test #7: global_enkf ...................... Passed 1239.77 sec 6/7 Test #5: hafs_4denvar_glbens ..............***Failed 1347.87 sec 7/7 Test #1: global_4denvar ................... Passed 1974.45 sec 71% tests passed, 2 tests failed out of 7 Total Test time (real) = 1974.91 sec The following tests FAILED: 3 - rrfs_3denvar_glbens (Failed) 5 - hafs_4denvar_glbens (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. (base) [Gang.Zhao@hfe11:build] (feature/windgust_in_3dvar_for_3drtma)$ ctest -R rrfs_3denvar_glbens Test project /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build Start 3: rrfs_3denvar_glbens 1/1 Test #3: rrfs_3denvar_glbens .............. Passed 430.52 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 430.55 sec (base) [Gang.Zhao@hfe11:build] (feature/windgust_in_3dvar_for_3drtma)$ ctest -R hafs_4denvar_glbens Test project /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build Start 5: hafs_4denvar_glbens 1/1 Test #5: hafs_4denvar_glbens .............. Passed 1225.37 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 1225.39 sec ~~~ **Note**: _When I was running the regression tests, GSI code was just updated to the latest commit [#b53740a7](https://github.com/GangZhao-NOAA/GSI/commit/f91f247dde65dea0b659c0b78d747433d9ca9559). Considering the frequent update in EMC GSI code recently and saving the time, after this PR has been reviewed and approved by peer-reviewers, I will update the code to latest EMC GSI commit, then re-run the regression tests for final approval. --- src/gsi/gsi_rfv3io_mod.f90 | 65 +++++++++++++++++++++++----- src/gsi/gsimod.F90 | 7 ++- src/gsi/m_berror_stats_reg.f90 | 29 ++++++++++++- src/gsi/rapidrefresh_cldsurf_mod.f90 | 50 +++++++++++++++++++++ src/gsi/read_prepbufr.f90 | 5 ++- 5 files changed, 140 insertions(+), 16 deletions(-) diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 652bad9a33..8e1c3ab98f 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -61,7 +61,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 + use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda implicit none public type_fv3regfilenameg @@ -147,7 +147,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,k_howv + public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv,k_gust 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 @@ -158,7 +158,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,k_howv + integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv,k_gust parameter( & k_f10m =1, & !fact10 k_stype=2, & !soil_type @@ -174,7 +174,8 @@ module gsi_rfv3io_mod k_q2m =12, & ! 2 m Q k_orog =13, & !terrain k_howv =14, & ! significant wave height (aka howv in GSI) - n2d=14 ) + k_gust =15, & ! wind gust (aka gust in GSI) + n2d=15 ) logical :: grid_reverse_flag character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_dynmetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields @@ -996,6 +997,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) 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_gust=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_ql=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qi=>NULL() @@ -1274,6 +1276,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) else if(trim(vartem)=='t2m') then else if(trim(vartem)=='q2m') then else if(trim(vartem)=='howv') then + else if(trim(vartem)=='gust') then else write(6,*)'the metvarname2 ',trim(vartem),' has not been considered yet, stop' call stop2(333) @@ -1294,7 +1297,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) 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") & - .or.(trim(vartem)=="howv"))) then ! z is treated separately + .or.(trim(vartem)=="howv").or.(trim(vartem)=="gust"))) then ! z is treated separately if (ifindstrloc(vardynvars,trim(vartem)) > 0) then jdynvar=jdynvar+1 fv3lam_io_dynmetvars2d_nouv(jdynvar)=trim(vartem) @@ -1557,6 +1560,12 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if (ier/=0) call die(trim(myname),'cannot get pointers for howv, ier=',ier) endif +!--- wind gust (gust) + if ( i_gust_3dda == 1 ) then + call GSI_BundleGetPointer(GSI_MetGuess_Bundle(it),'gust',ges_gust,istatus ); ier=ier+istatus + if (ier/=0) call die(trim(myname),'cannot get pointers for gust, 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)) @@ -1739,7 +1748,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif - call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m,ges_howv) + call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m, & + ges_howv,ges_gust) if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then ! Convert 2m guess mixing ratio to specific humidity @@ -1975,7 +1985,8 @@ 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,ges_howv) +subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m, & + ges_howv,ges_gust) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf2d_read @@ -2022,6 +2033,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) 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 + real(r_kind), intent(in),dimension(:,:),pointer::ges_gust type (type_fv3regfilenameg),intent(in) :: fv3filenamegin character(len=max_varname_length) :: name @@ -2036,8 +2048,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) 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 +! for checking the existence of howv/gust in firstguess file + integer(i_kind) id_howv, id_gust integer(i_kind) iret_bcast ! for io_layout > 1 @@ -2057,8 +2069,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) allocate(work(itotsub*n2d)) allocate( sfcn2d(lat2,lon2,n2d)) -!-- initialisation of the array for howv +!-- initialisation of the array for howv/gust sfcn2d(:,:,k_howv) = zero + sfcn2d(:,:,k_gust) = zero !-- initialisation of the array for sfc_var_exist sfc_var_exist = .false. @@ -2107,6 +2120,21 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) trim(sfcdata), ', iret, varid = ',iret, id_howv,' (on pe: ', mype,').' end if end if +!--- check the existence of wind gust (gust) in 2D FV3-LAM firstguess file +! (similar as done above for howv) + if ( i_gust_3dda == 1 ) then + iret = nf90_inq_varid(gfile_loc,'gust',id_gust) + if ( iret /= nf90_noerr ) then + iret = nf90_inq_varid(gfile_loc,'GUST',id_gust) ! double check with name in uppercase + end if + if ( iret /= nf90_noerr ) then + i_gust_3dda = 0 ! gust does not exist in firstguess, then stop GSI run. + call die('gsi_fv3ncdf2d_read','Warning: CANNOT find gust 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 gust in firstguess ', & + trim(sfcdata), ', iret, varid = ',iret, id_gust,' (on pe: ', mype,').' + end if + end if !!!!!!!!!!!! read in 2d variables !!!!!!!!!!!!!!!!!!!!!!!!!! do i=ndimensions+1,nvariables @@ -2150,6 +2178,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) else if( trim(name)=='HOWV'.or.trim(name)=='howv' ) then k=k_howv sfc_var_exist(k) = .true. + else if( trim(name)=='GUST'.or.trim(name)=='gust' ) then + k=k_gust + sfc_var_exist(k) = .true. else cycle endif @@ -2283,8 +2314,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) if(allocated(sfc_fulldomain)) deallocate (sfc_fulldomain) endif ! mype -!-- broadcast the updated i_howv_3dda to all tasks (!!!!) +!-- broadcast the updated i_howv_3dda, i_gust_3dda to all tasks (!!!!) call mpi_bcast(i_howv_3dda, 1, mpi_itype, mype_2d, mpi_comm_world, iret_bcast) + call mpi_bcast(i_gust_3dda, 1, mpi_itype, mype_2d, mpi_comm_world, iret_bcast) !-- broadcast the updated sfc_var_exist to all tasks (!!!!) call mpi_bcast(sfc_var_exist, n2d, mpi_itype, mype_2d, mpi_comm_world, iret_bcast) @@ -2313,6 +2345,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) if ( i_howv_3dda == 1 ) then if ( sfc_var_exist(k_howv) ) ges_howv(:,:)=sfcn2d(:,:,k_howv) endif + if ( i_gust_3dda == 1 ) then + if ( sfc_var_exist(k_gust) ) ges_gust(:,:)=sfcn2d(:,:,k_gust) + endif deallocate (sfcn2d,a) return end subroutine gsi_fv3ncdf2d_read @@ -3628,6 +3663,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) real(r_kind),pointer,dimension(:,: ):: ges_t2m =>NULL() real(r_kind),pointer,dimension(:,: ):: ges_q2m =>NULL() real(r_kind),pointer,dimension(:,: ):: ges_howv =>NULL() + real(r_kind),pointer,dimension(:,: ):: ges_gust =>NULL() integer(i_kind) i,k @@ -3750,6 +3786,9 @@ subroutine wrfv3_netcdf(fv3filenamegin) if ( i_howv_3dda == 1 ) then call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'howv',ges_howv,istatus); ier=ier+istatus endif + if ( i_gust_3dda == 1 ) then + call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'gust',ges_gust,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 @@ -3964,6 +4003,10 @@ subroutine wrfv3_netcdf(fv3filenamegin) if ( i_howv_3dda == 1 ) then call gsi_fv3ncdf_write_sfc(fv3filenamegin,'howv',ges_howv,add_saved) endif +!-- output analysis of gust + if ( i_gust_3dda == 1 ) then + call gsi_fv3ncdf_write_sfc(fv3filenamegin,'gust',ges_gust,add_saved) + endif if(allocated(g_prsi)) deallocate(g_prsi) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 8a1ce896bb..d7f5667252 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -184,7 +184,7 @@ module gsimod 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, & - corp_howv, hwllp_howv + corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust 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 @@ -1604,6 +1604,9 @@ module gsimod ! = 0.42 meters (default) ! hwllp_howv - real, background error de-correlation length scale of howv ! = 170,000.0 meters (default 170 km) +! corp_gust - real, static background error of gust (stddev error) +! hwllp_gust - real, background error de-correlation length scale of gust +! oerr_gust - real, observation error of gust ! namelist/rapidrefresh_cldsurf/dfi_radar_latent_heat_time_period, & metar_impact_radius,metar_impact_radius_lowcloud, & @@ -1625,7 +1628,7 @@ module gsimod 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, & - corp_howv, hwllp_howv + corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust ! 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 8730e56c3b..d7a30808e6 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -313,6 +313,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt use mpeu_util,only: getindex use radiance_mod, only: icloud_cv,n_clouds_fwd,cloud_names_fwd use chemmod, only: berror_fv3_cmaq_regional,berror_fv3_sd_regional + use rapidrefresh_cldsurf_mod, only: corp_gust, hwllp_gust, l_rtma3d implicit none @@ -825,11 +826,35 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt ! end if else if (n==nrf2_gust) then do i=1,mlat - corp(i,n)=three + corp(i,n)=three ! background error stddev of wind gust = 3 m/s (default: legacy code from 2DRTMA) end do do i=0,mlat+1 - hwllp(i,n)=hwll(i,1,nrf3_q) + hwllp(i,n)=hwll(i,1,nrf3_q) ! de-correlation length of bkgd error of gust is + ! same as the value of q at bottom level (default: legacy code from 2DRTMA) + ! for other DA apps, it is recommended to change it + ! by setting hwllp_gust in GSI namelist. end do + if ( l_rtma3d ) then ! For 3drtma only: allowing to change the stddev and + ! de-correlation length of bkgd error of gust: + ! corp_gust : set in namelist(if <=0, using default value above (3.0) + ! hwllp_gust: set in namelist(if <=0, using default value above (value of q) + if ( corp_gust .gt. 0.0_r_kind ) then + corp(1:mlat, n) = corp_gust + if (mype==0) write(6,'(1x,A,A,I5.5,A,F8.3)') & + myname_,"@pe=",mype," (3drtma) set b_error stddev of gust = ",corp_gust + else + if (mype==0) write(6,'(1x,A,A,I5.5,A,F8.3)') & + myname_,"@pe=",mype," (3drtma) set b_error stddev of gust (default) = ",three + end if + if ( hwllp_gust .gt. 0.0_r_kind ) then + hwllp(0:mlat+1,n) = hwllp_gust + if (mype==0) write(6,'(1x,A,A,I5.5,A,F12.3)') & + myname_,"@pe=",mype," (3drtma) set b_error de-corr length of gust = ",hwllp_gust + else + if (mype==0) write(6,'(1x,A,A,I5.5,A)') & + myname_,"@pe=",mype," (3drtma) set b_error de-corr length of gust is same as length of q." + end if + end if else if (n==nrf2_vis) then do i=1,mlat corp(i,n)=3.0_r_kind diff --git a/src/gsi/rapidrefresh_cldsurf_mod.f90 b/src/gsi/rapidrefresh_cldsurf_mod.f90 index 122d2872d0..475f44a9d3 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -197,6 +197,30 @@ module rapidrefresh_cldsurf_mod ! 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). +! corp_gust - namelist real, static BE of gust (standard error deviation) +! note: 1. initialised to be an arbitary negative value, in order to skip this +! negative value, instead to use value (3.0 m/s) set in subroutine +! berror_read_wgt_reg as default. +! 2. (3drtma only) if a user-specified value (e.g., 2.0 m/s) is preferred +! for corp_gust, in GSI namelist session "rapidrefresh_cldsurf", +! set "corp_gust=2.0," +! hwllp_gust - namelist real, static BE de-correlation length scale of gust +! note: 1. initialised to be an arbitary negative value, in order to skip this +! negative value, instead to use value (same value for q) set in +! subroutine berror_read_wgt_reg as default +! 2. (3drtma only) if a user-specified value (e.g., 100 km) is preferred +! for hwllp_gust, in GSI namelist session "rapidrefresh_cldsurf", +! set "hwllp_gust=100000.0," +! oerr_gust - namelist real, observation error of gust +! note: 1. initialised to be an arbitary negative value, in order to skip this +! negative value, instead to use value (1.0 m/s) set in read_prepbufr.f90 +! 2. (3drtma only) if a user-specified value (e.g., 1.5 m/s ) is preferred +! for oerr_gust, in GSI namelist session "rapidrefresh_cldsurf", +! set "oerr_gust=1.5," +! i_gust_3dda - integer, control the analysis of gust in 3D analysis (either var or hybrid) +! = 0 (gust-off: default) : no analysis of gust in 3D analysis. +! = 1 (gust-on) : if variable name "gust" is found in anavinfo, +! set it to be 1 to turn on analysis of gust; ! ! attributes: ! language: f90 @@ -270,6 +294,8 @@ module rapidrefresh_cldsurf_mod public :: i_precip_vertical_check public :: corp_howv, hwllp_howv public :: i_howv_3dda + public :: corp_gust, hwllp_gust, oerr_gust + public :: i_gust_3dda logical l_hydrometeor_bkio real(r_kind) dfi_radar_latent_heat_time_period @@ -330,6 +356,8 @@ module rapidrefresh_cldsurf_mod integer(i_kind) i_precip_vertical_check real(r_kind) :: corp_howv, hwllp_howv integer(i_kind) :: i_howv_3dda + real(r_kind) :: corp_gust, hwllp_gust, oerr_gust + integer(i_kind) :: i_gust_3dda contains @@ -447,6 +475,22 @@ subroutine init_rapidrefresh_cldsurf 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) + corp_gust = -1.50_r_kind ! initialised as negative & void to be skipped, in order to use + ! the value (3.0 m/s) set in sub berror_read_wgt_reg (as default). + ! If user-specified value is preferred, set it in session + ! "rapidrefresh_cldsurf" of GSI namelist file + + hwllp_gust = -90000.0_r_kind ! initialised as a value, in order to skip this negative value + ! and to use the value (used for q) set in sub berror_read_wgt_reg. + ! If user-specified value is preferred, set it in session + ! "rapidrefresh_cldsurf" of GSI namelist file + + oerr_gust = -2.5_r_kind ! initialised as a negative value, in order to skip this negative value + ! and to use the value (1.0 m/s) set in read_prepbufr.f90 + ! If user-specified value is preferred, set it in session + ! "rapidrefresh_cldsurf" of GSI namelist file + + i_gust_3dda = 0 ! no analysis of wind gust (gust) in 3D analysis (default) !-- searching for specific variable in state variable list (reading from anavinfo) do i2=1,ns2d @@ -456,6 +500,12 @@ subroutine init_rapidrefresh_cldsurf 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 + if ( trim(svars2d(i2))=='gust' .or. trim(svars2d(i2))=='GUST' ) then + i_gust_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_gust_3dda = ", i_gust_3dda + end if + end if end do ! i2 : looping over 2-D anasv return diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index f281573a4e..bed3b31db2 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -227,7 +227,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& use adjust_cloudobs_mod, only: adjust_convcldobs,adjust_goescldobs use mpimod, only: npe use rapidrefresh_cldsurf_mod, only: i_gsdsfc_uselist,i_gsdqc,i_ens_mean - use rapidrefresh_cldsurf_mod, only: l_rtma3d + use rapidrefresh_cldsurf_mod, only: l_rtma3d, oerr_gust use gsi_io, only: verbose use phil2, only: denest ! hilbert curve @@ -1825,6 +1825,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& gustqm=0 if (kx==188 .or. kx==288 .or. kx==195 .or. kx==295 ) & call get_gustqm(kx,c_station_id,c_prvstg,c_sprvstg,gustqm) + if ( l_rtma3d ) gustqm = 0 ! skipping get_gustqm for 3drtma run (missing list file) qm=gustqm else if(visob) then visqm=0 ! need to fix this later @@ -2556,6 +2557,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! need to find out gustoe ! gustoe=1.8 gustoe=1.0 + if ( l_rtma3d .and. oerr_gust > 0.0_r_kind ) gustoe = oerr_gust selev=stnelev oelev=obsdat(4,k) if(selev == oelev)oelev=r10+selev @@ -2577,6 +2579,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if ((kx==188).or.(kx==288) .or.(kx==195) .or.(kx==295)) then ! gustoe=2.5 gustoe=1.0 + if ( l_rtma3d .and. oerr_gust > 0.0_r_kind ) gustoe = oerr_gust windcorr=abs(obsdat(5,k))<1.0 .and. abs(obsdat(6,k))<1.0 .and. obsdat(8,k)>10.0 if (windcorr) gustoe=gustoe*1.5_r_kind