diff --git a/src/gsi/cplr_read_wrf_mass_guess.f90 b/src/gsi/cplr_read_wrf_mass_guess.f90 index 54db338ac..c10ee0d8b 100644 --- a/src/gsi/cplr_read_wrf_mass_guess.f90 +++ b/src/gsi/cplr_read_wrf_mass_guess.f90 @@ -1394,7 +1394,7 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) aerotot_guess,init_aerotot_guess,wrf_pm2_5,aero_ratios use rapidrefresh_cldsurf_mod, only: l_hydrometeor_bkio,l_gsd_soiltq_nudge use rapidrefresh_cldsurf_mod, only: i_use_2mq4b,i_use_2mt4b - use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda + use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda, i_sfcrough_fgs use wrf_mass_guess_mod, only: soil_temp_cld,isli_cld,ges_xlon,ges_xlat,ges_tten,create_cld_grids use gsi_bundlemod, only: GSI_BundleGetPointer use gsi_metguess_mod, only: gsi_metguess_get,GSI_MetGuess_Bundle @@ -1447,7 +1447,7 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) integer(i_kind) i_qc,i_qi,i_qr,i_qs,i_qg,i_qnr,i_qni,i_qnc,i_w,i_dbz integer(i_kind) kqc,kqi,kqr,kqs,kqg,kqnr,kqni,kqnc,i_xlon,i_xlat,i_tt,ktt integer(i_kind) i_th2,i_q2,i_soilt1,ksmois,ktslb - integer(i_kind) i_howv, i_gust + integer(i_kind) i_howv, i_gust, i_rough integer(i_kind) ier, istatus integer(i_kind) n_actual_clouds integer(i_kind) iv,n_gocart_var @@ -1566,6 +1566,7 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) if(i_use_2mq4b > 0 .and. i_use_2mt4b <=0 ) num_mass_fields=num_mass_fields + 1 if( i_howv_3dda > 0 ) num_mass_fields = num_mass_fields + 1 if( i_gust_3dda > 0 ) num_mass_fields = num_mass_fields + 1 + if( i_sfcrough_fgs > 0 ) num_mass_fields = num_mass_fields + 1 if (laeroana_gocart .and. wrf_pm2_5 ) then if(mype==0) write(6,*)'laeroana_gocart canoot be both true' @@ -1758,6 +1759,12 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) write(identity(i),'("record ",i3,"--gust")')i jsig_skip(i)=0 ; igtype(i)=1 end if + ! for surface roughness (sfcrough is after tsk/q2/soilt1/th2, and before cloud hydrometers) + if ( i_sfcrough_fgs >0 ) then + i=i+1 ; i_rough=i ! roughtness + write(identity(i),'("record ",i3,"--rough")')i + jsig_skip(i)=0 ; igtype(i)=1 + end if ! for cloud array if(l_hydrometeor_bkio .and. n_actual_clouds>0) then i_qc=i+1 @@ -2283,7 +2290,12 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) psfc_this=(psfc_this_dry-pt_ll)*q_integral(j,i)+pt_ll+q_integralc4h(j,i) ges_ps_it(j,i)=one_tenth*psfc_this ! convert from mb to cb sno(j,i,it)=real(all_loc(j,i,i_0+i_sno),r_kind) - sfc_rough(j,i,it)=rough_default + ! surface roughness (ZNT) + if ( i_sfcrough_fgs >0 ) then + sfc_rough(j,i,it) = real(all_loc(j,i,i_0+i_rough),r_kind) + else + sfc_rough(j,i,it) = rough_default + end if if(i_use_2mt4b > 0 ) then ges_t2m_it(j,i)=real(all_loc(j,i,i_0+i_th2),r_kind) ! convert from potential to sensible temperature diff --git a/src/gsi/cplr_regional_io.f90 b/src/gsi/cplr_regional_io.f90 index 62a27b8f9..7b2ffaa4c 100644 --- a/src/gsi/cplr_regional_io.f90 +++ b/src/gsi/cplr_regional_io.f90 @@ -89,7 +89,7 @@ subroutine convert_regional_guess_wrf(this,mype,ctph0,stph0,tlm0) use wrf_params_mod, only: update_pint,cold_start use gsi_io, only: verbose use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda - use rapidrefresh_cldsurf_mod, only: i_howv_mask + use rapidrefresh_cldsurf_mod, only: i_howv_mask, i_sfcrough_fgs implicit none @@ -150,11 +150,12 @@ subroutine convert_regional_guess_wrf(this,mype,ctph0,stph0,tlm0) call mpi_bcast(i_howv_3dda, 1, mpi_itype, 0, mpi_comm_world, ierror) call mpi_bcast(i_gust_3dda, 1, mpi_itype, 0, mpi_comm_world, ierror) call mpi_bcast(i_howv_mask, 1, mpi_itype, 0, mpi_comm_world, ierror) + call mpi_bcast(i_sfcrough_fgs, 1, mpi_itype, 0, mpi_comm_world, ierror) if(print_verbose)write(6,*)' in convert_regional_guess, for wrf arw binary input, byte_swap=',byte_swap if (mype <= 1 .and. print_verbose) then - write(6,'(1x,A,3(2x,I4),2x,A6,I6.6,A2)') & - ' in convert_regional_guess, i_howv_3dda i_gust_3dda i_howv_mask =', & - i_howv_3dda,i_gust_3dda,i_howv_mask,' (pe=',mype,').' + write(6,'(1x,A,4(2x,I4),2x,A6,I6.6,A2)') & + ' in convert_regional_guess, i_howv_3dda i_gust_3dda i_howv_mask i_sfcrough_fgs=', & + i_howv_3dda,i_gust_3dda,i_howv_mask,i_sfcrough_fgs,' (pe=',mype,').' end if elseif (cmaq_regional) then diff --git a/src/gsi/cplr_wrf_netcdf_interface.f90 b/src/gsi/cplr_wrf_netcdf_interface.f90 index 7cce1bd36..a64eb2874 100644 --- a/src/gsi/cplr_wrf_netcdf_interface.f90 +++ b/src/gsi/cplr_wrf_netcdf_interface.f90 @@ -67,7 +67,7 @@ subroutine convert_netcdf_mass_wrf(this) use rapidrefresh_cldsurf_mod, only: l_hydrometeor_bkio,l_gsd_soilTQ_nudge use rapidrefresh_cldsurf_mod, only: i_use_2mt4b,i_use_2mq4b use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda - use rapidrefresh_cldsurf_mod, only: i_howv_mask + use rapidrefresh_cldsurf_mod, only: i_howv_mask, i_sfcrough_fgs use gsi_metguess_mod, only: gsi_metguess_get use chemmod, only: laeroana_gocart, ppmv_conv,wrf_pm2_5 use gsi_chemguess_mod, only: gsi_chemguess_get @@ -1141,6 +1141,35 @@ subroutine convert_netcdf_mass_wrf(this) end if endif ! i_gust_3dda (reading 2D 10-m wind gust from netcdf-format background) +! Reading Surface Roughness (ZNT) in firstguess (it is dumped in sigf file after gust and before cloud variables.) + if ( i_sfcrough_fgs >= 1 ) then + rmse_var='ZNT' + call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, & + start_index,end_index, WrfType, ierr ) + if(print_verbose)then + write(6,*)' rmse_var = ',trim(rmse_var),' ndim1=',ndim1 + write(6,*)' WrfType = ',WrfType,' WRF_REAL=',WRF_REAL,'ierr = ',ierr !DEDE + write(6,*)' ordering = ',trim(ordering),' staggering = ',trim(staggering) + write(6,*)' start_index = ',start_index,' end_index = ',end_index + end if + if(ierr == 0) then + call ext_ncd_read_field(dh1,DateStr1,TRIM(rmse_var), & + field2,WRF_REAL,0,0,0,ordering, & + staggering, dimnames , & + start_index,end_index, & !dom + start_index,end_index, & !mem + start_index,end_index, & !pat + ierr ) + if(print_verbose)write(6,*)'convert_netcdf_mass_wrf:: max,min ZNT(sfc_rough)=',maxval(field2),minval(field2) + write(iunit)field2 !ZNT (Surface Roughness Length) + i_sfcrough_fgs = 1 + else + i_sfcrough_fgs = 0 + write(6,'(1x,A,1x,I4)') & + 'convert_netcdf_mass_wrf:: ZNT (sfc_roughness) is NOT in firstguess. Re-set i_sfcrough_fgs=',i_sfcrough_fgs + end if + end if ! i_sfcrough_fgs (reading surface roughness from netcdf-format background) + if(l_hydrometeor_bkio .and. n_actual_clouds>0) then rmse_var='QCLOUD' call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, & @@ -2556,7 +2585,7 @@ subroutine update_netcdf_mass_wrf(this) use constants, only: h300,tiny_single use rapidrefresh_cldsurf_mod, only: l_hydrometeor_bkio,l_gsd_soilTQ_nudge use rapidrefresh_cldsurf_mod, only: i_gsdcldanal_type - use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda + use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda, i_sfcrough_fgs use gsi_metguess_mod, only: gsi_metguess_get,GSI_MetGuess_Bundle use rapidrefresh_cldsurf_mod, only: i_use_2mt4b,i_use_2mq4b use gsi_bundlemod, only: GSI_BundleGetPointer @@ -3174,8 +3203,14 @@ subroutine update_netcdf_mass_wrf(this) start_index,end_index1, & !mem start_index,end_index1, & !pat ierr ) - endif ! i_howv_3dda (wave height) + endif ! i_gust_3dda (wind gust) +! Reading surface roughtness from binary siganl (only when i_sfcrough_fgs = 1), +! and no need to write it to netcdf analysis file, since it is not updated. + if ( i_sfcrough_fgs == 1 ) then + read(iunit) field2 !ZNT (surface roughness) + if(print_verbose)write(6,*)'update_netcdf_mass_wrf:: max,min ZNT(sfc_rough)=',maxval(field2),minval(field2) + end if ! i_sfcrough_fgs (surface roughness) if (l_hydrometeor_bkio .and. n_actual_clouds>0) then do k=1,nsig_regional diff --git a/src/gsi/cplr_wrwrfmassa.f90 b/src/gsi/cplr_wrwrfmassa.f90 index 0bebf1947..4b7268f6b 100644 --- a/src/gsi/cplr_wrwrfmassa.f90 +++ b/src/gsi/cplr_wrwrfmassa.f90 @@ -1877,7 +1877,7 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) use rapidrefresh_cldsurf_mod, only: l_hydrometeor_bkio,l_gsd_soilTQ_nudge,& i_use_2mq4b,i_use_2mt4b use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda - use rapidrefresh_cldsurf_mod, only: i_howv_mask + use rapidrefresh_cldsurf_mod, only: i_howv_mask, i_sfcrough_fgs use chemmod, only: laeroana_gocart,wrf_pm2_5 use gsi_bundlemod, only: GSI_BundleGetPointer @@ -2971,6 +2971,15 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) end if ! i_gust_3dda > 0 + ! ZNT: surface roughness (no updating to it) + if ( i_sfcrough_fgs > 0 ) then + if(mype == 0) then + ! This is surface roughness + read(lendian_in)temp1 + write(lendian_out)temp1 + end if + end if + ! ! for saving cloud analysis results if(l_hydrometeor_bkio .and. n_actual_clouds>0) then diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index d6d9a9447..ebef93e15 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -184,7 +184,8 @@ module gsimod cld_bld_coverage,cld_clr_coverage,& i_cloud_q_innovation,i_ens_mean,DTsTmax,& i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check, & - corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust, i_howv_mask + corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust, i_howv_mask, & + i_sfcrough_fgs 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 @@ -1619,6 +1620,10 @@ module gsimod ! corp_gust - real, static background error of gust (stddev error) ! hwllp_gust - real, background error de-correlation length scale of gust ! oerr_gust - real, observation error of gust +! i_sfcrough_fgs - integer, option to control the read-in of surface roughness from firstguess +! = 0 : do not read surface roughness from firstguess, +! and use the default value instead (default) +! = 1 : read surface roughness from firstguess and use it in analysis ! namelist/rapidrefresh_cldsurf/dfi_radar_latent_heat_time_period, & metar_impact_radius,metar_impact_radius_lowcloud, & @@ -1640,7 +1645,8 @@ module gsimod cld_bld_coverage,cld_clr_coverage,& i_cloud_q_innovation,i_ens_mean,DTsTmax, & i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check, & - corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust, i_howv_mask + corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust, i_howv_mask, & + i_sfcrough_fgs ! chem(options for gsi chem analysis) : ! berror_chem - .true. when background for chemical species that require diff --git a/src/gsi/rapidrefresh_cldsurf_mod.f90 b/src/gsi/rapidrefresh_cldsurf_mod.f90 index d0c0a0b8e..b9e1f52af 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -221,6 +221,10 @@ module rapidrefresh_cldsurf_mod ! = 0 (gust-off: default) : no analysis of gust in 3D analysis. ! = 1 (gust-on) : if variable name "gust" is found in anavinfo, ! set it to be 1 to turn on analysis of gust; +! i_sfcrough_fgs - integer, namelist option to control the read-in and usage of surface roughness in firstguess +! = 0 : do not read surface roughness from firstguess, +! and use the default value instead (default) +! = 1 : read surface roughness from firstguess and use it in analysis ! ! attributes: ! language: f90 @@ -297,6 +301,7 @@ module rapidrefresh_cldsurf_mod public :: i_howv_mask public :: corp_gust, hwllp_gust, oerr_gust public :: i_gust_3dda + public :: i_sfcrough_fgs logical l_hydrometeor_bkio real(r_kind) dfi_radar_latent_heat_time_period @@ -360,6 +365,7 @@ module rapidrefresh_cldsurf_mod integer(i_kind) :: i_howv_mask real(r_kind) :: corp_gust, hwllp_gust, oerr_gust integer(i_kind) :: i_gust_3dda + integer(i_kind) :: i_sfcrough_fgs contains @@ -495,6 +501,8 @@ subroutine init_rapidrefresh_cldsurf i_gust_3dda = 0 ! no analysis of wind gust (gust) in 3D analysis (default) + i_sfcrough_fgs = 0 ! do not read surface roughness from firstguess (default) + !-- searching for specific variable in state variable list (reading from anavinfo) do i2=1,ns2d if ( trim(svars2d(i2))=='howv' .or. trim(svars2d(i2))=='HOWV' ) then diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index e97c5e516..bbcad853f 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -452,7 +452,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& integer:: icase,klev,ikkk,tkk real:: diffhgt,diffuu,diffvv - integer,dimension(3)::kcount + integer,dimension(5)::kcount real(r_double),dimension(3,1500):: fcstdat logical print_verbose @@ -2066,6 +2066,12 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& call find_wind_height(c_prvstg,c_sprvstg,windsensht,kcount) endif endif + !retrieve wind sensor height for mesonet gustob only when running 3DRTMA + if (l_rtma3d) then + if ( gustob .and. (kx==188.or.kx==195) ) then + call find_wind_height(c_prvstg,c_sprvstg,windsensht,kcount) + endif + endif if(i_gsdqc==2) then ! filter bad 2-m dew point and 0 mesonet wind obs if (kx==288.or.kx==295) then ! for mesonet wind if(abs(obsdat(5,k))<0.01_r_kind .and. abs(obsdat(6,k))<0.01_r_kind) usage=115._r_kind @@ -2572,7 +2578,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if ((kx==280).or.(kx==180)) oelev=r20+selev if ((kx==299).or.(kx==199)) oelev=r20+selev if ((kx==282).or.(kx==182)) oelev=r20+selev - if (((kx==295).or.(kx==288).or.(kx==195).or.(kx==188)).and.twodvar_regional) then + if (((kx==295).or.(kx==288).or.(kx==195).or.(kx==188)).and.(twodvar_regional.or.l_rtma3d)) then !account for mesonet wind sensor height oelev=windsensht+selev end if @@ -3309,7 +3315,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& call destroy_aircraft_rjlists if(i_gsdsfc_uselist==1) call destroy_gsd_sfcuselist if (lhilbert) call destroy_hilbertcurve - if (twodvar_regional) then + if (twodvar_regional .or. l_rtma3d) then call destroy_ndfdgrid call destroy_windht_lists endif @@ -3323,6 +3329,10 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if (twodvar_regional .and. (uvob .or. gustob .or. spdob)) then write(6,*) 'kcount values from find wind height = ',kcount end if + if ( l_rtma3d .and. gustob ) then + write(6,'(1x,A,1x,A,A,5(1x,I9))') 'read_prepbufr:: for obstype=', & + trim(adjustl(obstype)), ': kcount values from find wind height = ',kcount + end if diff --git a/src/gsi/setupgust.f90 b/src/gsi/setupgust.f90 index c6b4aa260..c718d6ce0 100644 --- a/src/gsi/setupgust.f90 +++ b/src/gsi/setupgust.f90 @@ -62,7 +62,7 @@ subroutine setupgust(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag use mpeu_util, only: die,perr use kinds, only: r_kind,r_single,r_double,i_kind - use guess_grids, only: hrdifsig,nfldsig,ges_lnprsl, & + use guess_grids, only: hrdifsig,nfldsig,ges_lnprsl, ges_prsl, & geop_hgtl,sfcmod_gfs,sfcmod_mm5,comp_fact10 use m_obsdiagNode, only: obs_diag use m_obsdiagNode, only: obs_diags @@ -99,6 +99,7 @@ subroutine setupgust(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: l_rtma3d use aux2dvarflds, only: rtma_comp_fact10 implicit none @@ -174,6 +175,18 @@ subroutine setupgust(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag real(r_kind),allocatable,dimension(:,:,:) :: ges_ps real(r_kind),allocatable,dimension(:,:,:) :: ges_z real(r_kind),allocatable,dimension(:,:,:) :: ges_gust + +! arrays used in calcuation of factw in 3DRTMA run + real(r_kind),allocatable,dimension(:,:,:) :: ges_pres1 + real(r_kind),allocatable,dimension(:,:,:) :: ges_pres2 + real(r_kind),allocatable,dimension(:,:,:) :: ges_tv1 + real(r_kind),allocatable,dimension(:,:,:) :: ges_tv2 + real(r_kind),allocatable,dimension(:,:,:) :: ges_q1 + real(r_kind),allocatable,dimension(:,:,:) :: ges_q2 + real(r_kind),allocatable,dimension(:,:,:) :: ges_u1 + real(r_kind),allocatable,dimension(:,:,:) :: ges_v1 +! real(r_kind),allocatable,dimension(:,:,:) :: ges_hgt1 + type(obsLList),pointer,dimension(:):: gusthead gusthead => obsLL(:) @@ -348,6 +361,16 @@ subroutine setupgust(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag dpres=dpres-(dstn+fact*(zsges-dstn)) drpx=0.003*abs(dstn-zsges)*(one-fact) +! Note: caclculating geometric height (zges) at model layer +! mid-points at the observation location. +! But, for 2DRTMA, only use the height at the mid-point of +! bottom layer, and set the height zges(1) = 10 meters. +! Then in 2DRTMA run, there might be two potential issues: +! (1) array zges(1:nsig) is undefined except for zges(1)=10, +! this could cause potential problem in call grdcrd1, +! e.g., line 400 & line 483; +! (2) zges(1)=10 might cause erro of "devided by zero" in line +! 441; if (.not. twodvar_regional) then call tintrp2a1(geop_hgtl,zges,dlat,dlon,dtime,hrdifsig,& nsig,mype,nfldsig) @@ -405,14 +428,18 @@ subroutine setupgust(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag if (zob <= ten) then if(zob < ten)then - if (msonetob .and. twodvar_regional .and. use_similarity_2dvar) then + if (msonetob .and. (twodvar_regional .or. l_rtma3d) .and. use_similarity_2dvar) then if (neutral_stability_windfact_2dvar) then sfcr = data(isfcr,i) factw=log(max(sfcr,zob)/sfcr)/log(ten/sfcr) else sfcr = data(isfcr,i) skint = data(iskint,i) - call rtma_comp_fact10(dlat,dlon,dtime,zob,skint,sfcr,isli,mype,factw) + if(twodvar_regional) then + call rtma_comp_fact10(dlat,dlon,dtime,zob,skint,sfcr,isli,mype,factw) + else if(l_rtma3d) then + call rtma3d_comp_fact10_(dlat,dlon,dtime,zob,skint,sfcr,isli,mype,factw) + end if endif else term = max(zob,zero)/ten @@ -656,8 +683,10 @@ end subroutine check_vars_ subroutine init_vars_ real(r_kind),dimension(:,: ),pointer:: rank2=>NULL() + real(r_kind),dimension(:,:,:),pointer:: rank3=>NULL() character(len=5) :: varname integer(i_kind) ifld, istatus + integer(i_kind) n1, n2, n3 ! If require guess vars available, extract from bundle ... if(size(gsi_metguess_bundle)==nfldsig) then @@ -679,7 +708,7 @@ subroutine init_vars_ write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif -! get ps ... +! get ps ... (unit: cb) varname='ps' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) if (istatus==0) then @@ -715,6 +744,130 @@ subroutine init_vars_ write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif + + if ( l_rtma3d .and. use_similarity_2dvar .and. (.not. neutral_stability_windfact_2dvar) ) then +! get u ... + varname='u' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_u1))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif +! allocate(ges_u(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + n1=size(rank3,1) + n2=size(rank3,2) + n3=size(rank3,3) + allocate(ges_u1(n1,n2,nfldsig)) +! ges_u(:,:,:,1)=rank3 + ges_u1(:,:,1)=rank3(:,:,1) + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) +! ges_u(:,:,:,ifld)=rank3 + ges_u1(:,:,ifld)=rank3(:,:,1) + enddo +! deallocate(ges_u) + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get v ... + varname='v' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_v1))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif +! allocate(ges_v(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + n1=size(rank3,1) + n2=size(rank3,2) + n3=size(rank3,3) + allocate(ges_v1(n1,n2,nfldsig)) +! ges_v(:,:,:,1)=rank3 + ges_v1(:,:,1)=rank3(:,:,1) + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) +! ges_v(:,:,:,ifld)=rank3 + ges_v1(:,:,ifld)=rank3(:,:,1) + enddo +! deallocate(ges_v) + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get tv ... + varname='tv' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_tv1) .or. allocated(ges_tv2))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif +! allocate(ges_tv(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + n1=size(rank3,1) + n2=size(rank3,2) + n3=size(rank3,3) + allocate(ges_tv1(n1,n2,nfldsig)) + allocate(ges_tv2(n1,n2,nfldsig)) +! ges_tv(:,:,:,1)=rank3 + ges_tv1(:,:,1)=rank3(:,:,1) + ges_tv2(:,:,1)=rank3(:,:,2) + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) +! ges_tv(:,:,:,ifld)=rank3 + ges_tv1(:,:,ifld)=rank3(:,:,1) + ges_tv2(:,:,ifld)=rank3(:,:,2) + enddo +! deallocate(ges_tv) + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get q ... + varname='q' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_q1) .or. allocated(ges_q2))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif +! allocate(ges_q(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + n1=size(rank3,1) + n2=size(rank3,2) + n3=size(rank3,3) + allocate(ges_q1(n1,n2,nfldsig)) + allocate(ges_q2(n1,n2,nfldsig)) +! ges_q(:,:,:,1)=rank3 + ges_q1(:,:,1)=rank3(:,:,1) + ges_q2(:,:,1)=rank3(:,:,2) + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) +! ges_q(:,:,:,ifld)=rank3 + ges_q1(:,:,ifld)=rank3(:,:,1) + ges_q2(:,:,ifld)=rank3(:,:,2) + enddo +! deallocate(ges_q) + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get pressure ... ==> ges_prsl from module guess_grids (unit: cb) + varname='pres' + n1=size(ges_prsl,1) + n2=size(ges_prsl,2) + if(allocated(ges_pres1) .or. allocated(ges_pres2))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_pres1(n1,n2,nfldsig)) + allocate(ges_pres2(n1,n2,nfldsig)) + do ifld=1,nfldsig + ges_pres1(:,:,ifld)=ges_prsl(:,:,1,ifld) + ges_pres2(:,:,ifld)=ges_prsl(:,:,2,ifld) + end do +! get height ... ==> geop_hgtl from module guess_grids ==> zges + end if !if (l_rtma3d.and.use_similarity_2dvar.and. (.not.neutral_stability_windfact_2dvar)) else write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',& nfldsig,size(gsi_metguess_bundle) @@ -858,6 +1011,10 @@ subroutine contents_netcdf_diag_(odiag) r_sprvstg = data(isprvd,i) call nc_diag_metadata("Subprovider_Name", c_sprvstg ) + if ( twodvar_regional .or. l_rtma3d ) then + call nc_diag_metadata("Wind_Reduction_Factor_at_10m", factw ) + end if + if (lobsdiagsave) then do jj=1,miter if (odiag%muse(jj)) then @@ -879,7 +1036,126 @@ subroutine final_vars_ if(allocated(ges_z )) deallocate(ges_z ) if(allocated(ges_ps )) deallocate(ges_ps ) if(allocated(ges_gust)) deallocate(ges_gust) + if(allocated(ges_tv1 )) deallocate(ges_pres1 ) + if(allocated(ges_tv2 )) deallocate(ges_pres2 ) + if(allocated(ges_tv1 )) deallocate(ges_tv1 ) + if(allocated(ges_tv2 )) deallocate(ges_tv2 ) + if(allocated(ges_q1 )) deallocate(ges_q1 ) + if(allocated(ges_q2 )) deallocate(ges_q2 ) + if(allocated(ges_u1 )) deallocate(ges_u1 ) + if(allocated(ges_v1 )) deallocate(ges_v1 ) end subroutine final_vars_ +!=============================================================================! + subroutine rtma3d_comp_fact10_(dlat_i,dlon_i,dtime_i,zob_i,skint_i,sfcrough_i, & + islimsk_i,mype_i,factw_o) +!-----------------------------------------------------------------------------! +! !ROUTINE: rtma3d_comp_fact10 --- Compute the factw with similarity theory +! !REVISION HISTORY: +! 2025-04-16 zhao - code initial +! To be in consistency with rtma, this subroutine is mainly based +! on subroutine rtma_comp_fact10 developed by Xiaoyan Zhang & Manuel Pondeca. +!-----------------------------------------------------------------------------! +! +! !USES: + + use guess_grids, only: nfldsig,hrdifsig + use constants, only: zero,one,r100,fv + implicit none + +! Declare passed variables + real(r_kind) , intent(in ) :: dlat_i,dlon_i,dtime_i,skint_i,sfcrough_i,zob_i + integer(i_kind), intent(in ) :: mype_i,islimsk_i + real(r_kind) , intent( out) :: factw_o + + +! Declare local parameters + real(r_kind), parameter :: r0_001 = 0.001_r_kind + character(len=*), parameter :: myname='rtma3d_comput_fact10' + +! Declare external calls for code analysis + external:: intrp2a11,tintrp2a11 + external:: SFC_WTQ_FWD + +! Declare local variables + integer(i_kind) islimsk2 + +! for similarity theory + integer(i_kind) :: regime + logical iqtflg + real(r_kind) :: psfcges,tgges,roges + real(r_kind) :: tv1,tv2 + real(r_kind) :: pres1,pres2,qq1,qq2,uu1,vv1,hgt1 + real(r_kind) :: u10ges,v10ges,t2ges,q2ges + + iqtflg=.true. !.true. means output is virtual temperature + + islimsk2=islimsk_i + if(islimsk2 > 2)islimsk2=islimsk2-3 + + call tintrp2a11(ges_ps, psfcges,dlat_i,dlon_i,dtime_i,hrdifsig,& + mype_i,nfldsig) + call tintrp2a11(ges_pres1,pres1,dlat_i,dlon_i,dtime_i,hrdifsig,& + mype_i,nfldsig) + call tintrp2a11(ges_pres2,pres2,dlat_i,dlon_i,dtime_i,hrdifsig,& + mype_i,nfldsig) + call tintrp2a11(ges_tv1, tv1,dlat_i,dlon_i,dtime_i,hrdifsig,& + mype_i,nfldsig) + call tintrp2a11(ges_tv2, tv2,dlat_i,dlon_i,dtime_i,hrdifsig,& + mype_i,nfldsig) + call tintrp2a11(ges_q1, qq1,dlat_i,dlon_i,dtime_i,hrdifsig,& + mype_i,nfldsig) + call tintrp2a11(ges_q2, qq2,dlat_i,dlon_i,dtime_i,hrdifsig,& + mype_i,nfldsig) + call tintrp2a11(ges_u1, uu1,dlat_i,dlon_i,dtime_i,hrdifsig,& + mype_i,nfldsig) + call tintrp2a11(ges_v1, vv1,dlat_i,dlon_i,dtime_i,hrdifsig,& + mype_i,nfldsig) +! call tintrp2a11(geop_hgtl(:,:,1,:), hgt1,dlat_i,dlon_i,dtime_i,hrdifsig,& +! mype_i,nfldsig) + hgt1 = zges(1) ! <== subroutine setupgust + + !overwrite hgt1. Assume local height of model's first half-sigma level + ! to be zob, that is, the ob height + hgt1=zob_i + + tgges=skint_i + + !convert input pressure variables from Pa to cb (already in unit of cb, no need to convert). +! psfcges = r0_001*psfcges +! pres1 = r0_001*pres1 +! pres2 = r0_001*pres2 + + !convert sensible temperature to virtual temperature (tv1,tv2 are virtual temp already) +! tv1=tmp1*((one+fv*qq1)) +! tv2=tmp2*((one+fv*qq2)) + !!!tgges=tgges*((one+fv*qq1)) !no need for virtual temp conversion + + if (hgt1 <= sfcrough_i) then + factw_o = zero + else + + islimsk2=islimsk_i + if(islimsk2 > 2)islimsk2=islimsk2-3 + + !unit change: originally m --> cm + roges=sfcrough_i*r100 + + call SFC_WTQ_FWD (psfcges, tgges,& + pres1, tv1, qq1, uu1, vv1, & + pres2, tv2, qq2, hgt1, roges, islimsk2, & !why is "comp_fact10" using islimsk instead of islimsk2? + !output variables + factw_o, u10ges, v10ges, t2ges, q2ges, regime, iqtflg) + + if (factw_o > zero) factw_o = one/factw_o !you need this in the 2dvar application + !since the goal is to reduce the + !downscaled background 10-m wind + !to the zob elevation + endif + + return + + end subroutine rtma3d_comp_fact10_ + end subroutine setupgust end module gust_setup diff --git a/src/gsi/windht.f90 b/src/gsi/windht.f90 index bd685155c..7599632b5 100644 --- a/src/gsi/windht.f90 +++ b/src/gsi/windht.f90 @@ -149,7 +149,7 @@ subroutine find_wind_height(cprov,csubprov,finalheight,kcount) character(len=8),intent(in)::cprov,csubprov real(r_kind),intent(out)::finalheight - integer,dimension(3),intent(inout)::kcount + integer,dimension(5),intent(inout)::kcount !local vars integer(i_kind)::i @@ -193,12 +193,14 @@ subroutine find_wind_height(cprov,csubprov,finalheight,kcount) do i=1,nmax if (i>numprovs) then finalheight=r10 + kcount(4) = kcount(4) + 1 return else tmpprov=provlist(i)(1:8) tmpsubprov=provlist(i)(9:16) if (cprov==tmpprov.and.((tmpsubprov==csubprov).or.(tmpsubprov==allprov))) then finalheight=heightlist(i) + kcount(5) = kcount(5) + 1 return endif endif