From e117611b19181310d44f49aaf36d969e2aa6e68a Mon Sep 17 00:00:00 2001 From: Gang Zhao Date: Tue, 22 Aug 2023 13:29:42 +0000 Subject: [PATCH 1/7] 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). From f60375621fd6f43521d0359f0223bacb92a01dda Mon Sep 17 00:00:00 2001 From: Gang Zhao Date: Wed, 13 Sep 2023 03:02:07 +0000 Subject: [PATCH 2/7] some modifications to the code of analysis of significant wave height (howv) in 3D analysis, based on feedbacks from code reviewers. --- src/gsi/rapidrefresh_cldsurf_mod.f90 | 50 +++++++++++----------------- src/gsi/read_prepbufr.f90 | 11 +++--- src/gsi/setuphowv.f90 | 12 ------- 3 files changed, 24 insertions(+), 49 deletions(-) diff --git a/src/gsi/rapidrefresh_cldsurf_mod.f90 b/src/gsi/rapidrefresh_cldsurf_mod.f90 index 608e160541..736e9af341 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -28,13 +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): +! 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 +! i_howv_3dda: control the analysis of howv in 3D analysis ! ! Subroutines Included: ! sub init_rapidrefresh_cldsurf - initialize RR related variables to default values @@ -189,22 +187,16 @@ module rapidrefresh_cldsurf_mod ! 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). +! i_howv_in_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 @@ -277,8 +269,7 @@ module rapidrefresh_cldsurf_mod 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 + public :: i_howv_3dda logical l_hydrometeor_bkio real(r_kind) dfi_radar_latent_heat_time_period @@ -338,8 +329,7 @@ module rapidrefresh_cldsurf_mod 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 + integer(i_kind) :: i_howv_3dda contains @@ -456,16 +446,14 @@ subroutine init_rapidrefresh_cldsurf 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) + i_howv_3dda = 0 ! no analysis of significant wave height (howv) in 3D analysis (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 ( trim(svars2d(i2))=='howv' ) then + i_howv_3dda = 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." + 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 diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index ceeb0f150d..171682b803 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -149,8 +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) +! 2023-07-30 Zhao - added code to extract obs of significant wave height (howvob) for 3D analysis +! in prepbufr data file ! input argument list: ! infile - unit from which to read BUFR data @@ -224,7 +224,6 @@ 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 @@ -1135,9 +1134,9 @@ 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 +! Extract obs of howv in 3D Analysis +! (if-block is to avoid potential issue if reading same dataset 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 diff --git a/src/gsi/setuphowv.f90 b/src/gsi/setuphowv.f90 index c0e784e281..c2b1dfe3e9 100644 --- a/src/gsi/setuphowv.f90 +++ b/src/gsi/setuphowv.f90 @@ -33,9 +33,6 @@ 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 @@ -87,7 +84,6 @@ 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 @@ -328,14 +324,6 @@ 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). From 8432590ddd63a41a6405806a3f4c41b9cfa12832 Mon Sep 17 00:00:00 2001 From: Gang Zhao Date: Thu, 14 Sep 2023 16:09:58 +0000 Subject: [PATCH 3/7] based on feedback from code reviewers, modified the code of analysis for significant wave height: 1) combine two variables i_howv_in_anav and i_howv_in_data into one variable i_howv_3dda 2) if howv is required in analyis (as set howv in anavinfo), but if howv is not available in firstguess, then stop GSI run. 3) removing the code which tunes the static BE of howv in hybrid envar run. So in hybrid run, user needs to set the corp_howv in namelist if desired. 4) cleaning some unnecessary comments --- src/gsi/gsi_rfv3io_mod.f90 | 62 +++++----- src/gsi/gsimain.f90 | 1 + src/gsi/gsimod.F90 | 8 +- src/gsi/m_berror_stats_reg.f90 | 163 ++------------------------- src/gsi/rapidrefresh_cldsurf_mod.f90 | 6 +- src/gsi/read_prepbufr.f90 | 6 +- 6 files changed, 51 insertions(+), 195 deletions(-) diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index b5bde24adc..89388b7bb9 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -58,7 +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 + use rapidrefresh_cldsurf_mod, only: i_howv_3dda implicit none public type_fv3regfilenameg @@ -136,8 +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_howv ! index for the significant wave height (aka howv in GSI) + 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 @@ -148,8 +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_howv ! index for the significant wave height (aka howv in GSI) + 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 @@ -165,9 +163,7 @@ 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 ) ! It might be better if n2d is a variable - ! depending on the variables in the list of - ! anavinfo, instead of a constant parameter. + 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 @@ -826,7 +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() ! --> howv + real(r_kind),dimension(:,:),pointer::ges_howv=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_ql=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qi=>NULL() @@ -1104,7 +1100,6 @@ 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) @@ -1380,7 +1375,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif !--- significant wave height (howv) - if ( i_howv_in_anav == 1 ) then + 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 @@ -1566,7 +1561,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif - call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m,ges_howv) ! adding code to read howv + 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 @@ -1827,8 +1822,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) ! !$$$ 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 mpimod, only: ierror,mpi_comm_world,npe,mpi_rtype,mype,mpi_itype 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 @@ -1848,7 +1842,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) 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 + 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 @@ -1912,14 +1906,22 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) enddo !--- check the existence of significant wave height (howv) in 2D FV3-LAM firstguess file - if ( i_howv_in_anav == 1 ) then +! 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_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,')' + iret_howv = nf90_inq_varid(gfile_loc,'HOWV',id_howv) ! double check with name in uppercase + end if + if ( iret_howv /= nf90_noerr ) then + i_howv_3dda = 0 ! howv does not exist in firstguess, then stop GSI run. + 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,'). Stop running GSI!!!!' + call stop2(345) 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,')' + i_howv_3dda = 1 ! howv does exist in firstguess, running analysis with howv + 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 @@ -1950,7 +1952,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) 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.) + else if( trim(name)=='HOWV'.or.trim(name)=='howv' ) then k=k_howv else cycle @@ -2084,8 +2086,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) 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) +!-- 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,& @@ -2108,9 +2110,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) 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. + if ( i_howv_3dda == 1 ) then + ges_howv(:,:)=sfcn2d(:,:,k_howv) endif deallocate (sfcn2d,a) return @@ -3290,7 +3291,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 + real(r_kind),pointer,dimension(:,: ):: ges_howv =>NULL() integer(i_kind) i,k @@ -3407,8 +3408,7 @@ 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 + 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) @@ -3620,8 +3620,8 @@ 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 +!-- output analysis of howv + if ( i_howv_3dda == 1 ) then call gsi_fv3ncdf_write_sfc(fv3filenamegin,'howv',ges_howv,add_saved) endif diff --git a/src/gsi/gsimain.f90 b/src/gsi/gsimain.f90 index 5eaa38286a..7be4e3a646 100644 --- a/src/gsi/gsimain.f90 +++ b/src/gsi/gsimain.f90 @@ -530,6 +530,7 @@ program gsi ! = 342 - setuplight: failure to allocate obsdiags ! = 343 - setuplight: failure to allocate obsdiags ! = 344 - setuplight: index error +! = 345 - gsi_fv3ncdf2d_read: problem in reading variable in 2D firstguess of FV3-LAM ! = 899 - foto no longer available ! ! diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 249a653b65..70618120d0 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -177,7 +177,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, l_tuneBE_howv + 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 @@ -505,7 +505,7 @@ module gsimod ! 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 +! (aka howv in GSI code): corp_howv, hwllp_howv ! (in namelist session rapidrefresh_cldsurf) ! !EOP @@ -1568,8 +1568,6 @@ module gsimod ! = 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, & @@ -1591,7 +1589,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, l_tuneBE_howv + 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 cb92bb4009..5c7ff4ed45 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -313,7 +313,6 @@ 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 @@ -370,8 +369,6 @@ 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 ___________________________________________________________________ @@ -419,7 +416,6 @@ 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' @@ -875,7 +871,7 @@ 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,mype) ! static B for howv + 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) @@ -886,18 +882,6 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt 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 @@ -1139,20 +1123,16 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout,mype) 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' !-- first, assign the pre-defined values to corp and hwllp - if ( twodvar_regional ) then + 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 + 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 @@ -1167,141 +1147,18 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout,mype) 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.' + 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 - 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.' + 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 -!++++ -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 736e9af341..0d0dbb2930 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -32,7 +32,7 @@ module rapidrefresh_cldsurf_mod ! (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 +! 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 @@ -450,9 +450,9 @@ subroutine init_rapidrefresh_cldsurf !-- searching for speficif variable in state variable list (reading from anavinfo) do i2=1,ns2d - if ( trim(svars2d(i2))=='howv' ) then + if ( trim(svars2d(i2))=='howv' .or. trim(svars2d(i2))=='HOWV' ) then i_howv_3dda = 1 - if (mype == 0 ) then + 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 diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index ba6ad2ed57..c0347921d5 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -149,8 +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) for 3D analysis -! in prepbufr data file +! 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 @@ -1135,7 +1135,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& 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 reading same dataset twice in 2DRTMA run) +! (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 From 67b02a9cfc70064c0ce3c0fb66d1806f56b7b27a Mon Sep 17 00:00:00 2001 From: Gang Zhao Date: Fri, 15 Sep 2023 04:30:51 +0000 Subject: [PATCH 4/7] Minor modification in src/gsi/gsi_rfv3io_mod.f90 for reading howv in firstguess. --- src/gsi/gsi_rfv3io_mod.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 89388b7bb9..f2c50cde66 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -1100,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) From aa08628661c1e0d4c4bc0a2077b60cc9374b6814 Mon Sep 17 00:00:00 2001 From: Gang Zhao Date: Thu, 21 Sep 2023 14:31:31 +0000 Subject: [PATCH 5/7] minor modifications for analysis of significant wave height based on code review. --- src/gsi/gsi_rfv3io_mod.f90 | 19 +++++++++---------- src/gsi/gsimain.f90 | 1 - src/gsi/m_berror_stats_reg.f90 | 1 - 3 files changed, 9 insertions(+), 12 deletions(-) diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index f2c50cde66..c1b3eebdf0 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -1824,6 +1824,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) !$$$ end documentation block use kinds, only: r_kind,i_kind 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 @@ -1858,7 +1859,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) 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) id_howv integer(i_kind) iret_bcast ! for io_layout > 1 @@ -1910,19 +1911,17 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) ! 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_howv = nf90_inq_varid(gfile_loc,'howv',id_howv) - if ( iret_howv /= nf90_noerr ) then - iret_howv = nf90_inq_varid(gfile_loc,'HOWV',id_howv) ! double check with name in uppercase + 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_howv /= nf90_noerr ) then + if ( iret /= nf90_noerr ) then i_howv_3dda = 0 ! howv does not exist in firstguess, then stop GSI run. - 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,'). Stop running GSI!!!!' - call stop2(345) + call die('gsi_fv3ncdf2d_read','Warning: CANNOT find howv in firstguess, aborting..., iret = ', iret) else i_howv_3dda = 1 ! howv does exist in firstguess, running analysis with howv - 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,').' + 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 diff --git a/src/gsi/gsimain.f90 b/src/gsi/gsimain.f90 index 7be4e3a646..5eaa38286a 100644 --- a/src/gsi/gsimain.f90 +++ b/src/gsi/gsimain.f90 @@ -530,7 +530,6 @@ program gsi ! = 342 - setuplight: failure to allocate obsdiags ! = 343 - setuplight: failure to allocate obsdiags ! = 344 - setuplight: index error -! = 345 - gsi_fv3ncdf2d_read: problem in reading variable in 2D firstguess of FV3-LAM ! = 899 - foto no longer available ! ! diff --git a/src/gsi/m_berror_stats_reg.f90 b/src/gsi/m_berror_stats_reg.f90 index 5c7ff4ed45..601339e1ac 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -313,7 +313,6 @@ 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 hybrid_ensemble_parameters, only: l_hyb_ens implicit none From 78888315afe211b17aada291dc26e79e38caf461 Mon Sep 17 00:00:00 2001 From: Gang Zhao Date: Thu, 21 Sep 2023 15:42:41 +0000 Subject: [PATCH 6/7] trivial modification (fixing a typo) --- src/gsi/rapidrefresh_cldsurf_mod.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/rapidrefresh_cldsurf_mod.f90 b/src/gsi/rapidrefresh_cldsurf_mod.f90 index 0d0dbb2930..954a85d22e 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -187,7 +187,7 @@ module rapidrefresh_cldsurf_mod ! 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_3dda - integer, control the analysis of howv in 3D analysis (either var or hybrid) +! 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; From f6031bdb0d8e77259caa9fbc1e679a825813f94f Mon Sep 17 00:00:00 2001 From: Gang Zhao Date: Tue, 26 Sep 2023 20:49:44 +0000 Subject: [PATCH 7/7] minor modifications based on the feedback from code reviewers on the analysis of siginificant wave height. --- src/gsi/gsi_rfv3io_mod.f90 | 1 - src/gsi/rapidrefresh_cldsurf_mod.f90 | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index c1b3eebdf0..6d16be7c13 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -1919,7 +1919,6 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv) 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 - i_howv_3dda = 1 ! howv does exist in firstguess, running analysis with howv 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 diff --git a/src/gsi/rapidrefresh_cldsurf_mod.f90 b/src/gsi/rapidrefresh_cldsurf_mod.f90 index 954a85d22e..122d2872d0 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -448,7 +448,7 @@ subroutine init_rapidrefresh_cldsurf 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 speficif variable in state variable list (reading from anavinfo) +!-- 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