From 58d2180a6d66fbe89e20e0b640c028f8723d9a6c Mon Sep 17 00:00:00 2001
From: hongli-wang <hongli.wang@noaa.gov>
Date: Mon, 29 Aug 2022 15:42:17 +0000
Subject: [PATCH 1/8] 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 dc3f53a16..aa3c90a44 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 6f50d1386..47961ccce 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 859d4cedb..13b2c2f44 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,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,7 +1568,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 305215389..ecf8e7ccf 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 103522178..d2773d0d6 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 87f0f9a4f..15ca36f23 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
 

From e8d086e7a4168912cc8730a6f36c82be617842ee Mon Sep 17 00:00:00 2001
From: hongli-wang <hongli.wang@noaa.gov>
Date: Fri, 14 Oct 2022 21:37:24 +0000
Subject: [PATCH 2/8] Refine PM2.5 DA for RRFS-SMOKE

---
 src/gsi/intpm2_5.f90   | 22 ++++++++++++----------
 src/gsi/setuppm2_5.f90 | 37 +++++++++++++++++++++++++++++++------
 src/gsi/stppm2_5.f90   |  4 ++--
 3 files changed, 45 insertions(+), 18 deletions(-)

diff --git a/src/gsi/intpm2_5.f90 b/src/gsi/intpm2_5.f90
index ecf8e7ccf..2272fc184 100644
--- a/src/gsi/intpm2_5.f90
+++ b/src/gsi/intpm2_5.f90
@@ -342,6 +342,8 @@ subroutine intpm2_5_(pm2_5head,rval,sval)
                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= val*pm2_5ptr%pm25wc(iaero)
+
           nullify(spm2_5)
 !
           do iaero=2, naero_smoke_fv3
@@ -354,7 +356,7 @@ subroutine intpm2_5_(pm2_5head,rval,sval)
           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
+               w7* spm2_5(j7)+w8* spm2_5(j8))*pm2_5ptr%pm25wc(iaero) + val
 
           nullify(spm2_5)
           end do !iaero
@@ -396,7 +398,7 @@ subroutine intpm2_5_(pm2_5head,rval,sval)
 
 !       adjoint
              !aeroname='smoke'
-          do iaero=1,1 ! naero_smoke_fv3  ! Smoke only 
+          do iaero=1, naero_smoke_fv3   
              aeroname = aeronames_smoke_fv3(iaero)
              call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)
              if(istatus /= 0) then
@@ -405,14 +407,14 @@ subroutine intpm2_5_(pm2_5head,rval,sval)
                 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
+             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
           endif
diff --git a/src/gsi/setuppm2_5.f90 b/src/gsi/setuppm2_5.f90
index d2773d0d6..600eb4422 100644
--- a/src/gsi/setuppm2_5.f90
+++ b/src/gsi/setuppm2_5.f90
@@ -237,7 +237,8 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave)
      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
@@ -262,9 +263,12 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave)
            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
@@ -275,11 +279,13 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave)
            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
@@ -693,7 +699,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)
 
@@ -731,23 +737,42 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave)
                 mype,nfldsig)
            innov = conc - pm2_5ges
            if (laeroana_fv3smoke) then
-              if (abs(innov) >= 30.0 .or. conc >= 60.0) 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. &
+                   conc >= 100.0_r_kind ) then
                  innov = innov
               else 
-                 innov = 0.0 
+                 innov = 0.0_r_kind 
+              end if
+              if (tv_ges-273.15_r_kind < 0.0_r_kind) then
+                 innov = 0.0_r_kind
               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)
+              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
         else
           pm25wc_ges = 0.0_r_kind
         end if
diff --git a/src/gsi/stppm2_5.f90 b/src/gsi/stppm2_5.f90
index 15ca36f23..565e1c0e1 100644
--- a/src/gsi/stppm2_5.f90
+++ b/src/gsi/stppm2_5.f90
@@ -302,10 +302,10 @@ subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep)
                    call stop2(457)
                 endif
 
-                val= 1.0* &
+                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= 1.0* &
+                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
 

From 274e1be7bc3f999f6586a21193a62469923dbb36 Mon Sep 17 00:00:00 2001
From: hongli-wang <hongli.wang@noaa.gov>
Date: Tue, 18 Oct 2022 18:58:40 +0000
Subject: [PATCH 3/8] 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 aa3c90a44..14a90c818 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 47961ccce..27efb265c 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 13b2c2f44..ee4ca4fd9 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,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
@@ -1570,7 +1570,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 2272fc184..85ea82329 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 4531ec443..715d7aacd 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 600eb4422..ec390014c 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 565e1c0e1..cd88bfd5f 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)
 

From 08196ee7322f5697ac5f5b64e7752549f416be4a Mon Sep 17 00:00:00 2001
From: hongli-wang <hongli.wang@noaa.gov>
Date: Fri, 16 Dec 2022 20:51:24 +0000
Subject: [PATCH 4/8] Move getpointer out of obs do loop for rrfs_sd pm2.5
 assimilation

---
 src/gsi/intpm2_5.f90   | 102 ++++++++++++++++++++---------------------
 src/gsi/setuppm2_5.f90 |  14 ++++--
 src/gsi/stppm2_5.f90   |  67 +++++++++++++++------------
 3 files changed, 99 insertions(+), 84 deletions(-)

diff --git a/src/gsi/intpm2_5.f90 b/src/gsi/intpm2_5.f90
index 85ea82329..d1f3b87c6 100644
--- a/src/gsi/intpm2_5.f90
+++ b/src/gsi/intpm2_5.f90
@@ -92,8 +92,8 @@ subroutine intpm2_5_(pm2_5head,rval,sval)
     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
