From dfef936367734a8d4c406463a8a05c00c79bce57 Mon Sep 17 00:00:00 2001 From: hongli-wang Date: Tue, 18 Oct 2022 18:58:40 +0000 Subject: [PATCH] Refine use of PM2.5 for RRFS-SD; --- src/gsi/chemmod.f90 | 5 ++- src/gsi/gsi_rfv3io_mod.f90 | 16 +++++---- src/gsi/gsimod.F90 | 4 +-- src/gsi/intpm2_5.f90 | 72 +++++++++++++++++--------------------- src/gsi/m_pm2_5Node.F90 | 2 +- src/gsi/setuppm2_5.f90 | 42 +++++++++++----------- src/gsi/stppm2_5.f90 | 37 ++++++++++---------- 7 files changed, 89 insertions(+), 89 deletions(-) diff --git a/src/gsi/chemmod.f90 b/src/gsi/chemmod.f90 index aa3c90a448..14a90c818c 100644 --- a/src/gsi/chemmod.f90 +++ b/src/gsi/chemmod.f90 @@ -40,7 +40,7 @@ module chemmod public :: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3 ! fv3smoke - public :: naero_smoke_fv3,aeronames_smoke_fv3 + public :: naero_smoke_fv3,aeronames_smoke_fv3,pm2_5_innov_threshold public :: naero_gocart_wrf,aeronames_gocart_wrf @@ -79,6 +79,8 @@ module chemmod integer(i_kind) :: icvt_cmaq_fv3 real(r_kind) :: raod_radius_mean_scale,raod_radius_std_scale real(r_kind) :: ppmv_conv = 96.06_r_kind/28.964_r_kind*1.0e+3_r_kind + real(r_kind) :: pm2_5_innov_threshold + logical :: wrf_pm2_5 logical :: lread_ext_aerosol ! if true, will read in aerosols from aerfXX rather than from sigfXX @@ -305,6 +307,7 @@ subroutine init_chem laeroana_gocart = .false. laeroana_fv3cmaq = .false. ! .true. for performing aerosol analysis for regional FV3-CMAQ model(Please other parameters requred in gsimod.F90) laeroana_fv3smoke = .false. + pm2_5_innov_threshold = 20.0_r_kind l_aoderr_table = .false. icvt_cmaq_fv3 = 1 ! 1. Control variable is individual aerosol specie; 2: CV is total mass per I,J,K mode raod_radius_mean_scale = 1.0_r_kind ! Tune radius of particles when calculating AOD using CRTM diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 47961cccec..27efb265ca 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -1125,10 +1125,15 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) jtracer=jtracer+1 fv3lam_io_tracersmokevars3d_nouv(jtracer)=trim(vartem) write(6,*)'the chemvarname ',jtracer,vartem,' is found ' + if (trim(vartem) == "pm2_5")then + write(6,*)'the chemvarname ',vartem,' will be calculated.' + endif + else - write(6,*)'the chemvarname ',vartem,' is not in aeronames_smoke_fv3, !!!!!!!!!!' - !call flush(6) - !call stop2(333) + if (trim(vartem) /= "pm2_5")then + write(6,*)'the chemvarname ',vartem,' is not in aeronames_smoke_fv3 !!!' + call flush(6) + endif endif enddo ntracersmokeio3d=jtracer+1 @@ -1410,9 +1415,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) call GSI_BundleGetPointer ( gsibundle_fv3lam_tracersmoke_nouv,'smoke',ges_smoke,istatus );ier=ier+istatus call GSI_BundleGetPointer ( gsibundle_fv3lam_tracersmoke_nouv,'dust',ges_dust,istatus );ier=ier+istatus call GSI_BundleGetPointer ( gsibundle_fv3lam_tracersmoke_nouv,'pm2_5',ges_pm2_5,istatus );ier=ier+istatus - ! if (ier/=0) call die(trim(myname),'cannot get pointers for fv3 ! chem-fields, ier =',ier) - if (ier/=0) write(6,*),"cannot get pointers for fv3 smoke" -!! pm2_5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (ier/=0) call die(trim(myname),'cannot get pointers for fv3 chem-fields, ier =',ier) + ! Calculate pm2_5 do k=1,nsig do j=1,lon2 do i=1,lat2 diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 2a6d0e8136..e919342556 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -180,7 +180,7 @@ module gsimod oblon_chem,obpres_chem,diag_incr,elev_tolerance,tunable_error,& in_fname,out_fname,incr_fname, & laeroana_gocart, l_aoderr_table, aod_qa_limit, luse_deepblue, lread_ext_aerosol, & - laeroana_fv3cmaq,laeroana_fv3smoke,crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, & + laeroana_fv3cmaq,laeroana_fv3smoke,pm2_5_innov_threshold,crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, & icvt_cmaq_fv3, raod_radius_mean_scale,raod_radius_std_scale use chemmod, only : wrf_pm2_5,aero_ratios @@ -1540,7 +1540,7 @@ module gsimod in_fname,out_fname,incr_fname,& laeroana_gocart, laeroana_fv3cmaq,laeroana_fv3smoke,l_aoderr_table, aod_qa_limit, & crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, & - icvt_cmaq_fv3, & + icvt_cmaq_fv3,pm2_5_innov_threshold, & raod_radius_mean_scale,raod_radius_std_scale, luse_deepblue,& aero_ratios,wrf_pm2_5, lread_ext_aerosol diff --git a/src/gsi/intpm2_5.f90 b/src/gsi/intpm2_5.f90 index 2272fc1843..85ea823299 100644 --- a/src/gsi/intpm2_5.f90 +++ b/src/gsi/intpm2_5.f90 @@ -169,7 +169,7 @@ subroutine intpm2_5_(pm2_5head,rval,sval) if( ladtest_obs ) then grad = val else - grad = val*pm2_5ptr%raterr2*pm2_5ptr%err2 + grad = val*pm2_5ptr%raterr2*pm2_5ptr%err2 end if endif ! adjoint @@ -227,7 +227,7 @@ subroutine intpm2_5_(pm2_5head,rval,sval) w7* spm2_5(j7)+w8* spm2_5(j8) val=val *pm2_5ptr%pm25wc(imodes_cmaq_fv3(iaero)) nullify(spm2_5) -! + do iaero=2, naero_cmaq_fv3 aeroname=aeronames_cmaq_fv3(iaero) call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus) @@ -276,7 +276,7 @@ subroutine intpm2_5_(pm2_5head,rval,sval) if( ladtest_obs ) then grad = val else - grad = val*pm2_5ptr%raterr2*pm2_5ptr%err2 + grad = val*pm2_5ptr%raterr2*pm2_5ptr%err2 end if endif @@ -308,7 +308,6 @@ subroutine intpm2_5_(pm2_5head,rval,sval) end if -! if (laeroana_fv3smoke) then !pm2_5ptr => pm2_5head pm2_5ptr => pm2_5Node_typecast(pm2_5head) @@ -329,11 +328,11 @@ subroutine intpm2_5_(pm2_5head,rval,sval) w6=pm2_5ptr%wij(6) w7=pm2_5ptr%wij(7) w8=pm2_5ptr%wij(8) -!naero_smoke_fv3 + iaero=1 aeroname=aeronames_smoke_fv3(iaero) !'smoke' call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus) - if(istatus /= 0) then + if (istatus /= 0) then write(6,*) 'error gsi_bundlegetpointer in intpm2_5 for ', aeroname call stop2(454) endif @@ -345,11 +344,11 @@ subroutine intpm2_5_(pm2_5head,rval,sval) val= val*pm2_5ptr%pm25wc(iaero) nullify(spm2_5) -! + do iaero=2, naero_smoke_fv3 aeroname=aeronames_smoke_fv3(iaero) call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus) - if(istatus /= 0) then + if (istatus /= 0) then write(6,*) 'error gsi_bundlegetpointer in intpm2_5 for ',aeroname call stop2(454) endif @@ -384,39 +383,37 @@ subroutine intpm2_5_(pm2_5head,rval,sval) pm2_5_pg=pm2_5ptr%pg*varqc_iter cg_pm2_5=cg_term/pm2_5ptr%b wnotgross= one-pm2_5_pg - wgross =pm2_5_pg*cg_pm2_5/wnotgross ! wgross isi gama in the reference by enderson + wgross=pm2_5_pg*cg_pm2_5/wnotgross ! wgross isi gama in the reference by enderson p0=wgross/(wgross+exp(-half*pm2_5ptr%err2*val**2)) ! p0 is p in the reference by enderson - val=val*(one-p0) ! term is wqc in the referenc by enderson + val=val*(one-p0) ! term is wqc in the referenc by enderson endif - if( ladtest_obs ) then + if ( ladtest_obs ) then grad = val else - grad = val*pm2_5ptr%raterr2*pm2_5ptr%err2 + grad = val*pm2_5ptr%raterr2*pm2_5ptr%err2 end if endif -! adjoint - !aeroname='smoke' - do iaero=1, naero_smoke_fv3 - aeroname = aeronames_smoke_fv3(iaero) - call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus) - if(istatus /= 0) then - write(6,*) 'error gsi_bundlegetpointer in intpm2_5 for ',& - aeroname - call stop2(455) - endif - - rpm2_5(j1)=rpm2_5(j1)+w1*grad*pm2_5ptr%pm25wc(iaero) - rpm2_5(j2)=rpm2_5(j2)+w2*grad*pm2_5ptr%pm25wc(iaero) - rpm2_5(j3)=rpm2_5(j3)+w3*grad*pm2_5ptr%pm25wc(iaero) - rpm2_5(j4)=rpm2_5(j4)+w4*grad*pm2_5ptr%pm25wc(iaero) - rpm2_5(j5)=rpm2_5(j5)+w5*grad*pm2_5ptr%pm25wc(iaero) - rpm2_5(j6)=rpm2_5(j6)+w6*grad*pm2_5ptr%pm25wc(iaero) - rpm2_5(j7)=rpm2_5(j7)+w7*grad*pm2_5ptr%pm25wc(iaero) - rpm2_5(j8)=rpm2_5(j8)+w8*grad*pm2_5ptr%pm25wc(iaero) - nullify(rpm2_5) - end do +! adjoint + do iaero=1, naero_smoke_fv3 + aeroname = aeronames_smoke_fv3(iaero) + call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus) + if (istatus /= 0) then + write(6,*) 'error gsi_bundlegetpointer in intpm2_5 for ',aeroname + call stop2(455) + endif + + rpm2_5(j1)=rpm2_5(j1)+w1*grad*pm2_5ptr%pm25wc(iaero) + rpm2_5(j2)=rpm2_5(j2)+w2*grad*pm2_5ptr%pm25wc(iaero) + rpm2_5(j3)=rpm2_5(j3)+w3*grad*pm2_5ptr%pm25wc(iaero) + rpm2_5(j4)=rpm2_5(j4)+w4*grad*pm2_5ptr%pm25wc(iaero) + rpm2_5(j5)=rpm2_5(j5)+w5*grad*pm2_5ptr%pm25wc(iaero) + rpm2_5(j6)=rpm2_5(j6)+w6*grad*pm2_5ptr%pm25wc(iaero) + rpm2_5(j7)=rpm2_5(j7)+w7*grad*pm2_5ptr%pm25wc(iaero) + rpm2_5(j8)=rpm2_5(j8)+w8*grad*pm2_5ptr%pm25wc(iaero) + nullify(rpm2_5) + end do ! iaero endif pm2_5ptr => pm2_5Node_nextcast(pm2_5ptr) @@ -425,9 +422,6 @@ subroutine intpm2_5_(pm2_5head,rval,sval) end if -! - -! if (wrf_mass_regional .and. laeroana_gocart) then @@ -635,15 +629,15 @@ subroutine intpm2_5_(pm2_5head,rval,sval) pm2_5_pg=pm2_5ptr%pg*varqc_iter cg_pm2_5=cg_term/pm2_5ptr%b wnotgross= one-pm2_5_pg - wgross =pm2_5_pg*cg_pm2_5/wnotgross ! wgross is gama in the reference by enderson + wgross=pm2_5_pg*cg_pm2_5/wnotgross ! wgross is gama in the reference by enderson p0=wgross/(wgross+exp(-half*pm2_5ptr%err2*val**2)) ! p0 is p in the reference by enderson - val=val*(one-p0) ! term is wqc in the referenc by enderson + val=val*(one-p0) ! term is wqc in the referenc by enderson endif if( ladtest_obs ) then grad = val else - grad = val*pm2_5ptr%raterr2*pm2_5ptr%err2 + grad = val*pm2_5ptr%raterr2*pm2_5ptr%err2 end if endif diff --git a/src/gsi/m_pm2_5Node.F90 b/src/gsi/m_pm2_5Node.F90 index 4531ec4433..715d7aacdd 100644 --- a/src/gsi/m_pm2_5Node.F90 +++ b/src/gsi/m_pm2_5Node.F90 @@ -51,7 +51,7 @@ module m_pm2_5Node real(r_kind) :: pg ! variational quality control parameter real(r_kind) :: wij(8) ! horizontal interpolation weights integer(i_kind) :: ij(8) ! horizontal locations - real(r_kind) :: pm25wc(3) ! weight coes at i,j,k modes + real(r_kind) :: pm25wc(3) ! weights !logical :: luse ! flag indicating if ob is used in pen. !integer(i_kind) :: idv,iob ! device id and obs index for sorting !real (r_kind) :: elat, elon ! earth lat-lon for redistribution diff --git a/src/gsi/setuppm2_5.f90 b/src/gsi/setuppm2_5.f90 index 600eb44226..ec390014c2 100644 --- a/src/gsi/setuppm2_5.f90 +++ b/src/gsi/setuppm2_5.f90 @@ -38,6 +38,10 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) ! 2017-02-09 guo - Remove m_alloc, n_alloc. ! . Remove my_node with corrected typecast(). ! 2022-04-19 h.wang - add code for fv3_cmaq_regional +! - fv3_cmaq_regional=.true. : model is regional FV3-CMAQ +! - laeroana_fv3cmaq=.true. : produce the analysis for regional FV3-CMAQ +! 2022-08-10 h.Wang - add code for regional FV3-SD (RRFS-SMOKE/DUST) model +! - laeroana_fv3smoke=.true. : produce the analysis for RRFS-SD ! ! input argument list: ! lunin - unit from which to read observations @@ -111,7 +115,7 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) use chemmod, only: naero_gocart_wrf,aeronames_gocart_wrf,& upper2lower,lower2upper,laeroana_gocart,wrf_pm2_5 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 chemmod, only: naero_smoke_fv3,aeronames_smoke_fv3,laeroana_fv3smoke,pm2_5_innov_threshold use gridmod, only : cmaq_regional,wrf_mass_regional,fv3_cmaq_regional implicit none @@ -233,12 +237,11 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) endif if (laeroana_fv3smoke)then -!check if aerosol species in control +! check if aerosol species in control call gsi_chemguess_get ( 'aerosols::3d', n_smoke_var, ier ) -!n_smoke_var ges vars in anainfo; naero_smoke_fv3 in chemmod -!if naero_smoke_fv3 is greater than 3, change the dimension of -!pm25wc_ges(naero_smoke_fv3) +! n_smoke_var ges vars in anainfo; naero_smoke_fv3 in chemmod +! if naero_smoke_fv3 is greater than 3, change the dimension of pm25wc_ges(naero_smoke_fv3) if (n_smoke_var /= naero_smoke_fv3) then if (n_smoke_var < naero_smoke_fv3) then write(6,*) 'setuppm2_5: not all smoke aerosols in anavinfo',n_smoke_var,naero_smoke_fv3 @@ -274,7 +277,7 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) write(6,*) 'setuppm2_5: ',trim(aeroname),' not found in chembundle,ier= ',ier call stop2(453) endif - !!! + do i=2,naero_smoke_fv3 aeroname=trim(aeronames_smoke_fv3(i)) call gsi_bundlegetpointer(gsi_chemguess_bundle(1),trim(aeroname),& @@ -737,23 +740,25 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) mype,nfldsig) innov = conc - pm2_5ges if (laeroana_fv3smoke) then - if ( -1.0*innov >= 5.0_r_kind .or. & - (innov > 15.0_r_kind .and. pm2_5ges >=1.0_r_kind).or. & - (conc >= 60.0_r_kind .and. pm2_5ges >=1.0_r_kind).or. & + if ( -1.0*innov >= pm2_5_innov_threshold .or. & + (innov > pm2_5_innov_threshold .and. pm2_5ges >=1.0_r_kind).or. & + (conc >= 40.0_r_kind .and. pm2_5ges >=1.0_r_kind).or. & conc >= 100.0_r_kind ) then innov = innov else innov = 0.0_r_kind + muse(i)=.false. end if - if (tv_ges-273.15_r_kind < 0.0_r_kind) then + if (tv_ges-273.15_r_kind < 5.0_r_kind) then innov = 0.0_r_kind + muse(i)=.false. end if end if end if if ( fv3_cmaq_regional .and. laeroana_fv3cmaq) then - ! interpoloate pm25ac + ! interpoloate pm25ac call tintrp2a11(pm25wc(:,:,:,1,nfldsig),pm25wc_ges(1),dlat,dlon,dtime,hrdifsig,& mype,nfldsig) call tintrp2a11(pm25wc(:,:,:,2,nfldsig),pm25wc_ges(2),dlat,dlon,dtime,hrdifsig,& @@ -761,18 +766,13 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) call tintrp2a11(pm25wc(:,:,:,3,nfldsig),pm25wc_ges(3),dlat,dlon,dtime,hrdifsig,& mype,nfldsig) elseif (laeroana_fv3smoke) then - call tintrp2a11(pm25wc(:,:,:,1,nfldsig),pm25wc_ges(1),dlat,dlon,dtime,hrdifsig,& + call tintrp2a11(pm25wc(:,:,:,1,nfldsig),pm25wc_ges(1),dlat,dlon,dtime,hrdifsig,& mype,nfldsig) - call tintrp2a11(pm25wc(:,:,:,2,nfldsig),pm25wc_ges(2),dlat,dlon,dtime,hrdifsig,& + call tintrp2a11(pm25wc(:,:,:,2,nfldsig),pm25wc_ges(2),dlat,dlon,dtime,hrdifsig,& mype,nfldsig) - print*,"tv_ges= ",tv_ges,conc,pm2_5ges,pm25wc_ges(1:2) - if ( sum(pm25wc_ges(1:2)) >= 1.0_r_kind) then - pm25wc_ges(1)=pm25wc_ges(1)/sum(pm25wc_ges(1:2)) - pm25wc_ges(2)=1.0_r_kind-pm25wc_ges(1) - else - pm25wc_ges(1)=1.0_r_kind - pm25wc_ges(2)=0.0_r_kind - end if + pm25wc_ges = 0.0_r_kind + if (pm25wc_ges(1) >= 1.0_r_kind) pm25wc_ges(1)=1.0_r_kind + if (pm25wc_ges(2) >= 1.0_r_kind) pm25wc_ges(2)=1.0_r_kind else pm25wc_ges = 0.0_r_kind end if diff --git a/src/gsi/stppm2_5.f90 b/src/gsi/stppm2_5.f90 index 565e1c0e1d..cd88bfd5f4 100644 --- a/src/gsi/stppm2_5.f90 +++ b/src/gsi/stppm2_5.f90 @@ -284,41 +284,40 @@ subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep) val2=-pm2_5ptr%res val=zero - !!! - do iaero=1,naero_smoke_fv3 - aeroname=aeronames_smoke_fv3(iaero) - call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus) - if(istatus /= 0) then - write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',& + do iaero=1,naero_smoke_fv3 + aeroname=aeronames_smoke_fv3(iaero) + call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus) + if (istatus /= 0) then + write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',& aeroname - call stop2(456) - endif + call stop2(456) + endif - call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus) + call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus) - if(istatus /= 0) then - write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',& + if (istatus /= 0) then + write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',& aeroname - call stop2(457) - endif + call stop2(457) + endif - val= pm2_5ptr%pm25wc(iaero)* & + val= pm2_5ptr%pm25wc(iaero)* & (w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4*rpm2_5(j4)+& w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7*rpm2_5(j7)+w8*rpm2_5(j8))+val - val2= pm2_5ptr%pm25wc(iaero)* & + val2= pm2_5ptr%pm25wc(iaero)* & (w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4*spm2_5(j4)+& w5* spm2_5(j5)+w6* spm2_5(j6)+w7*spm2_5(j7)+w8*spm2_5(j8))+val2 - nullify(spm2_5,rpm2_5) - end do + nullify(spm2_5) + end do ! iaero do kk=1,nstep qq=val2+sges(kk)*val pen(kk)=qq*qq*pm2_5ptr%err2 end do else pen(1)=pm2_5ptr%res*pm2_5ptr%res*pm2_5ptr%err2 - end if + end if !nstep ! modify penalty term if nonlinear qc @@ -337,7 +336,7 @@ subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep) do kk=2,nstep out(kk) = out(kk)+(pen(kk)-pen(1))*pm2_5ptr%raterr2 end do - end if + end if ! pm2_5ptr%luse pm2_5ptr => pm2_5Node_nextcast(pm2_5ptr)