diff --git a/src/gsi/chemmod.f90 b/src/gsi/chemmod.f90 index dc3f53a160..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 @@ -98,7 +100,7 @@ module chemmod code_pm25_anowbufr=102,code_pm10_ncbufr=code_pm25_ncbufr,& code_pm10_anowbufr=code_pm25_anowbufr - real(r_kind),parameter :: pm2_5_teom_max=100.0_r_kind !ug/m3 + real(r_kind),parameter :: pm2_5_teom_max=900.0_r_kind !ug/m3 !some parameters need to be put here since convinfo file won't !accomodate, stands for maximum realistic value of surface pm2.5 real(r_kind),parameter :: pm10_teom_max=150.0_r_kind !ug/m3 @@ -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 6f50d13861..eb7a86160f 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -21,7 +21,8 @@ module gsi_rfv3io_mod ! 2022-03-15 Hu - add code to read/write 2m T and Q for they will be ! used as background for surface observation operator ! 2022-04-15 Wang - add IO for regional FV3-CMAQ (RRFS-CMAQ) model - +! 2022-08-10 Wang - add IO for regional FV3-SMOKE (RRFS-SMOKE) model +! ! subroutines included: ! sub gsi_rfv3io_get_grid_specs ! sub read_fv3_files @@ -54,6 +55,8 @@ module gsi_rfv3io_mod use guess_grids, only: nfldsig,ntguessig,ifilesig use rapidrefresh_cldsurf_mod, only: i_use_2mq4b,i_use_2mt4b use chemmod, only: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3,laeroana_fv3cmaq + use chemmod, only: naero_smoke_fv3,aeronames_smoke_fv3,laeroana_fv3smoke + implicit none public type_fv3regfilenameg public bg_fv3regfilenameg @@ -87,22 +90,23 @@ module gsi_rfv3io_mod type(sub2grid_info) :: grd_fv3lam_dynvar_ionouv type(sub2grid_info) :: grd_fv3lam_tracer_ionouv type(sub2grid_info) :: grd_fv3lam_tracerchem_ionouv + type(sub2grid_info) :: grd_fv3lam_tracersmoke_ionouv type(sub2grid_info) :: grd_fv3lam_uv integer(i_kind) ,parameter:: ndynvarslist=13, ntracerslist=8 character(len=max_varname_length), dimension(ndynvarslist), parameter :: & vardynvars = [character(len=max_varname_length) :: & "u","v","u_w","u_s","v_w","v_s","t","tv","tsen","w","delp","ps","delzinc"] - character(len=max_varname_length), dimension(ntracerslist+naero_cmaq_fv3+7), parameter :: & + character(len=max_varname_length), dimension(ntracerslist+naero_cmaq_fv3+7+naero_smoke_fv3), parameter :: & vartracers = [character(len=max_varname_length) :: & - 'q','oz','ql','qi','qr','qs','qg','qnr',aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk'] - character(len=max_varname_length), dimension(15+naero_cmaq_fv3+7), parameter :: & + 'q','oz','ql','qi','qr','qs','qg','qnr',aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk',aeronames_smoke_fv3] + character(len=max_varname_length), dimension(15+naero_cmaq_fv3+7+naero_smoke_fv3), parameter :: & varfv3name = [character(len=max_varname_length) :: & 'u','v','W','T','delp','sphum','o3mr','liq_wat','ice_wat','rainwat','snowwat','graupel','rain_nc','ps','DZ', & - aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk'], & + aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk',aeronames_smoke_fv3], & vgsiname = [character(len=max_varname_length) :: & 'u','v','w','tsen','delp','q','oz','ql','qi','qr','qs','qg','qnr','ps','delzinc', & - aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk'] + aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk',aeronames_smoke_fv3] character(len=max_varname_length),dimension(:),allocatable:: name_metvars2d character(len=max_varname_length),dimension(:),allocatable:: name_metvars3d character(len=max_varname_length),dimension(:),allocatable:: name_chemvars3d @@ -125,7 +129,8 @@ module gsi_rfv3io_mod public :: k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m public :: ijns,ijns2d,displss,displss2d,ijnz,displsz_g - public :: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv,fv3lam_io_tracerchemvars3d_nouv + public :: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv + public :: fv3lam_io_tracerchemvars3d_nouv,fv3lam_io_tracersmokevars3d_nouv public :: fv3lam_io_dynmetvars2d_nouv,fv3lam_io_tracermetvars2d_nouv integer(i_kind) mype_u,mype_v,mype_t,mype_q,mype_p,mype_delz,mype_oz,mype_ql @@ -154,6 +159,7 @@ module gsi_rfv3io_mod character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_tracermetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_tracerchemvars3d_nouv + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_tracersmokevars3d_nouv ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_dynmetvars2d_nouv ! copy of cvars3d excluding uv 3-d fields @@ -166,7 +172,8 @@ module gsi_rfv3io_mod type(gsi_bundle):: gsibundle_fv3lam_dynvar_nouv type(gsi_bundle):: gsibundle_fv3lam_tracer_nouv type(gsi_bundle):: gsibundle_fv3lam_tracerchem_nouv - + type(gsi_bundle):: gsibundle_fv3lam_tracersmoke_nouv + contains subroutine fv3regfilename_init(this,it) implicit None @@ -873,6 +880,9 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) real(r_kind),dimension(:,:,:),pointer::ges_amassj=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_amassk=>NULL() + real(r_kind),dimension(:,:,:),pointer::ges_smoke=>NULL() + real(r_kind),dimension(:,:,:),pointer::ges_dust=>NULL() + character(len=max_varname_length) :: vartem="" character(len=64),dimension(:,:),allocatable:: names !to be same as in the grid the dummy sub2grid_info @@ -881,7 +891,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) integer(i_kind),dimension(:,:),allocatable:: uvlnames integer(i_kind):: inner_vars,numfields integer(i_kind):: ndynvario2d,ntracerio2d,ilev,jdynvar,jtracer - integer(r_kind):: iuv,ndynvario3d,ntracerio3d,ntracerchemio3d + integer(i_kind):: iuv,ndynvario3d,ntracerio3d + integer(i_kind):: ntracerchemio3d,ntracersmokeio3d integer(i_kind):: loc_id,ncfmt !clt this block is still maintained for they would be needed for a certain 2d fields IO @@ -915,7 +926,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) allocate( name_metvars3d(GSI_MetGuess_Bundle(it)%n3d)) end if - if (laeroana_fv3cmaq) then + if (laeroana_fv3cmaq.or.laeroana_fv3smoke) then if (.not.allocated(name_chemvars3d)) then allocate( name_chemvars3d(GSI_ChemGuess_Bundle(it)%n3d)) endif @@ -923,7 +934,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) call gsi_bundleinquire (GSI_MetGuess_Bundle(it),'shortnames::2d', name_metvars2d,istatus) call gsi_bundleinquire (GSI_MetGuess_Bundle(it),'shortnames::3d', name_metvars3d,istatus) - if (laeroana_fv3cmaq) then + if (laeroana_fv3cmaq.or.laeroana_fv3smoke) then call gsi_bundleinquire (GSI_ChemGuess_Bundle(it),'shortnames::3d', name_chemvars3d,istatus) endif @@ -934,7 +945,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) do i=1,GSI_MetGuess_Bundle(it)%n3d write(6,*)'metvardeb333-3d name ', trim(name_metvars3d(i)) enddo - if (laeroana_fv3cmaq) then + if (laeroana_fv3cmaq.or.laeroana_fv3smoke) then do i=1,GSI_ChemGuess_Bundle(it)%n3d write(6,*)'chemvardeb333-3d name ', trim(name_chemvars3d(i)) enddo @@ -980,6 +991,10 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if (laeroana_fv3cmaq) then allocate(fv3lam_io_tracerchemvars3d_nouv(naero_cmaq_fv3+7)) endif + + if (laeroana_fv3smoke) then + allocate(fv3lam_io_tracersmokevars3d_nouv(naero_smoke_fv3+1)) + endif jdynvar=0 jtracer=0 @@ -1101,6 +1116,36 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif !laeroana_fv3cmaq + if (laeroana_fv3smoke) then + jtracer = 0 +!name_chemvars3d chemguess from anainfo + do i=1,size(name_chemvars3d) + vartem=trim(name_chemvars3d(i)) + if (ifindstrloc(aeronames_smoke_fv3,trim(vartem)) > 0) then + 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 + 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 + fv3lam_io_tracersmokevars3d_nouv(jtracer+1)="pm2_5" + + if (mype == 0) then + write(6,*) ' fv3lam_io_tracersmokevars3d_nouv is',(trim(fv3lam_io_tracersmokevars3d_nouv(i)),i=1,ntracersmokeio3d) + endif + + endif !laeroana_fv3smoke + + if (allocated(fv3lam_io_dynmetvars2d_nouv) ) then call gsi_bundlecreate(gsibundle_fv3lam_dynvar_nouv,GSI_MetGuess_Bundle(it)%grid,'gsibundle_fv3lam_dynvar_nouv',istatus,& names2d=fv3lam_io_dynmetvars2d_nouv,names3d=fv3lam_io_dynmetvars3d_nouv) @@ -1127,6 +1172,14 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif endif + if (laeroana_fv3smoke) then + if (allocated(fv3lam_io_tracersmokevars3d_nouv) ) then + call gsi_bundlecreate(gsibundle_fv3lam_tracersmoke_nouv,GSI_ChemGuess_Bundle(it)%grid,'gsibundle_fv3lam_tracersmoke_nouv',istatus,& + names3d=fv3lam_io_tracersmokevars3d_nouv) + endif + endif + + inner_vars=1 numfields=inner_vars*(ndynvario3d*grd_a%nsig+ndynvario2d) allocate(lnames(1,numfields),names(1,numfields)) @@ -1183,6 +1236,25 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif + if (laeroana_fv3smoke) then + inner_vars=1 + numfields=inner_vars*(ntracersmokeio3d*grd_a%nsig) + deallocate(lnames,names) + allocate(lnames(1,numfields),names(1,numfields)) + ilev=1 + do i=1,ntracersmokeio3d + do k=1,grd_a%nsig + lnames(1,ilev)=k + names(1,ilev)=trim(fv3lam_io_tracersmokevars3d_nouv(i)) + ilev=ilev+1 + enddo + enddo + call general_sub2grid_create_info(grd_fv3lam_tracersmoke_ionouv,inner_vars,grd_a%nlat,& + grd_a%nlon,grd_a%nsig,numfields,regional,names=names,lnames=lnames) + + endif + + inner_vars=2 numfields=grd_a%nsig allocate(uvlnames(inner_vars,numfields),uvnames(inner_vars,numfields)) @@ -1254,6 +1326,12 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if (ier/=0) call die(trim(myname),'cannot get pointers for fv3chem-fields, ier =',ier) end if + if (laeroana_fv3smoke) then + call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it),'smoke',ges_smoke,istatus );ier=ier+istatus + call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it),'dust', ges_dust,istatus );ier=ier+istatus + call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it),'pm2_5',ges_pm2_5,istatus );ier=ier+istatus + end if + if( fv3sar_bg_opt == 0) then call gsi_fv3ncdf_readuv(grd_fv3lam_uv,ges_u,ges_v,fv3filenamegin(it)) else @@ -1269,6 +1347,10 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) call gsi_fv3ncdf_read(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv & & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) endif + if (laeroana_fv3smoke) then + call gsi_fv3ncdf_read(grd_fv3lam_tracersmoke_ionouv,gsibundle_fv3lam_tracersmoke_nouv & + & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) + endif else call gsi_fv3ncdf_read_v1(grd_fv3lam_dynvar_ionouv,gsibundle_fv3lam_dynvar_nouv & & ,fv3filenamegin(it)%dynvars,fv3filenamegin(it)) @@ -1278,6 +1360,10 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) call gsi_fv3ncdf_read_v1(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv & & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) endif + if (laeroana_fv3smoke) then + call gsi_fv3ncdf_read_v1(grd_fv3lam_tracersmoke_ionouv,gsibundle_fv3lam_tracersmoke_nouv & + & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) + endif endif if (laeroana_fv3cmaq) then @@ -1324,12 +1410,35 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) enddo end if ! laeroana_fv3cmaq + if (laeroana_fv3smoke) then + ier=0 + 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) + ! Calculate pm2_5 + do k=1,nsig + do j=1,lon2 + do i=1,lat2 + ges_pm2_5(i,j,k)=ges_smoke(i,j,k) + ges_dust(i,j,k) + enddo + enddo + enddo + endif !laeroana_fv3smoke + if( fv3sar_bg_opt == 0) then call GSI_BundleGetPointer ( gsibundle_fv3lam_dynvar_nouv, 'delp' ,ges_delp ,istatus );ier=ier+istatus if(istatus==0) ges_delp=ges_delp*0.001_r_kind endif call gsi_copy_bundle(gsibundle_fv3lam_dynvar_nouv,GSI_MetGuess_Bundle(it)) call gsi_copy_bundle(gsibundle_fv3lam_tracer_nouv,GSI_MetGuess_Bundle(it)) + if (laeroana_fv3cmaq) then + call gsi_copy_bundle(gsibundle_fv3lam_tracerchem_nouv,GSI_ChemGuess_Bundle(it)) + endif + if (laeroana_fv3smoke) then + call gsi_copy_bundle(gsibundle_fv3lam_tracersmoke_nouv,GSI_ChemGuess_Bundle(it)) + endif + call GSI_BundleGetPointer ( gsibundle_fv3lam_dynvar_nouv, 'tsen' ,ges_tsen_readin ,istatus );ier=ier+istatus !! tsen2tv !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do k=1,nsig @@ -1866,8 +1975,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m) end do end do - deallocate (dim_id,sfc,sfc1,dim) - deallocate (sfc_fulldomain) + if(allocated(sfc1) .and. allocated(sfc))deallocate (dim_id,sfc,sfc1,dim) + if(allocated(sfc_fulldomain)) deallocate (sfc_fulldomain) endif ! mype @@ -2825,6 +2934,10 @@ subroutine wrfv3_netcdf(fv3filenamegin) if (laeroana_fv3cmaq) then call gsi_copy_bundle(GSI_ChemGuess_Bundle(it),gsibundle_fv3lam_tracerchem_nouv) end if + if (laeroana_fv3smoke) then + call gsi_copy_bundle(GSI_ChemGuess_Bundle(it),gsibundle_fv3lam_tracersmoke_nouv) + end if + call gsi_bundleputvar (gsibundle_fv3lam_dynvar_nouv,'tsen',ges_tsen(:,:,:,it),istatus) if( fv3sar_bg_opt == 0) then @@ -2871,6 +2984,10 @@ subroutine wrfv3_netcdf(fv3filenamegin) call gsi_fv3ncdf_write(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv, & add_saved,fv3filenamegin%tracers,fv3filenamegin) endif + if (laeroana_fv3smoke) then + call gsi_fv3ncdf_write(grd_fv3lam_tracersmoke_ionouv,gsibundle_fv3lam_tracersmoke_nouv,& + add_saved,fv3filenamegin%tracers,fv3filenamegin) + endif else call gsi_fv3ncdf_write_v1(grd_fv3lam_dynvar_ionouv,gsibundle_fv3lam_dynvar_nouv,& @@ -2882,6 +2999,11 @@ subroutine wrfv3_netcdf(fv3filenamegin) call gsi_fv3ncdf_write_v1(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv,& add_saved,fv3filenamegin%tracers,fv3filenamegin) endif + if (laeroana_fv3smoke) then + call gsi_fv3ncdf_write_v1(grd_fv3lam_tracersmoke_ionouv,gsibundle_fv3lam_tracersmoke_nouv,& + add_saved,fv3filenamegin%tracers,fv3filenamegin) + endif + endif if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 8bc149a3a1..55d9298614 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -183,7 +183,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,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 @@ -1568,9 +1568,9 @@ module gsimod oneob_type_chem,oblat_chem,oblon_chem,obpres_chem,& diag_incr,elev_tolerance,tunable_error,& in_fname,out_fname,incr_fname,& - laeroana_gocart, laeroana_fv3cmaq,l_aoderr_table, aod_qa_limit, & + 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 305215389f..42c05ad3ff 100644 --- a/src/gsi/intpm2_5.f90 +++ b/src/gsi/intpm2_5.f90 @@ -78,6 +78,7 @@ subroutine intpm2_5_(pm2_5head,rval,sval) use gridmod, only: cmaq_regional,wrf_mass_regional,fv3_cmaq_regional use chemmod, only: s_2_5,d_2_5,nh4_mfac,oc_mfac,laeroana_gocart 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 implicit none @@ -89,10 +90,9 @@ subroutine intpm2_5_(pm2_5head,rval,sval) ! declare local variables integer(i_kind) j1,j2,j3,j4,j5,j6,j7,j8,ier,istatus real(r_kind) w1,w2,w3,w4,w5,w6,w7,w8 -! real(r_kind) penalty real(r_kind) cg_pm2_5,val,p0,grad,wnotgross,wgross,pm2_5_pg - real(r_kind),pointer,dimension(:) :: spm2_5 - real(r_kind),pointer,dimension(:) :: rpm2_5 + real(r_kind),pointer,dimension(:) :: spm2_5,sdust + real(r_kind),pointer,dimension(:) :: rpm2_5,rdust type(pm2_5Node), pointer :: pm2_5ptr character(len=max_varname_length) :: aeroname @@ -168,7 +168,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 @@ -226,7 +226,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) @@ -275,7 +275,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 @@ -307,6 +307,119 @@ subroutine intpm2_5_(pm2_5head,rval,sval) end if + if (laeroana_fv3smoke) then + !pm2_5ptr => pm2_5head + ier=0 + iaero=1 + aeroname=aeronames_smoke_fv3(iaero) !'smoke' + call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus);ier=istatus+ier + iaero=2 + aeroname=aeronames_smoke_fv3(iaero) !'dust' + call gsi_bundlegetpointer(sval,trim(aeroname),sdust,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,trim(aeroname),rdust,istatus);ier=istatus+ier + + if (istatus /= 0) then + write(6,*) 'error gsi_bundlegetpointer in intpm2_5 for ',aeronames_smoke_fv3(1:2) + call stop2(454) + endif + + pm2_5ptr => pm2_5Node_typecast(pm2_5head) + do while (associated(pm2_5ptr)) + j1=pm2_5ptr%ij(1) + j2=pm2_5ptr%ij(2) + j3=pm2_5ptr%ij(3) + j4=pm2_5ptr%ij(4) + j5=pm2_5ptr%ij(5) + j6=pm2_5ptr%ij(6) + j7=pm2_5ptr%ij(7) + j8=pm2_5ptr%ij(8) + w1=pm2_5ptr%wij(1) + w2=pm2_5ptr%wij(2) + w3=pm2_5ptr%wij(3) + w4=pm2_5ptr%wij(4) + w5=pm2_5ptr%wij(5) + w6=pm2_5ptr%wij(6) + w7=pm2_5ptr%wij(7) + w8=pm2_5ptr%wij(8) + + iaero=1 + val= (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))*pm2_5ptr%pm25wc(iaero) + iaero=2 + val= (w1*sdust(j1)+w2*sdust(j2)+ & + w3*sdust(j3)+w4*sdust(j4)+ & + w5*sdust(j5)+w6*sdust(j6)+ & + w7*sdust(j7)+w8*sdust(j8))*pm2_5ptr%pm25wc(iaero)+ val + + if(luse_obsdiag)then + if (lsaveobsens) then + call obsdiagNode_set(pm2_5ptr%diags,jiter=jiter,obssen=val*pm2_5ptr%raterr2*pm2_5ptr%err2) + else + if (pm2_5ptr%luse) call obsdiagNode_set(pm2_5ptr%diags,jiter=jiter,tldepart=val) + endif + endif + + if (l_do_adjoint) then + if (lsaveobsens) then + call obsdiagNode_get(pm2_5ptr%diags,jiter=jiter,obssen=grad) + + else + if( .not. ladtest_obs ) val=val-pm2_5ptr%res + +! gradient of nonlinear operator + + if (nlnqc_iter .and. pm2_5ptr%pg > tiny_r_kind .and. & + pm2_5ptr%b > tiny_r_kind) then + 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 + 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 + endif + + if ( ladtest_obs ) then + grad = val + else + grad = val*pm2_5ptr%raterr2*pm2_5ptr%err2 + end if + endif + +! adjoint + iaero=1 + 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) + iaero=2 + rdust(j1)=rdust(j1)+w1*grad*pm2_5ptr%pm25wc(iaero) + rdust(j2)=rdust(j2)+w2*grad*pm2_5ptr%pm25wc(iaero) + rdust(j3)=rdust(j3)+w3*grad*pm2_5ptr%pm25wc(iaero) + rdust(j4)=rdust(j4)+w4*grad*pm2_5ptr%pm25wc(iaero) + rdust(j5)=rdust(j5)+w5*grad*pm2_5ptr%pm25wc(iaero) + rdust(j6)=rdust(j6)+w6*grad*pm2_5ptr%pm25wc(iaero) + rdust(j7)=rdust(j7)+w7*grad*pm2_5ptr%pm25wc(iaero) + rdust(j8)=rdust(j8)+w8*grad*pm2_5ptr%pm25wc(iaero) + endif ! l_do_adjoint + + pm2_5ptr => pm2_5Node_nextcast(pm2_5ptr) + + end do + + nullify(spm2_5) + nullify(rpm2_5) + nullify(sdust) + nullify(rdust) + + end if ! laeroana_fv3smoke + if (wrf_mass_regional .and. laeroana_gocart) then !pm2_5ptr => pm2_5head @@ -513,15 +626,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 1035221786..79e6129cd8 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 @@ -110,7 +114,8 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) use chemmod, only : s_2_5,d_2_5,nh4_mfac,oc_mfac 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_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3,laeroana_fv3cmaq + 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 @@ -178,7 +183,7 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) integer(i_kind) :: ipm2_5,n_gocart_var - integer(i_kind) :: n_cmaq_var + integer(i_kind) :: n_cmaq_var,n_smoke_var real(r_kind),allocatable,dimension(:,:,:,:,:) :: pm25wc real(r_kind) :: pm25wc_ges(3) @@ -201,8 +206,7 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) !********************************************************************************* ! get pointer to pm2_5 guess state, if not present return - if ( (fv3_cmaq_regional .and. .not. laeroana_fv3cmaq) .or. cmaq_regional .or. (wrf_mass_regional .and. wrf_pm2_5) ) then -!for fv3cmaq, ges_pm2_5 is calculated in read_fv3 aeros + if ( cmaq_regional .or. (wrf_mass_regional .and. wrf_pm2_5) ) then call gsi_chemguess_get ('var::pm2_5', ipm2_5, ier ) if (ipm2_5 <= 0) then @@ -232,6 +236,72 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) endif + if (laeroana_fv3smoke)then +! 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) + 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 + call stop2(451) + endif + endif + + do i=1,naero_smoke_fv3 + aeroname=aeronames_smoke_fv3(i) + call gsi_chemguess_get ('var::'//trim(aeroname), ipm2_5, ier ) + if (ier > 0 .or. ipm2_5 <= 0) then + write(6,*) 'setuppm2_5: ',trim(aeroname),' missing in anavinfo' + call stop2(452) + endif + enddo + + if (size(gsi_chemguess_bundle)==nfldsig) then + aeroname='smoke' + call gsi_bundlegetpointer(gsi_chemguess_bundle(1),trim(aeroname),& + rank3,ier) + if (ier==0) then + allocate(ges_pm2_5(size(rank3,1),size(rank3,2),size(rank3,3),& + nfldsig)) + ges_pm2_5(:,:,:,1)=rank3 + allocate(pm25wc(size(rank3,1),size(rank3,2),size(rank3,3),naero_smoke_fv3,nfldsig)) + pm25wc(:,:,:,1,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_chemguess_bundle(ifld),trim(aeroname),rank3,ier) + ges_pm2_5(:,:,:,ifld)=rank3 + pm25wc(:,:,:,1,ifld)=rank3 + enddo + else + 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),& + rank3,ier) + pm25wc(:,:,:,i,1)=rank3 + if (ier==0) then + ges_pm2_5(:,:,:,1)=ges_pm2_5(:,:,:,1)+rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_chemguess_bundle(ifld),trim(aeroname),rank3,ier) + ges_pm2_5(:,:,:,ifld)=ges_pm2_5(:,:,:,ifld)+rank3 + pm25wc(:,:,:,i,ifld)=rank3 + enddo + else + write(6,*) 'setuppm2_5: ',trim(aeroname),' not found in chembundle,ier= ',ier + call stop2(453) + end if + end do + else + write(6,*) 'setuppm2_5: size(gsi_chemguess_bundle)/=nfldsig ges_pm2_5 not setup !!!' + call stop2(454) + end if ! eq. nfldsig + + endif + if (fv3_cmaq_regional .and. laeroana_fv3cmaq) then !check if pm25at, pm25ac and pm25co are in ges call gsi_chemguess_get ('var::pm25at', ipm2_5, ier ) @@ -632,7 +702,7 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) !convert for cmaq as well - if (wrf_mass_regional .or. fv3_cmaq_regional) then + if (wrf_mass_regional .or. fv3_cmaq_regional .or. laeroana_fv3smoke) then call tintrp2a11(ges_ps,ps_ges,dlat,dlon,dtime,hrdifsig,& mype,nfldsig) @@ -669,17 +739,48 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave) call tintrp2a11(ges_pm2_5,pm2_5ges,dlat,dlon,dtime,hrdifsig,& mype,nfldsig) innov = conc - pm2_5ges + if (laeroana_fv3smoke) then + 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 < 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,& + mype,nfldsig) + call tintrp2a11(pm25wc(:,:,:,2,nfldsig),pm25wc_ges(2),dlat,dlon,dtime,hrdifsig,& + mype,nfldsig) + if (pm25wc_ges(1) >= 1.0_r_kind) then + pm25wc_ges(1)=1.0_r_kind + else + pm25wc_ges(2)=0.0_r_kind + end if + if (pm25wc_ges(2) >= 1.0_r_kind) then + pm25wc_ges(2)=1.0_r_kind + else + pm25wc_ges(2)=0.0_r_kind + end if else pm25wc_ges = 0.0_r_kind end if diff --git a/src/gsi/stppm2_5.f90 b/src/gsi/stppm2_5.f90 index 87f0f9a4f5..3483563b9e 100644 --- a/src/gsi/stppm2_5.f90 +++ b/src/gsi/stppm2_5.f90 @@ -66,6 +66,7 @@ subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep) use gridmod, only: cmaq_regional,wrf_mass_regional,fv3_cmaq_regional USE chemmod, ONLY: d_2_5,s_2_5,nh4_mfac,oc_mfac,laeroana_gocart 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 implicit none @@ -82,7 +83,7 @@ subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep) real(r_kind) cg_pm2_5,val,val2,wgross,wnotgross,pm2_5_pg real(r_kind),dimension(max(1,nstep))::pen real(r_kind) w1,w2,w3,w4,w5,w6,w7,w8,qq - real(r_kind),pointer,dimension(:):: rpm2_5,spm2_5 + real(r_kind),pointer,dimension(:):: rpm2_5,spm2_5,rdust,sdust type(pm2_5Node), pointer :: pm2_5ptr character(len=max_varname_length) :: aeroname @@ -255,6 +256,105 @@ subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep) end if + if ( laeroana_fv3smoke) then + pm2_5ptr => pm2_5Node_typecast(pm2_5head) + + ier=0 + iaero=1 + aeroname=aeronames_smoke_fv3(iaero) !'smoke' + call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus);ier=istatus+ier + iaero=2 + aeroname=aeronames_smoke_fv3(iaero) !'dust' + call gsi_bundlegetpointer(sval,trim(aeroname),sdust,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,trim(aeroname),rdust,istatus);ier=istatus+ier + + if (istatus /= 0) then + write(6,*) 'error gsi_bundlegetpointer in intpm2_5 for ',aeronames_smoke_fv3(1:2) + call stop2(454) + endif + + do while (associated(pm2_5ptr)) + if(pm2_5ptr%luse)then + if(nstep > 0)then + j1=pm2_5ptr%ij(1) + j2=pm2_5ptr%ij(2) + j3=pm2_5ptr%ij(3) + j4=pm2_5ptr%ij(4) + j5=pm2_5ptr%ij(5) + j6=pm2_5ptr%ij(6) + j7=pm2_5ptr%ij(7) + j8=pm2_5ptr%ij(8) + w1=pm2_5ptr%wij(1) + w2=pm2_5ptr%wij(2) + w3=pm2_5ptr%wij(3) + w4=pm2_5ptr%wij(4) + w5=pm2_5ptr%wij(5) + w6=pm2_5ptr%wij(6) + w7=pm2_5ptr%wij(7) + w8=pm2_5ptr%wij(8) + + istatus=0 + + val2=-pm2_5ptr%res + val=zero + + iaero=1 + 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 + iaero=2 + val= pm2_5ptr%pm25wc(iaero)* & + (w1*rdust(j1)+w2*rdust(j2)+w3*rdust(j3)+w4*rdust(j4)+& + w5*rdust(j5)+w6*rdust(j6)+w7*rdust(j7)+w8*rdust(j8))+val + + iaero=1 + 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 + iaero=2 + val2= pm2_5ptr%pm25wc(iaero)* & + (w1*sdust(j1)+w2*sdust(j2)+w3*sdust(j3)+w4*sdust(j4)+& + w5*sdust(j5)+w6*sdust(j6)+w7*sdust(j7)+w8*sdust(j8))+val2 + + 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 !nstep + +! modify penalty term if nonlinear qc + + if (nlnqc_iter .and. pm2_5ptr%pg > tiny_r_kind .and. & + pm2_5ptr%b > tiny_r_kind) then + 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 + do kk=1,max(1,nstep) + pen(kk)= -two*log((exp(-half*pen(kk))+wgross)/(one+wgross)) + end do + endif + + out(1) = out(1)+pen(1)*pm2_5ptr%raterr2 + do kk=2,nstep + out(kk) = out(kk)+(pen(kk)-pen(1))*pm2_5ptr%raterr2 + end do + end if ! pm2_5ptr%luse + + pm2_5ptr => pm2_5Node_nextcast(pm2_5ptr) + + end do + + nullify(spm2_5) + nullify(rpm2_5) + nullify(sdust) + nullify(rdust) + + end if ! laeroana_fv3smoke + if (wrf_mass_regional .and. laeroana_gocart) then