diff --git a/src/GSD/gsdcloud/cloudCover_Surface.f90 b/src/GSD/gsdcloud/cloudCover_Surface.f90 index 40554e9e4c..f56e43590d 100644 --- a/src/GSD/gsdcloud/cloudCover_Surface.f90 +++ b/src/GSD/gsdcloud/cloudCover_Surface.f90 @@ -1,4 +1,4 @@ -SUBROUTINE cloudCover_Surface(mype,nlat,nlon,nsig,r_radius,thunderRadius,& +SUBROUTINE cloudCover_Surface(mype,nlat,nlon,nsig,thunderRadius,& cld_bld_hgt,t_bk,p_bk,q,h_bk,zh, & mxst_p,NVARCLD_P,numsao,OI,OJ,OCLD,OWX,Oelvtn,Odist,& cld_cover_3d,cld_type_3d,wthr_type,pcp_type_3d, & @@ -23,7 +23,6 @@ SUBROUTINE cloudCover_Surface(mype,nlat,nlon,nsig,r_radius,thunderRadius,& ! nlon - no. of lons on subdomain (buffer points on ends) ! nlat - no. of lats on subdomain (buffer points on ends) ! nsig - no. of levels -! r_radius - influence radius of the cloud observation ! thunderRadius - ! cld_bld_hgt - Height below which cloud building is done ! @@ -71,7 +70,6 @@ SUBROUTINE cloudCover_Surface(mype,nlat,nlon,nsig,r_radius,thunderRadius,& implicit none integer(i_kind),intent(in) :: mype - REAL(r_single), intent(in) :: r_radius integer(i_kind),intent(in) :: nlat,nlon,nsig real(r_single), intent(in) :: thunderRadius real(r_kind), intent(in) :: cld_bld_hgt diff --git a/src/gsi/gsdcloudanalysis.F90 b/src/gsi/gsdcloudanalysis.F90 index 599208ed8c..21fc21b8a2 100644 --- a/src/gsi/gsdcloudanalysis.F90 +++ b/src/gsi/gsdcloudanalysis.F90 @@ -23,7 +23,21 @@ subroutine gsdcloudanalysis(mype) ! 2015-01-13 ladwig - add code for Qni and Qnc (cloud ice and water number concentration) ! 2017-03-23 Hu - add code to use hybrid vertical coodinate in WRF MASS ! core -! +! 2019-10-10 Zhao - add code to check and adjust Qr/qs/qg and Qnr for +! each vertical profile to reduce the background +! reflectivity ghost in final analysis. (for RTMA3D +! only now, option l_precip_vertical_check) +! 2020-04-16 Zhao - modifications to the code which checks and adjusts the vertical +! profile of Qg/Qr/Qs/Qnr retrieved through cloud analysis to +! alleviate the background ghost reflectivity in analysis. +! Modifications includes: +! 1. change option l_precip_vertical_check to i_precip_vertical_check +! 2. i_precip_vertical_check: +! = 0(no adjustment, default) +! = 1(Clean off Qg only, where dbz_obs_max<=35dbz in the profile) +! = 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); ! ! input argument list: ! mype - processor ID that does this IO @@ -66,7 +80,6 @@ subroutine gsdcloudanalysis(mype) use mpimod, only: mpi_comm_world,ierror,mpi_real4 use rapidrefresh_cldsurf_mod, only: dfi_radar_latent_heat_time_period, & metar_impact_radius, & - metar_impact_radius_lowCloud, & l_cleanSnow_WarmTs,l_conserve_thetaV,& r_cleanSnow_WarmTs_threshold, & i_conserve_thetaV_iternum, & @@ -76,16 +89,18 @@ subroutine gsdcloudanalysis(mype) iclean_hydro_withRef, iclean_hydro_withRef_allcol, & l_use_hydroretrieval_all, & i_lightpcp, l_numconc, qv_max_inc,ioption, & - l_precip_clear_only,l_fog_off,cld_bld_coverage,cld_clr_coverage + l_precip_clear_only,l_fog_off,cld_bld_coverage,cld_clr_coverage,& + i_T_Q_adjust,l_saturate_bkCloud,i_precip_vertical_check,l_rtma3d use gsi_metguess_mod, only: GSI_MetGuess_Bundle use gsi_bundlemod, only: gsi_bundlegetpointer + use gsi_io, only: verbose +#ifdef RR_CLOUDANALYSIS implicit none ! Declare passed variables integer(i_kind),intent(in):: mype -#ifdef RR_CLOUDANALYSIS ! ! background ! @@ -182,8 +197,6 @@ subroutine gsdcloudanalysis(mype) real(r_single),allocatable :: vis2qc(:,:) ! fog real(r_kind) :: thunderRadius=2.5_r_kind - real(r_single) :: r_radius ! influence radius of cloud based on METAR obs - real(r_single) :: r_radius_lowCloud ! influence radius of low cloud to cloud top pressure integer(i_kind) :: miss_obs_int real(r_kind) :: miss_obs_real parameter ( miss_obs_int = -99999 ) @@ -235,7 +248,10 @@ subroutine gsdcloudanalysis(mype) logical :: ifindomain integer(i_kind) :: imaxlvl_ref real(r_kind) :: max_retrieved_qrqs,max_bk_qrqs,ratio_hyd_bk2obs + real(r_kind) :: qrqs_retrieved real(r_kind) :: qrlimit,qrlimit_lightpcp + real(r_kind) :: qnr_limit + real(r_kind) :: dbz_clean_graupel character(10) :: obstype integer(i_kind) :: lunin, is, ier, istatus integer(i_kind) :: nreal,nchanl,ilat1s,ilon1s @@ -253,6 +269,9 @@ subroutine gsdcloudanalysis(mype) real(r_kind),parameter :: rho_w = 999.97_r_kind, rho_a = 1.2_r_kind real(r_kind),parameter :: cldDiameter = 10.0E3_r_kind +! local variables used for adjustment of qr/qs for RTMA_3D to alleviate ghost reflectivity + logical :: print_verbose + integer(i_kind) :: k_cap ! highest level when adjument is done (used for adjust qr/qs for RTMA_3D) ! ! @@ -289,8 +308,6 @@ subroutine gsdcloudanalysis(mype) ! ! krad_bot=7.0_r_single - r_radius=metar_impact_radius - r_radius_lowCloud=metar_impact_radius_lowCloud opt_hydrometeor_retri=3 ! 1=Kessler 2=Lin 3=Thompson opt_cloudtemperature=3 ! 3=latent heat, 4,5,6 = adiabat profile @@ -310,6 +327,9 @@ subroutine gsdcloudanalysis(mype) istat_lightning=0 istat_nasalarc=0 + print_verbose=.false. + if (verbose) print_verbose=.true. + call load_gsdpbl_hgt(mype) ! ! check consistency of the options @@ -694,7 +714,7 @@ subroutine gsdcloudanalysis(mype) ! ! if(istat_surface == 1) then - call cloudCover_surface(mype,lat2,lon2,nsig,r_radius,thunderRadius, & + call cloudCover_surface(mype,lat2,lon2,nsig,thunderRadius, & cld_bld_hgt,t_bk,p_bk,q_bk,h_bk,zh, & numsao,nvarcld_p,numsao,oi,oj,ocld,owx,oelvtn,odist, & cld_cover_3d,cld_type_3d,wthr_type_2d,pcp_type_3d, & @@ -869,7 +889,7 @@ subroutine gsdcloudanalysis(mype) ! 2013) ! - if(l_use_hydroretrieval_all) then !RTMA + if(l_use_hydroretrieval_all .or. l_rtma3d) then !RTMA qrlimit=15.0_r_kind*0.001_r_kind do k=1,nsig do j=2,lat2-1 @@ -907,6 +927,103 @@ subroutine gsdcloudanalysis(mype) end do end do end do + +! ---- verical check and adjustment to the analysis of precipitation +! in order to remove/reduce the backround reflectivity "ghost" in +! analysis. +! Note: here rain_3d, snow_3d have been already changed into unit of kg/kg. + if(i_precip_vertical_check > 0) then + + if(print_verbose) then + write(6,'(1x,A,I4.4,A)')"SUB gsdcloudanalysis: precip_vertical_check start... (for pe=",mype,")." + else + if(mype == 0) then + write(6,'(1x,A,I4.4,A)')"SUB gsdcloudanalysis: precip_vertical_check start ... (only print for pe=",mype,")." + end if + end if + + qnr_limit=200000_r_kind + dbz_clean_graupel=35.0_r_kind + + do j=2,lat2-1 + do i=2,lon2-1 + +! 1. search the max reflectivity obs for veritcal profile at each grid +! point (same code used in hydrometeor anlysis for RAP forecast) + refmax=-999.0_r_kind + imaxlvl_ref=0 + do k=1,nsig + if(ref_mos_3d(i,j,k) > refmax) then + imaxlvl_ref=k + refmax=ref_mos_3d(i,j,k) + endif + enddo +! 2. check and adjustment along the profile at each grid point + if( refmax > 0 .and. (imaxlvl_ref > 0 .and. imaxlvl_ref < nsig ) ) then + ! cleaning the Graupel, if refmax <= dbz_clean_graupel (35dbz) + ! because graupel is copied from background, not retrieved in cloud analysis. + ! (as seen above, graupel_3d(i,j,k) = ges_qg(j,i,k) ) + if( refmax <= dbz_clean_graupel ) graupel_3d(i,j,:) = zero + + ! adjusting hydrometeors based on maximum reflectivity level + select case (i_precip_vertical_check) + case(2) ! adjust each level along the profile (1:nsig) + max_retrieved_qrqs=snow_3d(i,j,imaxlvl_ref)+rain_3d(i,j,imaxlvl_ref) + do k=1,nsig + qrqs_retrieved=snow_3d(i,j,k)+rain_3d(i,j,k) + if(qrqs_retrieved > max_retrieved_qrqs .and. qrqs_retrieved > 0.0001_r_kind) then + ratio_hyd_bk2obs=max(min(max_retrieved_qrqs/qrqs_retrieved,1.0_r_kind),0.0_r_kind) + if(rain_3d(i,j,k) > zero) then + rain_3d(i,j,k) = rain_3d(i,j,k)*ratio_hyd_bk2obs + nrain_3d(i,j,k)= min(nrain_3d(i,j,k)/ratio_hyd_bk2obs*2.5_r_kind,qnr_limit) + endif + if(snow_3d(i,j,k) > zero) then + snow_3d(i,j,k) = snow_3d(i,j,k)*ratio_hyd_bk2obs + end if + end if + end do + case(3) ! adjust the dbz-obs-missed levels below max-dbz layer (1:kcap) + ! based on the qr+qs on max-refl level + ! keep the retrieved cloud analysis as much as possible + max_retrieved_qrqs=snow_3d(i,j,imaxlvl_ref)+rain_3d(i,j,imaxlvl_ref) + k_cap=min(imaxlvl_ref,nsig) + do k=k_cap,1,-1 + if( ref_mos_3d(i,j,k) <= -100.0_r_kind ) then ! dbz-obs-missing level + qrqs_retrieved=snow_3d(i,j,k)+rain_3d(i,j,k) + if(qrqs_retrieved > max_retrieved_qrqs .and. qrqs_retrieved > 0.0001_r_kind) then + ratio_hyd_bk2obs=max(min(max_retrieved_qrqs/qrqs_retrieved,1.0_r_kind),0.0_r_kind) + if(rain_3d(i,j,k) > zero) then + rain_3d(i,j,k) = rain_3d(i,j,k)*ratio_hyd_bk2obs + ! for nrain_3d: 2.5(old) or 1.0(new4) or 1.5(new5/6) 2.5(old, new7) 2.0(new8) + nrain_3d(i,j,k)= min(nrain_3d(i,j,k)/ratio_hyd_bk2obs*2.5_r_kind,qnr_limit) + endif + if(snow_3d(i,j,k) > zero) then + snow_3d(i,j,k) = snow_3d(i,j,k)*ratio_hyd_bk2obs + end if + end if + end if + end do + case default + rain_3d(i,j,k) = rain_3d(i,j,k) + nrain_3d(i,j,k)= nrain_3d(i,j,k) + snow_3d(i,j,k) = snow_3d(i,j,k) + end select + + end if + + end do + end do + + if(print_verbose) then + write(6,'(1x,A,I4.4,A)')"SUB gsdcloudanalysis: precip_vertical_check is done ... (for pe=",mype,")." + else + if(mype == 0) then + write(6,'(1x,A,I4.4,A)')"SUB gsdcloudanalysis: precip_vertical_check is done ... (only print for pe=",mype,")." + end if + end if + + end if + elseif(l_precip_clear_only) then !only clear for HRRRE do k=1,nsig do j=2,lat2-1 @@ -1099,7 +1216,7 @@ subroutine gsdcloudanalysis(mype) ! call cloud_saturation(mype,l_conserve_thetaV,i_conserve_thetaV_iternum, & lat2,lon2,nsig,q_bk,t_bk,p_bk, & - cld_cover_3d,wthr_type_2d,cldwater_3d,cldice_3d,sumqci,qv_max_inc) + cld_cover_3d,wthr_type_2d,cldwater_3d,cldice_3d,sumqci,qv_max_inc, l_saturate_bkCloud) ! @@ -1136,24 +1253,50 @@ subroutine gsdcloudanalysis(mype) do k=1,nsig do j=1,lat2 do i=1,lon2 - if(l_conserve_thetaV) then + ! T/Q update + ! =0 no T/Q adjustment + if(i_T_Q_adjust==0) then + if(mype==0) then + write(6,*) 'gsdcloudanalysis: no T/Q adjustment',mype + endif + ! =1 default T/Q adjustment + elseif(i_T_Q_adjust==1) then if(.not.twodvar_regional .or. .not.tsensible) then - ges_tv(j,i,k)=t_bk(i,j,k)*(p_bk(i,j,k)/h1000)**rd_over_cp * & ! t_bk is potential T - (one+fv*q_bk(i,j,k)) ! convert T to virtual T + ! t_bk is potential T, convert to virtual T + ges_tv(j,i,k)=t_bk(i,j,k)*(p_bk(i,j,k)/h1000)**rd_over_cp * (one+fv*q_bk(i,j,k)) ges_tsen(j,i,k,itsig) = ges_tv(j,i,k)/(one+fv*q_bk(i,j,k)) else - ges_tsen(j,i,k,itsig)=t_bk(i,j,k)*(p_bk(i,j,k)/h1000)**rd_over_cp ! t_bk is potential T - ges_tv(j,i,k) = ges_tsen(j,i,k,itsig)*(one+fv*q_bk(i,j,k)) ! convert virtual T to T + ! t_bk is potential T, convert virtual T to T + ges_tsen(j,i,k,itsig)=t_bk(i,j,k)*(p_bk(i,j,k)/h1000)**rd_over_cp + ges_tv(j,i,k) = ges_tsen(j,i,k,itsig)*(one+fv*q_bk(i,j,k)) endif - endif ! l_conserve_thetaV - ges_q(j,i,k)=q_bk(i,j,k)/(1+q_bk(i,j,k)) ! Here q is mixing ratio kg/kg, - ! need to convert to specific humidity + ! Here q is mixing ratio kg/kg, need to convert to specific humidity + ges_q(j,i,k)=q_bk(i,j,k)/(1+q_bk(i,j,k)) + ! =2 T/Q adjustment only for case of clearing + elseif(i_T_Q_adjust==2) then + if(.not.twodvar_regional .or. .not.tsensible) then + ! t_bk is potential T, convert to virtual T + ges_tv(j,i,k)=max(ges_tv(j,i,k),t_bk(i,j,k)*(p_bk(i,j,k)/h1000)**rd_over_cp * (one+fv*q_bk(i,j,k))) + ges_tsen(j,i,k,itsig) = max(ges_tsen(j,i,k,itsig),ges_tv(j,i,k)/(one+fv*q_bk(i,j,k))) + else + ! t_bk is potential T, convert virtual T to T + ges_tsen(j,i,k,itsig)= max(ges_tsen(j,i,k,itsig),t_bk(i,j,k)*(p_bk(i,j,k)/h1000)**rd_over_cp) + ges_tv(j,i,k) = max(ges_tv(j,i,k),ges_tsen(j,i,k,itsig)*(one+fv*q_bk(i,j,k))) + endif + ! Here q is mixing ratio kg/kg, need to convert to specific humidity + ges_q(j,i,k)=min(ges_q(j,i,k),q_bk(i,j,k)/(1+q_bk(i,j,k))) + else + write(6,*) 'gsdcloudanalysis: WARNING no T/Q adjustment, check i_T_Q_adjust value',mype + endif + + ! hydrometeor update ges_qr(j,i,k)=rain_3d(i,j,k) ges_qs(j,i,k)=snow_3d(i,j,k) ges_qg(j,i,k)=graupel_3d(i,j,k) ges_ql(j,i,k)=cldwater_3d(i,j,k) ges_qi(j,i,k)=cldice_3d(i,j,k) ges_qnr(j,i,k)=nrain_3d(i,j,k) + ! cloud number concentration update if( l_numconc ) then ges_qni(j,i,k)=nice_3d(i,j,k) ges_qnc(j,i,k)=nwater_3d(i,j,k) @@ -1193,9 +1336,13 @@ subroutine gsdcloudanalysis(mype) endif #else /* Start no RR cloud analysis library block */ + implicit none + +! Declare passed variables + integer(i_kind),intent(in):: mype +! if( mype == 0) write(6,*)'gsdcloudanalysis: dummy routine, does nothing!' #endif /* End no RR cloud analysis library block */ - end subroutine gsdcloudanalysis diff --git a/src/gsi/gsdcloudanalysis4gfs.F90 b/src/gsi/gsdcloudanalysis4gfs.F90 index c778471120..294ae32cf4 100644 --- a/src/gsi/gsdcloudanalysis4gfs.F90 +++ b/src/gsi/gsdcloudanalysis4gfs.F90 @@ -78,11 +78,11 @@ subroutine gsdcloudanalysis4gfs(mype) use gsi_metguess_mod, only: GSI_MetGuess_Bundle use gsi_bundlemod, only: gsi_bundlegetpointer +#ifdef RR_CLOUDANALYSIS implicit none ! Declare passed variables integer(i_kind),intent(in):: mype -#ifdef RR_CLOUDANALYSIS ! ! background ! @@ -931,10 +931,15 @@ subroutine gsdcloudanalysis4gfs(mype) write(6,*) 'gsdcloudanalysis: generalized cloud analysis finished:',mype write(6,*) '========================================' endif + #else /* Start no RR cloud analysis library block */ + implicit none + +! Declare passed variables + integer(i_kind),intent(in):: mype +! if( mype == 0) write(6,*)'gsdcloudanalysis: dummy routine, does nothing!' #endif /* End no RR cloud analysis library block */ - end subroutine gsdcloudanalysis4gfs diff --git a/src/gsi/gsdcloudlib_pseudoq_mod.f90 b/src/gsi/gsdcloudlib_pseudoq_mod.f90 index 5afba70452..b7544a860c 100644 --- a/src/gsi/gsdcloudlib_pseudoq_mod.f90 +++ b/src/gsi/gsdcloudlib_pseudoq_mod.f90 @@ -32,7 +32,8 @@ module gsdcloudlib_pseudoq_mod contains -SUBROUTINE cloudCover_Surface_col(mype,nsig,& +SUBROUTINE cloudCover_Surface_col(mype,nsig, & + i_cloud_q_innovation,& cld_bld_hgt,h_bk,zh, & NVARCLD_P,ocld,Oelvtn,& wthr_type,pcp_type_obs, & @@ -57,6 +58,7 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& ! input argument list: ! mype - processor ID ! nsig - no. of levels +! i_cloud_q_innovation - flag to control building/clearing/both ! cld_bld_hgt - Height below which cloud building is done ! ! h_bk - 3D background height (m) @@ -96,6 +98,7 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& integer(i_kind),intent(in) :: mype integer(i_kind),intent(in) :: nsig + integer(i_kind),intent(in) :: i_cloud_q_innovation real(r_kind), intent(in) :: cld_bld_hgt ! ! surface observation @@ -124,7 +127,7 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& INTEGER(i_kind) :: k INTEGER(i_kind) :: ic integer(i_kind) :: firstcloud,cl_base_broken_k,obused - integer(i_kind) :: kcld + integer(i_kind) :: kcld,kclr real(r_single) :: underlim REAL(r_kind) :: zdiff REAL(r_kind) :: zlev_clr,cloud_dz,cl_base_ista,betav @@ -140,6 +143,7 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& firstcloud = 0 obused =0 kcld=-9 + kclr=99 ! !***************************************************************** ! analysis of surface/METAR cloud observations @@ -161,6 +165,15 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& endif enddo + ! cloud clearing obs + if(i_cloud_q_innovation==20 .or. i_cloud_q_innovation==22) then + do k=3,nsig,5 + if (h_bk(k) < zlev_clr) then + cld_cover_obs(k)=0.0_r_single + endif + enddo + endif + ! -- Now consider non-clear obs ! -------------------------- else @@ -203,10 +216,12 @@ SUBROUTINE cloudCover_Surface_col(mype,nsig,& if(k==8) underlim=95.0_r_kind ! 3000 feet if(k>=9 .and. k= 1.0_r_kind .and. (firstcloud==0 .or. abs(zdiff)= 3) cld_cover_obs(kclr)=0.0_r_single + endif - endif ! end if ocld valid endif ! obused enddo ! end IC loop endif ! end if cloudy ob diff --git a/src/gsi/gsi_cldtotOper.F90 b/src/gsi/gsi_cldtotOper.F90 index 5a63e24765..11516daac3 100644 --- a/src/gsi/gsi_cldtotOper.F90 +++ b/src/gsi/gsi_cldtotOper.F90 @@ -107,7 +107,7 @@ subroutine setup_(self, lunin, mype, is, nobs, init_pass,last_pass) diagsave = write_diag(jiter) .and. diag_conv select case(i_cloud_q_innovation) - case(2) + case(20, 21, 22) call setup(self%obsLL(:), self%odiagLL(:), & lunin,mype,bwork,awork(:,iwork),nele,nobs,is,diagsave) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 34fbd38ea0..95ff7f79b4 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -152,6 +152,9 @@ module gsimod use rapidrefresh_cldsurf_mod, only: init_rapidrefresh_cldsurf, & dfi_radar_latent_heat_time_period,metar_impact_radius,& metar_impact_radius_lowcloud,l_gsd_terrain_match_surftobs, & + l_metar_impact_radius_change, & + metar_impact_radius_max,metar_impact_radius_min,& + metar_impact_radius_max_height,metar_impact_radius_min_height,& l_sfcobserror_ramp_t, l_sfcobserror_ramp_q, & l_pbl_pseudo_surfobst,l_pbl_pseudo_surfobsq,l_pbl_pseudo_surfobsuv, & pblh_ration,pps_press_incr,l_gsd_limit_ocean_q, & @@ -165,7 +168,8 @@ module gsimod i_lightpcp,i_sfct_gross,l_use_hydroretrieval_all,l_numconc,l_closeobs,& 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_cloud_q_innovation,i_ens_mean,DTsTmax,& + i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check 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 @@ -468,10 +472,18 @@ module gsimod ! 2021-01-05 x.zhang/lei - add code for updating delz analysis in regional da ! 09-07-2020 CAPS Add options for directDA_radaruse_mod to use direct radar DA capabilities ! 02-09-2021 CAPS(J. Park) Add vad_near_analtime flag (obsqc) to assimilate newvad obs around analysis time only +! 10-10-2019 Zhao added options l_rtma3d and l_precip_vertical_check +! (adjustment to the cloud-analysis retrieved profile of +! Qg/Qs/Qr/QnrQto to alleviate the reflectivity ghost in +! RTMA3D.) +! 04-16-2020 Zhao change option l_precip_vertical_check to i_precip_vertical_check +! 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) ! 01-07-2022 Hu Add fv3_io_layout_y to let fv3lam interface read/write subdomain restart ! files. The fv3_io_layout_y needs to match fv3lam model ! option io_layout(2). -! +! !EOP !------------------------------------------------------------------------- @@ -1328,6 +1340,17 @@ module gsimod ! enhancement for RR appilcation ): ! dfi_radar_latent_heat_time_period - DFI forward integration window in minutes ! metar_impact_radius - metar low cloud observation impact radius in grid number +! l_metar_impact_radius_change - if .true. the impact radius will change +! with height that set up with the metar_impact_radius_max, min, +! max_height, min_height, (default:false) +! metar_impact_radius_max - The max impact radius of metar cloud observation +! in meter (default: 100000 m). +! metar_impact_radius_min - The min impact radius of metar cloud observation +! in meter (default: 10000 m). +! metar_impact_radius_max_height - The hight above which metar_impact_radius_max apply +! in meter (default: 1200m). +! metar_impact_radius_min_height - The hight below which metar_impact_radius_min apply +! in meter (default: 200m). ! l_gsd_terrain_match_surftobs - if .true., GSD terrain match for surface temperature observation ! l_sfcobserror_ramp_t - namelist logical for adjusting surface temperature observation error ! l_sfcobserror_ramp_q - namelist logical for adjusting surface moisture observation error @@ -1421,7 +1444,9 @@ module gsimod ! i_cloud_q_innovation - integer to choose if and how cloud obs are used ! 0= no innovations ! 1= cloud total innovations -! 2= water vapor innovations +! 20= cloud build/clear derived water vapor innovations +! 21= cloud build derived water vapor innovations +! 22= cloud clear derived water vapor innovations ! 3= cloud total & water vapor innovations ! i_ens_mean - integer for setupcldtot behavior ! 0=single model run @@ -1429,10 +1454,28 @@ module gsimod ! 2=ensemble members ! DTsTmax - maximum allowed difference between Tskin and the first ! level T. This is to safety guard soil T adjustment. +! i_T_Q_adjust - =0 no temperature and moisture adjustment in hydrometeor analyis +! =1 (default) temperature and moisture are adjusted in hydrometeor analyis +! =2 temperature and moisture only adjusted for clearing (warmer, drier) +! l_saturate_bkCloud - if .true. ensure saturation for all cloud 3-d points in background +! where observed cloud cover is missing (default:true). +! l_rtma3d - logical option for turning on configuration for RTMA3D +! (default is .FALSE.) +! i_precip_vertical_check - integer option for checking and adjusting +! Qr/Qs/Qg and Qnr after cloud analysis +! to reduce the background reflectivity ghost in +! analysis. (default is 0) +! = 0(no adjustment) +! = 1(Clean off Qg only, where dbz_obs_max<=35dbz in the profile) +! = 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); ! namelist/rapidrefresh_cldsurf/dfi_radar_latent_heat_time_period, & metar_impact_radius,metar_impact_radius_lowcloud, & - l_gsd_terrain_match_surftobs, & + l_metar_impact_radius_change,metar_impact_radius_max,& + metar_impact_radius_min,metar_impact_radius_max_height,& + metar_impact_radius_min_height,l_gsd_terrain_match_surftobs, & l_sfcobserror_ramp_t,l_sfcobserror_ramp_q, & l_pbl_pseudo_surfobst,l_pbl_pseudo_surfobsq,l_pbl_pseudo_surfobsuv, & pblh_ration,pps_press_incr,l_gsd_limit_ocean_q, & @@ -1446,7 +1489,8 @@ module gsimod i_lightpcp,i_sfct_gross,l_use_hydroretrieval_all,l_numconc,l_closeobs,& 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_cloud_q_innovation,i_ens_mean,DTsTmax, & + i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check ! 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 60c2bd84e6..1ee35fffba 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -22,6 +22,13 @@ module rapidrefresh_cldsurf_mod ! 04-01-2017 Hu added option i_gsdqc to turn on special observation qc ! from GSD (for RAP/HRRR application) ! 2018-09-12 Ladwig added options l_precip_clear_only +! 2019-10-10 Zhao added options l_rtma3d and l_precip_vertical_check (for +! RTMA3D only now) +! 2020-04-16 Zhao change option l_precip_vertical_check to i_precip_vertical_check +! 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) +! ! ! Subroutines Included: ! sub init_rapidrefresh_cldsurf - initialize RR related variables to default values @@ -33,6 +40,20 @@ module rapidrefresh_cldsurf_mod ! def metar_impact_radius - impact radius for METAR cloud observation ! def metar_impact_radius_lowCloud - impact radius for METAR cloud observation ! that indicate low cloud base +! def l_metar_impact_radius_change - if .true. the impact radius will change +! with height that set up with the +! metar_impact_radius_max, min, max_height, +! min_height, (default:false) +! def metar_impact_radius_max - The max impact radius of metar cloud +! observation in meter (default: 100000 m). +! def metar_impact_radius_min - The min impact radius of metar cloud +! observation in meter (default: 10000 m). +! def metar_impact_radius_max_height - The hight above which +! metar_impact_radius_max apply +! in meter (default: 1200m). +! def metar_impact_radius_min_height - The hight below which +! metar_impact_radius_min apply +! in meter (default: 200m). ! def l_gsd_terrain_match_surfTobs - namelist logical for GSD terrain ! match for surface temperature observation ! def l_sfcobserror_ramp_t - namelist logical for adjusting surface temperature observation error @@ -135,12 +156,31 @@ module rapidrefresh_cldsurf_mod ! i_cloud_q_innovation - integer to choose if and how cloud obs are used ! 0= no innovations ! 1= cloud total innovations -! 2= water vapor innovations +! 20= cloud build/clear derived water vapor innovations +! 21= cloud build derived water vapor innovations +! 22= cloud clear derived water vapor innovations ! 3= cloud total & water vapor innovations ! i_ens_mean - integer for setupcldtot behavior ! 0=single model run ! 1=ensemble mean ! 2=ensemble members +! i_T_Q_adjust - =0 no temperature and moisture adjustment in hydrometeor analyis +! =1 (default) temperature and moisture are adjusted in hydrometeor analyis +! =2 temperature and moisture only adjusted for clearing (warmer, drier) +! l_saturate_bkCloud - if .true. ensure saturation for all cloud 3-d points +! in background where observed cloud cover is missing +! (default:true). +! l_rtma3d - logical option for turning on configuration for RTMA3D +! (default is .FALSE.) +! i_precip_vertical_check - integer option for checking and adjusting +! Qr/Qs/Qg and Qnr after precipitation analysis +! to reduce the background reflectivity ghost in +! analysis. (default is 0) +! = 0(no adjustment) +! = 1(Clean off Qg only, where dbz_obs_max<=35dbz in the profile) +! = 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); ! ! attributes: ! language: f90 @@ -159,6 +199,11 @@ module rapidrefresh_cldsurf_mod public :: dfi_radar_latent_heat_time_period public :: metar_impact_radius public :: metar_impact_radius_lowCloud + public :: l_metar_impact_radius_change + public :: metar_impact_radius_max + public :: metar_impact_radius_min + public :: metar_impact_radius_max_height + public :: metar_impact_radius_min_height public :: l_gsd_terrain_match_surfTobs public :: l_sfcobserror_ramp_t public :: l_sfcobserror_ramp_q @@ -203,11 +248,20 @@ module rapidrefresh_cldsurf_mod public :: i_cloud_q_innovation public :: i_ens_mean public :: DTsTmax + public :: i_T_Q_adjust + public :: l_saturate_bkCloud + public :: l_rtma3d + public :: i_precip_vertical_check logical l_hydrometeor_bkio real(r_kind) dfi_radar_latent_heat_time_period real(r_kind) metar_impact_radius real(r_kind) metar_impact_radius_lowCloud + logical l_metar_impact_radius_change + real(r_kind) metar_impact_radius_max + real(r_kind) metar_impact_radius_min + real(r_kind) metar_impact_radius_max_height + real(r_kind) metar_impact_radius_min_height logical l_gsd_terrain_match_surfTobs logical l_sfcobserror_ramp_t logical l_sfcobserror_ramp_q @@ -252,6 +306,10 @@ module rapidrefresh_cldsurf_mod integer(i_kind) i_cloud_q_innovation integer(i_kind) i_ens_mean real(r_kind) DTsTmax + integer(i_kind) i_T_Q_adjust + logical l_saturate_bkCloud + logical l_rtma3d + integer(i_kind) i_precip_vertical_check contains @@ -293,6 +351,11 @@ subroutine init_rapidrefresh_cldsurf dfi_radar_latent_heat_time_period = 30.0_r_kind ! in minutes metar_impact_radius = 10.0_r_kind ! in grid metar_impact_radius_lowCloud = 4.0_r_kind ! in grid + l_metar_impact_radius_change = .false. ! .true. =radius change vertically + metar_impact_radius_max = 50000.0_r_kind ! in meter + metar_impact_radius_min = 20000.0_r_kind ! in meter + metar_impact_radius_max_height = 3000.0_r_kind ! in meter + metar_impact_radius_min_height = 200.0_r_kind ! in meter l_gsd_terrain_match_surfTobs = .false. ! .true. = turn on GSD terrain ! match for surface ! temperature observation @@ -351,6 +414,11 @@ subroutine init_rapidrefresh_cldsurf i_cloud_q_innovation = 0 ! 0 = no increments from cloud obs i_ens_mean = 0 ! typical ob behavior DTsTmax = 20.0_r_kind ! maximum allowed difference between Ts and T 1st level + i_T_Q_adjust= 1 ! temperature and moisture are adjusted + 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) + return end subroutine init_rapidrefresh_cldsurf diff --git a/src/gsi/read_Lightning.f90 b/src/gsi/read_Lightning.f90 index 28f2814c9b..00dbcc3f69 100644 --- a/src/gsi/read_Lightning.f90 +++ b/src/gsi/read_Lightning.f90 @@ -402,11 +402,6 @@ subroutine read_lightning_grid(nread,ndata,infile,obstype,lunout,twind,sis,nobs) call stop2(50) endif loop_report: do while (ireadsb(lunin) == 0) - ntb = ntb+1 - if (ntb>maxobs) then - write(6,*)'read_lightning: reports exceed maximum ',maxobs - call stop2(50) - endif ! Extract type, date, and location information call ufbint(lunin,hdr,5,1,iret,hdrstr) @@ -426,6 +421,12 @@ subroutine read_lightning_grid(nread,ndata,infile,obstype,lunout,twind,sis,nobs) cycle loop_report endif endif +! + ntb = ntb+1 + if (ntb>maxobs) then + write(6,*)'read_lightning: reports exceed maximum ',maxobs + call stop2(50) + endif ! read in observations call ufbint(lunin,obs,1,35,iret,obsstr) diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index b08ea7f094..626952bf9d 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -1562,7 +1562,7 @@ subroutine read_obs(ndata,mype) if(i_gsdcldanal_type==2) then call read_NASA_LaRC_cloud(nread,npuse,nouse,infile,obstype,lunout,sis,nobs_sub1(1,i)) else if(i_gsdcldanal_type==1 .or. i_gsdcldanal_type==6 & - .or. i_gsdcldanal_type==3 .or. i_gsdcldanal_type==7) then + .or. i_gsdcldanal_type==3 .or. i_gsdcldanal_type==7 .or. i_gsdcldanal_type==99) then call read_nasa_larc(nread,npuse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i)) end if string='READ_NASA_LaRC' diff --git a/src/gsi/reorg_metar_cloud.f90 b/src/gsi/reorg_metar_cloud.f90 index f21400910f..2d947f2ef3 100644 --- a/src/gsi/reorg_metar_cloud.f90 +++ b/src/gsi/reorg_metar_cloud.f90 @@ -42,8 +42,12 @@ subroutine reorg_metar_cloud(cdata,nreal,ndata,cdata_all,maxobs,ngrid) use kinds, only: r_kind,i_kind,r_double use gridmod, only: nlon,nlat - use constants, only: one,half + use gridmod, only: region_dx + use constants, only: one,half,zero use rapidrefresh_cldsurf_mod, only: metar_impact_radius + use rapidrefresh_cldsurf_mod, only: l_metar_impact_radius_change, & + metar_impact_radius_max,metar_impact_radius_min,& + metar_impact_radius_max_height,metar_impact_radius_min_height implicit none @@ -59,6 +63,8 @@ subroutine reorg_metar_cloud(cdata,nreal,ndata,cdata_all,maxobs,ngrid) ! integer(i_kind) :: ista_prev,ista_prev2,ista_save + real(r_kind),dimension(nreal) :: cdata_temp + real(r_kind),dimension(12) :: cloudlevel_temp INTEGER(i_kind),allocatable :: first_sta(:,:) INTEGER(i_kind),allocatable :: next_sta(:) INTEGER(i_kind) :: null_p @@ -84,10 +90,18 @@ subroutine reorg_metar_cloud(cdata,nreal,ndata,cdata_all,maxobs,ngrid) real(r_kind) :: spval_p parameter (spval_p = 99999._r_kind) + real(r_kind) :: mindx,dxij,delat_radius,delta_height,radiusij + integer(i_kind) :: isprdij ! ! isprd=int(metar_impact_radius + half) + if(l_metar_impact_radius_change) then + mindx=minval(region_dx) + isprd=int(metar_impact_radius_max/mindx + half) + delat_radius=metar_impact_radius_max-metar_impact_radius_min + delta_height=metar_impact_radius_max_height-metar_impact_radius_min_height + endif if(isprd <= 0) return if(ndata <= 0) return ngrid = 0 @@ -260,19 +274,54 @@ subroutine reorg_metar_cloud(cdata,nreal,ndata,cdata_all,maxobs,ngrid) if (min_dist < 1.e9_r_kind) then if (i1 > 1 .and. i1 < nlon .and. j1 > 1 .and. j1 < nlat) then - iout = iout + 1 - if(iout > maxobs) then - write(6,*)'reorg_metar_cloud: ***Error*** ndata > maxobs ' - call stop2(50) - end if + min_dist=sqrt(min_dist) do k=1,nreal - cdata_all(k,iout) = cdata(k,ista_min) + cdata_temp(k)=cdata(k,ista_min) enddo - cdata_all(24,iout) = cdata_all(2,iout) ! save observaion station i - cdata_all(25,iout) = cdata_all(3,iout) ! save observaion station j - cdata_all(2,iout) = float(i1) ! grid index i - cdata_all(3,iout) = float(j1) ! grid index j - cdata_all(23,iout)=sqrt(min_dist) ! distance from station + if(l_metar_impact_radius_change) then + dxij=region_dx(j1,i1) + + cloudlevel_temp=-99999.0_r_kind + ic=0 + do k=1,6 + if(cdata_temp(5+k) > zero .and. cdata_temp(11+k) > zero ) then + radiusij=metar_impact_radius_min + delat_radius* & + max(min((cdata_temp(11+k)-metar_impact_radius_min_height)/delta_height,one),zero) + isprdij=int(radiusij/dxij+half) + if(min_dist <= isprdij) then + ic=ic+1 + cloudlevel_temp(0+ic)=cdata_temp(5+k) + cloudlevel_temp(6+ic)=cdata_temp(11+k) + endif + endif + enddo + if(ic==0) then + if(abs(cdata_temp(6)) < 0.0001_r_kind) ic=1 ! clear + else + do k=1,6 + cdata_temp(5+k)=cloudlevel_temp(0+k) + cdata_temp(11+k)=cloudlevel_temp(6+k) + enddo + endif + else + ic=1 + endif + + if(ic > 0) then + iout = iout + 1 + if(iout > maxobs) then + write(6,*)'reorg_metar_cloud: ***Error*** ndata > maxobs ' + call stop2(50) + end if + do k=1,nreal + cdata_all(k,iout) = cdata_temp(k) + enddo + cdata_all(24,iout) = cdata_all(2,iout) ! save observaion station i + cdata_all(25,iout) = cdata_all(3,iout) ! save observaion station j + cdata_all(2,iout) = float(i1) ! grid index i + cdata_all(3,iout) = float(j1) ! grid index j + cdata_all(23,iout)= min_dist ! distance from station + endif endif endif enddo ! j1 diff --git a/src/gsi/setupcldtot.F90 b/src/gsi/setupcldtot.F90 index 565075bec3..3d899d1a82 100755 --- a/src/gsi/setupcldtot.F90 +++ b/src/gsi/setupcldtot.F90 @@ -12,7 +12,7 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di !! . . . . ! subprogram: setupcldtot compute rhs of oi for pseudo moisture observations from ! METAR and Satellite cloud observations -! prgmmr: Ladwag org: GSD date: 2019-06-01 +! prgmmr: Ladwig org: GSD date: 2019-06-01 ! ! abstract: For moisture observations, this routine ! a) reads obs assigned to given mpi task (geographic region), @@ -272,7 +272,7 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di allocate(cdiagbuf(nobs*nsig),rdiagbuf(nreal,nobs*nsig)) rdiagbuf=zero endif - if (i_cloud_q_innovation == 2 .or. i_cloud_q_innovation == 3) then + if (i_cloud_q_innovation >= 20 .or. i_cloud_q_innovation == 3) then iip=0 allocate(cdiagbufp(nobs*nsig),rdiagbufp(nreal,nobs*nsig)) cdiagbufp="EMPTY" @@ -460,7 +460,7 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di cycle endif - call cloudCover_surface_col(mype,nsig,cld_bld_hgt,h_bk,z_bk, & + call cloudCover_surface_col(mype,nsig,i_cloud_q_innovation,cld_bld_hgt,h_bk,z_bk, & nvarcld_p,ocld,oelvtn,wthr_type,pcp_type_obs,vis2qc,cld_cover_obs) @@ -515,8 +515,8 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di muse(i)=.true. !******************************************************************************* - if (i_cloud_q_innovation /= 2) then - write(*,*) "Warning - setupcldtot: this code version is only designed for i_cloud_q_innovation == 2" + if (i_cloud_q_innovation < 20 .or. i_cloud_q_innovation > 22 ) then + write(*,*) "Warning - setupcldtot: this code version is only designed for i_cloud_q_innovation == 20,21,22" return else @@ -565,6 +565,8 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di ddiff=qv_ob-q_bk(k) q_build0_count=q_build0_count+1 endif + ! build error = 80% + error=one/(cloudqvis*8.E-01_r_kind) elseif (qob > -0.000001_r_single) then @@ -577,13 +579,15 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di ddiff=qv_ob-q_bk(k) q_clear0_count=q_clear0_count+1 endif + ! clear error = 30% + error=one/(cloudqvis*3.E-01_r_kind) else cycle endif q_obcount=q_obcount+1 - error=one/(cloudqvis*3.E-01_r_kind) + ! all obs errors = 30% ratio_errors=1.0_r_kind val = error*ddiff @@ -711,7 +715,7 @@ subroutine setupcldtot(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_di !! Write information to diagnostic file if(conv_diagsave)then - if (i_cloud_q_innovation == 2 .and. iip>0) then + if (i_cloud_q_innovation >= 20 .and. iip>0) then if(netcdf_diag) call nc_diag_write if(binary_diag)then write(7)' q',nchar,nreal,iip,mype,ioff0