From e117611b19181310d44f49aaf36d969e2aa6e68a Mon Sep 17 00:00:00 2001 From: Gang Zhao Date: Tue, 22 Aug 2023 13:29:42 +0000 Subject: [PATCH] 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) --- src/gsi/gsi_rfv3io_mod.f90 | 73 ++++++++++- src/gsi/gsimod.F90 | 15 ++- src/gsi/m_berror_stats_reg.f90 | 185 +++++++++++++++++++++++++-- src/gsi/rapidrefresh_cldsurf_mod.f90 | 53 +++++++- src/gsi/read_prepbufr.f90 | 8 ++ src/gsi/setuphowv.f90 | 12 ++ 6 files changed, 329 insertions(+), 17 deletions(-) diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 4fcb2aba1d..921cef7b99 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_in_anav, i_howv_in_data implicit none public type_fv3regfilenameg @@ -134,6 +137,7 @@ module gsi_rfv3io_mod 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_howv ! index for the significant wave height (aka howv in GSI) 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 @@ -145,6 +149,7 @@ module gsi_rfv3io_mod 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_howv ! index for the significant wave height (aka howv in GSI) parameter( & k_f10m =1, & !fact10 k_stype=2, & !soil_type @@ -159,7 +164,10 @@ 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 ) ! It might be better if n2d is a variable + ! depending on the variables in the list of + ! anavinfo, instead of a constant parameter. 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 +775,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 +826,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() ! --> howv real(r_kind),dimension(:,:,:),pointer::ges_ql=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qi=>NULL() @@ -1093,6 +1104,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 +1122,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 +1378,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_in_anav == 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 +1566,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) ! adding code to read howv if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then ! Convert 2m guess mixing ratio to specific humidity @@ -1782,7 +1802,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 +1812,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 @@ -1806,6 +1828,7 @@ 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: mpi_itype ! used to broadcast an integer 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 +1836,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 +1848,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 ! --> howv type (type_fv3regfilenameg),intent(in) :: fv3filenamegin character(len=max_varname_length) :: name integer(i_kind),allocatable,dimension(:):: dim @@ -1835,6 +1862,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, iret_howv + integer(i_kind) iret_bcast ! for io_layout > 1 real(r_kind),allocatable,dimension(:,:):: sfc_fulldomain @@ -1850,6 +1880,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 +1910,19 @@ 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 ( i_howv_in_anav == 1 ) then + iret_howv = nf90_inq_varid(gfile_loc,'howv',id_howv) + if ( iret_howv /= nf90_noerr ) then + i_howv_in_data = 0 ! howv does not exist in firstguess + write(6,'(1x,A,1x,A,1x,A,1x,I4,1x,A,1x,I4.4,A)') 'subroutine gsi_fv3ncdf2d_read:: howv is NOT found in firstguess ', trim(sfcdata), ', iret = ',iret_howv, ' (on pe: ', mype,')' + else + i_howv_in_data = 1 ! howv is found in firstguess + write(6,'(1x,A,1x,A,1x,A,1x,I4,1x,I4,1x,A,1x,I4.4,A)') 'subroutine gsi_fv3ncdf2d_read:: Found howv in firstguess ', trim(sfcdata), ', iret, varid = ',iret_howv, 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 +1950,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 ! howv (read if howv exists in firstguess, but no check on its existence.) + k=k_howv else cycle endif @@ -2036,6 +2084,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) if(allocated(sfc_fulldomain)) deallocate (sfc_fulldomain) endif ! mype +!-- broadcast i_howv_in_data to all tasks (!!!!) + call mpi_bcast(i_howv_in_data, 1, mpi_itype, mype_2d, mpi_comm_world, iret_bcast) !!!!!!! scatter !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call mpi_scatterv(work,ijns2d,displss2d,mpi_rtype,& @@ -2058,6 +2108,10 @@ 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_in_anav == 1 ) then + ges_howv(:,:)=sfcn2d(:,:,k_howv) !howv: even no howv in firstguess, sfcn2d(:,:,k_howv) + ! is still allocated and filled with initial values. + endif deallocate (sfcn2d,a) return end subroutine gsi_fv3ncdf2d_read @@ -3131,6 +3185,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: ! @@ -3173,6 +3229,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() ! --> howv integer(i_kind) i,k @@ -3289,6 +3346,10 @@ 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 +!-- howv + if ( i_howv_in_anav == 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 @@ -3498,6 +3559,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 only if howv is already in firstguess, i.e. i_howv_in_data = 1 + if ( i_howv_in_anav == 1 .and. i_howv_in_data == 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 cf885c2b64..1b6c252b04 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -175,7 +175,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, l_tuneBE_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 @@ -502,6 +503,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, l_tuneBE_howv +! (in namelist session rapidrefresh_cldsurf) ! !EOP !------------------------------------------------------------------------- @@ -1558,6 +1562,12 @@ 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) +! l_tuneBE_howv - logical, on/off the tuning of static BE of howv +! = .false. (default: turn off tuning the BE of howv) ! namelist/rapidrefresh_cldsurf/dfi_radar_latent_heat_time_period, & metar_impact_radius,metar_impact_radius_lowcloud, & @@ -1578,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, l_tuneBE_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..cb92bb4009 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -313,6 +313,8 @@ 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 + use rapidrefresh_cldsurf_mod, only: l_tuneBE_howv + use hybrid_ensemble_parameters, only: l_hyb_ens implicit none @@ -368,6 +370,8 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt ! vertical length scale from 1/layer to 1/sigma ! 2022-05-24 ESRL(H.Wang) - Add B for reginal FV3-CMAQ ! (berror_fv3_cmaq_regional=.true.) . +! 2023-07-30 Zhao - added code for tuning the static BE of howv in hybrid +! run (since ensemble of howv is not available yet). ! !EOP ___________________________________________________________________ @@ -415,6 +419,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt !real(r_kind), parameter :: corz_default=one,hwll_default=27000.00000000,vz_default=one*10 logical :: print_verbose real(r_kind) ,dimension(mlat,1,2) :: cov_dum + real(r_kind) :: b_howv ! character(256) :: filename ! filename = 'howv_var_berr.bin' @@ -870,16 +875,29 @@ 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) ! static B for howv 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 +!--- tuning the static BE of howv in hybrid run + if ( (.not. twodvar_regional) .and. l_hyb_ens ) then + b_howv = one + call get_factor_tuneBE_howv(mype, b_howv) + if ( l_tuneBE_howv ) then + corp(:,n) = corp(:,n) * b_howv + if (mype==0) then + write(6,'(1x,A,A,A,A,1x,F7.4)') myname_, ' Tuning factor of static BE corp( :,n) (for ', trim(adjustl(cvars2d(n))), ')= ', b_howv + write(6,*) myname_, ' Tuned static BE corp( :,n) (for ', trim(adjustl(cvars2d(n))), ')= ', corp(:,n) + end if + end if + end if ! corp(:,n)=cov_dum(:,1) !do i=1,mlat ! corp(i,n)=0.4_r_kind !#ww3 @@ -1055,7 +1073,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 +1108,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 +1123,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 @@ -1114,15 +1139,25 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) integer(i_kind), parameter :: lun34=34 character(len=256),dimension(npar) :: filename integer(i_kind), parameter :: dp1 = kind(1D0) + logical :: print_verbose +! + print_verbose=.false. + if(verbose)print_verbose=.true. ! 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 read in namelist + arrout(:,:,2) = hwllp_howv ! 17000.0_r_kind used in 3dvar (default) if not read 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,11 +1167,141 @@ 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 (print_verbose .and. mype .eq. 0) then + print*,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 (print_verbose .and. mype .eq. 0) then + print*,trim(adjustl(myname)), trim(filename(i_npar)) // ' does not exist for static BE of howv.' + print*,trim(adjustl(myname)), '::using pre-defined corp_howv=',corp_howv, ' hwllp_howv=',hwllp_howv, ' for howv.' + end if end if end do end subroutine read_howv_stats +!++++ +subroutine get_factor_tuneBE_howv(mype, b_factor) +!$$$ subprogram documentation block +! . . . . +! subprogram: get_factor_tuneBE_howv +! prgmmr: Zhao, Gang org: ncep/emc date: 2023-07-30 +! +! abstract: +! read in the weights of static BE(beta_s0, or beta_s) and ensemble +! BE (beta_e) used in hybrid envar, then compute the b_factor which +! is used to tune the static BE of howv. +! Because there is no ensemble of howv by far, so the contribution of +! ensemble to total BE of howv is zero, then due to the weight of +! static BE (beta_s) < 1.0, then the total BE of howv in hybrid run is +! smaller than the BE of howv used in pure 3dvar run. To make the +! analysis of howv in hybrid run is similar to the analysis of howv in +! pure 3dvar run, in hybrid run the staic BE of howv needs to be +! tuned (actually inflated since beta_s < 1.0). +! note: +! the name of file with ensemble weights and localisation +! information is 'hybens_info'. If it is changed to other name, user +! need to change it here, following the file name definition in +! subroutine hybens_localization_setup in hybrid_ensemble_isotropic.F90. +! +! +! program history log: +! 2023-07-30 Zhao - initialize the code. It adopts some lines of code from +! the subroutine hybens_localization_setup in +! hybrid_ensemble_isotropic.F90 to read in the beta_s +! and beta_e from file "hybens_info". +! +! input argument list: +! mype - pe number of this task run +! output argument list: +! b_factor - the factor used to tune the static BE of howv +! +! attributes: +! language: f90 +! machine: theia/gyre +! +!$$$ end documentation block +! + use hybrid_ensemble_parameters, only: readin_beta, readin_localization + use hybrid_ensemble_parameters, only: beta_s0,beta_e0 + use rapidrefresh_cldsurf_mod, only: l_tuneBE_howv + use gsi_io, only : verbose +! + implicit none +! Declare passed variables + integer(i_kind), intent(in ) :: mype ! "my" processor ID + real(r_kind), intent(inout) :: b_factor +! Declare local variables + integer(i_kind),parameter :: lunin = 49 + character(len=40),parameter :: fname = 'hybens_info' + integer(i_kind) :: istat, msig + real(r_kind) :: beta + real(r_kind) :: b_s, b_e, s_hv, s_vv + + logical :: lexist + logical :: print_verbose +!---------------------------------------------------- + print_verbose=.false. + if(verbose)print_verbose=.true. + + if ( readin_localization .or. readin_beta ) then ! read info from file + + inquire(file=trim(fname),exist=lexist) + if ( lexist ) then + open(lunin,file=trim(fname),form='formatted') + rewind(lunin) + read(lunin,100,iostat=istat) msig + if ( istat /= 0 ) then + l_tuneBE_howv = .false. + beta = one + if (mype == 0) then + write(6,'(1x,A,A,A,I4)') 'get_factor_tuneBE_howv: error reading file ',trim(fname),' iostat = ',istat + write(6,'(1x,A)') 'get_factor_tuneBE_howv: skipping reading beta_s and No tuning static BE of howv' + end if + else + read(lunin,101) s_hv, s_vv, b_s, b_e ! read the value for lowest level only + beta = b_s/(b_s + b_e) ! in case b_s + b_e /= 1 + if(mype==0 .and. print_verbose) then + write(6,*) "get_factor_tuneBE_howv: read in BETA_S, BETA_E at level 1: ", b_s, b_e + write(6,*) "get_factor_tuneBE_howv: used beta= ", beta + end if + end if + else + l_tuneBE_howv = .false. + beta = one + if (mype == 0) then + write(6,'(1x,A,A)') 'get_factor_tuneBE_howv: could NOT find file ',trim(fname) + write(6,'(1x,A)') 'get_factor_tuneBE_howv: skipping reading beta_s and No tuning static BE of howv' + end if + end if + + else + + beta = beta_s0 + if(mype==0 .and. print_verbose) then + write(6,*) "get_factor_tuneBE_howv: the universal BETA_S0 : ", beta_s0, ' (no localisation file)' + write(6,*) "get_factor_tuneBE_howv: used beta= ", beta + end if + + end if + + close(lunin) + + if ( beta >= 0.001_r_kind ) then + b_factor = sqrt(one/beta) + else + b_factor = one ! actually in this case, b_factor should be a huge value + ! since it is using pure ensemble in hybrid envar run + l_tuneBE_howv = .false. ! for pure ensemble run, no tuning of static BE of howv + end if + +100 format(I4) +!101 format(F8.1,3x,F5.1,2(3x,F8.4)) +101 format(F8.1,3x,F8.3,F8.4,3x,F8.4) + +300 continue + + return + +end subroutine get_factor_tuneBE_howv end module m_berror_stats_reg diff --git a/src/gsi/rapidrefresh_cldsurf_mod.f90 b/src/gsi/rapidrefresh_cldsurf_mod.f90 index 1ee35fffba..608e160541 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -28,7 +28,13 @@ 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 +! l_tuneBE_howv: on/off the tuning of static BE of howv (in hyrid run) +! i_howv_in_anav: check if howv is in anavinfo +! i_howv_in_data: check if howv is in firstguess file ! ! Subroutines Included: ! sub init_rapidrefresh_cldsurf - initialize RR related variables to default values @@ -181,6 +187,24 @@ 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_in_anav - integer, check if howv is found in anavinfo as state/analysis variable +! = 0 if "howv" is not in anavinfo as state/analysis variable +! = 1 if "howv" is found in anavinfo as state/analysis variable +! i_howv_in_data - integer, check if "howv" is found in firstguess data file +! = 0 if "howv" is not in firstguess data file +! = 1 if "howv" is found in firstgues file +! l_tunBE_howv - namelist logical, on/off the tuning of static BE of howv ONLY in hybrid run +! = .true., (default) tuning the static BE of howv background error of howv +! note: in hybrid envar run, the static BE is redueced by beta_s (<1.0), +! since there is no howv ensemble, there is no +! ensemble contribution to the total BE of howv, +! then the total BE of howv is just the reduced +! static BE of howv. 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 +276,9 @@ module rapidrefresh_cldsurf_mod public :: l_saturate_bkCloud public :: l_rtma3d public :: i_precip_vertical_check + public :: corp_howv, hwllp_howv + public :: i_howv_in_anav, i_howv_in_data + public :: l_tuneBE_howv logical l_hydrometeor_bkio real(r_kind) dfi_radar_latent_heat_time_period @@ -310,6 +337,9 @@ 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_in_anav, i_howv_in_data + logical :: l_tuneBE_howv contains @@ -325,6 +355,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 +369,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 +454,21 @@ 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) + l_tuneBE_howv = .false. ! no tuning the static BE of howv in hybrid run (default) + i_howv_in_data = 0 ! no significant wave height (howv) in firstguess data file (default) + i_howv_in_anav = 0 ! no howv in anavinfo (default) + +!-- searching for speficif 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_in_anav = 1 + if (mype == 0 ) then + write(6,'(1x,A,1x,A8,1x,A)')"init_rapidrefresh_cldsurf: anavinfo svars2d (state variable): ",trim(adjustl(svars2d(i2))), " is found in anavinfo." + 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 355441e209..d2f71480c1 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) +! in prepbufr data file (i_howv_in_anav=1 and i_howv_in_data=1) ! input argument list: ! infile - unit from which to read BUFR data @@ -222,6 +224,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: i_howv_in_anav, i_howv_in_data use gsi_io, only: verbose use phil2, only: denest ! hilbert curve @@ -1132,6 +1135,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 (significant wave height, aka howv, +! and howv must be listed in anavinfo and in firstguess) + if ( i_howv_in_anav == 1 .and. i_howv_in_data == 1 ) 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 diff --git a/src/gsi/setuphowv.f90 b/src/gsi/setuphowv.f90 index c2b1dfe3e9..c0e784e281 100644 --- a/src/gsi/setuphowv.f90 +++ b/src/gsi/setuphowv.f90 @@ -33,6 +33,9 @@ subroutine setuphowv(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag ! . Remove my_node with corrected typecast(). ! 2018-01-08 pondeca - addd option l_closeobs to use closest obs to analysis ! time in analysis +! 2023-07-30 Zhao - using the option i_howv_in_data to control the usage of +! obs of howv in setuphowv only if honw is available +! in figrstguess (i_howv_in_data=1). (not applied for 2DRTMA) ! ! input argument list: ! lunin - unit from which to read observations @@ -84,6 +87,7 @@ subroutine setuphowv(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag use gsi_bundlemod, only : gsi_bundlegetpointer use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle use rapidrefresh_cldsurf_mod, only: l_closeobs + use rapidrefresh_cldsurf_mod, only: i_howv_in_data implicit none ! Declare passed variables @@ -324,6 +328,14 @@ subroutine setuphowv(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag muse(i) = .true. endif +!-- Note: +! The following if-block must NOT be applied for 2DRTMA (ie. twodvar_regional=.false.), +! because 2DRTMA does NOT check the existence of howv in firstguess, +! so no defnition of i_howv_in_data in 2DRTMA. +! In the analysis of non-2DRTMA, if there is no howv data in firstguess, +! (i_howv_in_data /= 1), then do not use the obs of howv. (set muse = .false.) + if ( (.not. twodvar_regional) .and. (i_howv_in_data /= 1) ) muse(i) = .false. + if (nobskeep>0 .and. luse_obsdiag) call obsdiagNode_get(my_diag, jiter=nobskeep,muse=muse(i)) ! Compute penalty terms (linear & nonlinear qc).