Skip to content

Commit

Permalink
Refine use of PM2.5 for RRFS-SD;
Browse files Browse the repository at this point in the history
  • Loading branch information
@hongli-wang
hongli-wang committed Dec 12, 2022
1 parent cbe5180 commit dfef936
Show file tree
Hide file tree
Showing 7 changed files with 89 additions and 89 deletions.
5 changes: 4 additions & 1 deletion src/gsi/chemmod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
16 changes: 10 additions & 6 deletions src/gsi/gsi_rfv3io_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/gsi/gsimod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
72 changes: 33 additions & 39 deletions src/gsi/intpm2_5.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -425,9 +422,6 @@ subroutine intpm2_5_(pm2_5head,rval,sval)


end if
!

!

if (wrf_mass_regional .and. laeroana_gocart) then

Expand Down Expand Up @@ -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

Expand Down
2 changes: 1 addition & 1 deletion src/gsi/m_pm2_5Node.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
42 changes: 21 additions & 21 deletions src/gsi/setuppm2_5.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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),&
Expand Down Expand Up @@ -737,42 +740,39 @@ 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,&
mype,nfldsig)
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
Expand Down
Loading

0 comments on commit dfef936

Please sign in to comment.