From 5572c26cba2a773eefa91a555b28996c55b7e424 Mon Sep 17 00:00:00 2001 From: hongli-wang Date: Mon, 29 Aug 2022 15:42:17 +0000 Subject: [PATCH] 1. add copy ges for rrfs-cmaq 2. add pm2.5 DA for rrfs-smoke --- src/gsi/chemmod.f90 | 2 +- src/gsi/gsi_rfv3io_mod.f90 | 146 +++++++++++++++++++++++++++++++++---- src/gsi/gsimod.F90 | 4 +- src/gsi/intpm2_5.f90 | 120 ++++++++++++++++++++++++++++++ src/gsi/setuppm2_5.f90 | 76 ++++++++++++++++++- src/gsi/stppm2_5.f90 | 92 +++++++++++++++++++++++ 6 files changed, 419 insertions(+), 21 deletions(-) diff --git a/src/gsi/chemmod.f90 b/src/gsi/chemmod.f90 index dc3f53a160..aa3c90a448 100644 --- a/src/gsi/chemmod.f90 +++ b/src/gsi/chemmod.f90 @@ -98,7 +98,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 diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 6f50d13861..47961cccec 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(r_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,31 @@ 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 ' + else + write(6,*)'the chemvarname ',vartem,' is not in aeronames_smoke_fv3, !!!!!!!!!!' + !call flush(6) + !call stop2(333) + 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 +1167,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 +1231,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 +1321,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 +1342,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 +1355,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 +1405,36 @@ 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) + if (ier/=0) write(6,*),"cannot get pointers for fv3 smoke" +!! 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 +1971,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 +2930,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 +2980,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 +2995,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 eb10a7bf16..2a6d0e8136 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,crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, & + laeroana_fv3cmaq,laeroana_fv3smoke,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 @@ -1538,7 +1538,7 @@ 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, & raod_radius_mean_scale,raod_radius_std_scale, luse_deepblue,& diff --git a/src/gsi/intpm2_5.f90 b/src/gsi/intpm2_5.f90 index 305215389f..ecf8e7ccff 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 @@ -307,6 +308,125 @@ subroutine intpm2_5_(pm2_5head,rval,sval) end if +! + if (laeroana_fv3smoke) then + !pm2_5ptr => pm2_5head + 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) +!naero_smoke_fv3 + iaero=1 + aeroname=aeronames_smoke_fv3(iaero) !'smoke' + call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus) + if(istatus /= 0) then + write(6,*) 'error gsi_bundlegetpointer in intpm2_5 for ', aeroname + call stop2(454) + endif + + 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) + 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 + write(6,*) 'error gsi_bundlegetpointer in intpm2_5 for ',aeroname + call stop2(454) + endif + 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)) + val + + nullify(spm2_5) + end do !iaero + + + 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 + !aeroname='smoke' + do iaero=1,1 ! naero_smoke_fv3 ! Smoke only + 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 + rpm2_5(j2)=rpm2_5(j2)+w2*grad + rpm2_5(j3)=rpm2_5(j3)+w3*grad + rpm2_5(j4)=rpm2_5(j4)+w4*grad + rpm2_5(j5)=rpm2_5(j5)+w5*grad + rpm2_5(j6)=rpm2_5(j6)+w6*grad + rpm2_5(j7)=rpm2_5(j7)+w7*grad + rpm2_5(j8)=rpm2_5(j8)+w8*grad + nullify(rpm2_5) + end do + endif + + pm2_5ptr => pm2_5Node_nextcast(pm2_5ptr) + + end do + + + end if +! + +! + if (wrf_mass_regional .and. laeroana_gocart) then !pm2_5ptr => pm2_5head diff --git a/src/gsi/setuppm2_5.f90 b/src/gsi/setuppm2_5.f90 index 1035221786..d2773d0d62 100644 --- a/src/gsi/setuppm2_5.f90 +++ b/src/gsi/setuppm2_5.f90 @@ -110,7 +110,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 use gridmod, only : cmaq_regional,wrf_mass_regional,fv3_cmaq_regional implicit none @@ -178,7 +179,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 +202,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 +232,67 @@ 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 (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 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_chemguess_bundle(ifld),trim(aeroname),rank3,ier) + ges_pm2_5(:,:,:,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) + 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 + 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 ) @@ -669,6 +730,13 @@ 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 (abs(innov) >= 30.0 .or. conc >= 60.0) then + innov = innov + else + innov = 0.0 + end if + end if end if if ( fv3_cmaq_regional .and. laeroana_fv3cmaq) then diff --git a/src/gsi/stppm2_5.f90 b/src/gsi/stppm2_5.f90 index 87f0f9a4f5..15ca36f239 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 @@ -255,6 +256,97 @@ subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep) end if + if ( laeroana_fv3smoke) then + pm2_5ptr => pm2_5Node_typecast(pm2_5head) + + 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 + + !!! + 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 gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus) + + if(istatus /= 0) then + write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',& + aeroname + call stop2(457) + endif + + val= 1.0* & + (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= 1.0* & + (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 + 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 + +! 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 => pm2_5Node_nextcast(pm2_5ptr) + + end do + + + + end if ! laeroana_fv3smoke + if (wrf_mass_regional .and. laeroana_gocart) then