@@ -310,6 +310,21 @@ subroutine intpm2_5_(pm2_5head,rval,sval)
 
    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)
@@ -330,36 +345,15 @@ subroutine intpm2_5_(pm2_5head,rval,sval)
           w8=pm2_5ptr%wij(8)
 
           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)
-          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
-             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))*pm2_5ptr%pm25wc(iaero) + val
-
-          nullify(spm2_5)
-          end do !iaero
-
+          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
@@ -396,32 +390,36 @@ subroutine intpm2_5_(pm2_5head,rval,sval)
              endif
 
 !           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
+             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
+    end if ! laeroana_fv3smoke
 
     if (wrf_mass_regional .and. laeroana_gocart) then
 
diff --git a/src/gsi/setuppm2_5.f90 b/src/gsi/setuppm2_5.f90
index ec390014c..79e6129cd 100644
--- a/src/gsi/setuppm2_5.f90
+++ b/src/gsi/setuppm2_5.f90
@@ -770,9 +770,17 @@ subroutine setuppm2_5(obsLL,odiagLL,lunin,mype,nreal,nobs,isis,is,conv_diagsave)
                 mype,nfldsig)
           call tintrp2a11(pm25wc(:,:,:,2,nfldsig),pm25wc_ges(2),dlat,dlon,dtime,hrdifsig,&
                 mype,nfldsig)
-          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
+
+          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 cd88bfd5f..3483563b9 100644
--- a/src/gsi/stppm2_5.f90
+++ b/src/gsi/stppm2_5.f90
@@ -83,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    
@@ -259,6 +259,21 @@ subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep)
     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
@@ -283,34 +298,25 @@ subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep)
 
                 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 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= 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)* &
-                    (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)
-                end do ! iaero
                 do kk=1,nstep
                    qq=val2+sges(kk)*val
                    pen(kk)=qq*qq*pm2_5ptr%err2
@@ -342,7 +348,10 @@ subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep)
 
        end do
 
-
+       nullify(spm2_5)
+       nullify(rpm2_5)
+       nullify(sdust)
+       nullify(rdust)
 
     end if ! laeroana_fv3smoke
 

From 20ac75e2e86314feedd7e6d6e4688f0eb714d569 Mon Sep 17 00:00:00 2001
From: hongli-wang <hongli.wang@noaa.gov>
Date: Wed, 21 Dec 2022 16:22:05 +0000
Subject: [PATCH 5/8] rebase to the latest EMC GSI develop branch

---
 fix | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fix b/fix
index 023f81d4b..e47bcf18c 160000
--- a/fix
+++ b/fix
@@ -1 +1 @@
-Subproject commit 023f81d4ba631d235156172e04a529a4bf273617
+Subproject commit e47bcf18cc9cab37a1f3ccc198b8e6bab578c9b9

From 34b96be0937cce62ec9a5118023cd577e8a51a3a Mon Sep 17 00:00:00 2001
From: Cory Martin <cory.r.martin@noaa.gov>
Date: Mon, 23 Jan 2023 11:28:31 -0500
Subject: [PATCH 6/8] change submodule for fix to use github

---
 .gitmodules | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.gitmodules b/.gitmodules
index 641f07f69..ea1af1223 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,3 +1,4 @@
 [submodule "fix"]
 	path = fix
-	url = gerrit:GSI-fix
+	url = https://github.com/NOAA-EMC/GSI-fix
+	branch = develop 

From b6ab52aa52b6464d79bb4d1f1c16d75f71a6e67b Mon Sep 17 00:00:00 2001
From: Cory Martin <cory.r.martin@noaa.gov>
Date: Mon, 23 Jan 2023 12:18:10 -0500
Subject: [PATCH 7/8] update submodule to resolve change

---
 fix | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fix b/fix
index e47bcf18c..0be26971f 160000
--- a/fix
+++ b/fix
@@ -1 +1 @@
-Subproject commit e47bcf18cc9cab37a1f3ccc198b8e6bab578c9b9
+Subproject commit 0be26971f834fe9b1d5b118e1e0ffed53facf671

From 08b212f4cc12c326aba3e4bd45cbac282bb61138 Mon Sep 17 00:00:00 2001
From: hongli-wang <hongli.wang@noaa.gov>
Date: Tue, 24 Jan 2023 16:04:00 +0000
Subject: [PATCH 8/8] Change variables(iuv,ndynvario3d,ntracerio3d) to i_kind
 and clear code.

---
 src/gsi/gsi_rfv3io_mod.f90 | 2 +-
 src/gsi/intpm2_5.f90       | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90
index 27efb265c..eb7a86160 100644
--- a/src/gsi/gsi_rfv3io_mod.f90
+++ b/src/gsi/gsi_rfv3io_mod.f90
@@ -891,7 +891,7 @@ 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
+    integer(i_kind):: iuv,ndynvario3d,ntracerio3d
     integer(i_kind):: ntracerchemio3d,ntracersmokeio3d
     integer(i_kind):: loc_id,ncfmt
 
diff --git a/src/gsi/intpm2_5.f90 b/src/gsi/intpm2_5.f90
index d1f3b87c6..42c05ad3f 100644
--- a/src/gsi/intpm2_5.f90
+++ b/src/gsi/intpm2_5.f90
@@ -90,7 +90,6 @@ 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,sdust
     real(r_kind),pointer,dimension(:) :: rpm2_5,rdust