Refine PM2.5 DA for the RRFS_SD model (#609)

**Description**

Refine the PM2.5 DA for the RRFS_SD model by use of veg_type, which is
used to decide whether the obs is in urban area or not. Different
thresholds for innovations outside/inside urban areas will be used. Add
new namelist parameters, such as threshold for innovations, anowbufr
type
Read in station terrain height, PM10 et al if extended BUFR format for
anow air quality data is used

Fixes #606 


**Type of change**

Please delete options that are not relevant.

- [ ] Bug fix (non-breaking change which fixes an issue)
- [X ] New feature (non-breaking change which adds functionality)
- [ ] Breaking change (fix or feature that would cause existing
functionality to not work as expected)
- [ ] This change requires a documentation update

**How Has This Been Tested?**

<!-- This change will only affect PM2.5 DA for RRFS_SD, so it will not
affect any other tests -->

  
**Checklist**

- [X ] My code follows the style guidelines of this project
- [X] I have performed a self-review of my own code
- [X] I have commented my code, particularly in hard-to-understand areas
- [ ] New and existing tests pass with my changes
- [ ] Any dependent changes have been merged and published

**DUE DATE for this PR is 9/22/2023.** If this PR is not merged into
develop by this date, the PR will be closed and returned to the
developer.

src/gsi/chemmod.f90

@@ -40,21 +40,23 @@ module chemmod
   public :: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3
 
 ! fv3smoke
-  public :: naero_smoke_fv3,aeronames_smoke_fv3,pm2_5_innov_threshold 
+  public :: naero_smoke_fv3,aeronames_smoke_fv3
+  public :: pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_threshold 
+  public :: pm10_innov_threshold,pm10_urban_innov_threshold,pm10_bg_threshold,pm10_obs_threshold
 
   public :: naero_gocart_wrf,aeronames_gocart_wrf
 
   public :: pm2_5_guess,init_pm2_5_guess,&
        aerotot_guess,init_aerotot_guess
   public :: init_chem
-  public :: berror_chem,berror_fv3_cmaq_regional,oneobtest_chem,maginnov_chem,magoberr_chem,oneob_type_chem,conconeobs
+  public :: berror_chem,berror_fv3_cmaq_regional,berror_fv3_sd_regional,oneobtest_chem,maginnov_chem,magoberr_chem,oneob_type_chem,conconeobs
   public :: oblat_chem,oblon_chem,obpres_chem,diag_incr,oneobschem
   public :: site_scale,nsites
   public :: tunable_error
   public :: in_fname,out_fname,incr_fname,maxstr
   public :: code_pm25_ncbufr,code_pm25_anowbufr
   public :: code_pm10_ncbufr,code_pm10_anowbufr
-
+  public :: anowbufr_ext
   public :: l_aoderr_table
 
   public :: laeroana_gocart,laeroana_fv3cmaq,laeroana_fv3smoke,crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, &
@@ -79,7 +81,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 
+  real(r_kind)    :: pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_threshold
+  real(r_kind)    :: pm10_innov_threshold,pm10_urban_innov_threshold,pm10_bg_threshold,pm10_obs_threshold
 
   logical :: wrf_pm2_5
 
@@ -90,7 +93,9 @@ module chemmod
 
   logical :: aero_ratios
 
-  logical :: oneobtest_chem,diag_incr,berror_chem,berror_fv3_cmaq_regional
+  logical :: oneobtest_chem,diag_incr,berror_chem
+  logical :: berror_fv3_cmaq_regional,berror_fv3_sd_regional
+  logical :: anowbufr_ext
   character(len=max_varname_length) :: oneob_type_chem
   integer(i_kind), parameter :: maxstr=256
   real(r_kind) :: maginnov_chem,magoberr_chem,conconeobs,&
@@ -103,7 +108,7 @@ module chemmod
   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
+  real(r_kind),parameter :: pm10_teom_max=3000.0_r_kind !ug/m3
 
 
   real(r_kind),parameter :: elev_missing=-9999.0_r_kind
@@ -157,10 +162,10 @@ module chemmod
         'AOLGAJ', 'AISO1J', 'AISO2J', 'AISO3J',   'ATRP1J', 'ATRP2J',&
         'ASQTJ',  'AOLGBJ', 'AORGCJ']
 ! fv3smoke
-  integer(i_kind), parameter ::  naero_smoke_fv3=2 
+  integer(i_kind), parameter ::  naero_smoke_fv3=3 
 
   character(len=max_varname_length), dimension(naero_smoke_fv3), parameter :: &
-        aeronames_smoke_fv3=[character(len=max_varname_length) ::  'smoke','dust' ]
+        aeronames_smoke_fv3=[character(len=max_varname_length) ::  'smoke','dust','coarsepm']
 
 ! FV3CMAQ
   integer(i_kind), parameter :: naero_cmaq_fv3=70 !  !number of cmaq aerosol species aero6
@@ -286,8 +291,15 @@ subroutine init_chem
 !initialiazes default values to &CHEM namelist parameters
 
     berror_chem=.false.
-    berror_fv3_cmaq_regional=.false.  ! Set .true. to use berror for fv3_cmaq_regional, whose cv has 10 characters
+    berror_fv3_cmaq_regional=.false.  ! .False. : Dont perform aerosal DA for the online RRFS_CMAQ model so dont need to read in B for RRFS_CMAQ.  
+                                      ! .true.  : Use berror for fv3_cmaq_regional, whose cv has 10 characters
+    berror_fv3_sd_regional=.false.    ! .False. : Dont perform aerosal DA for the RRFS_SD model so dont need to read in B for RRFS_SD.  
+                                      ! .true. to use berror for rrfs_sd model, whose cv has 10 characters
     oneobtest_chem=.false.
+    anowbufr_ext=.false.              ! .False. : use default anowbufr data
+                                      ! .True.  : use the extented bufr data
+                                      ! that includes PM10, station elevation
+                                      ! etal in addition to pm2.5. 
     maginnov_chem=30_r_kind
     magoberr_chem=2_r_kind
     oneob_type_chem='pm2_5'
@@ -307,9 +319,15 @@ 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
+    pm2_5_innov_threshold = 15.0_r_kind
+    pm2_5_urban_innov_threshold = 30.0_r_kind
+    pm2_5_bg_threshold = 2.0_r_kind
+    pm10_innov_threshold = 15.0_r_kind
+    pm10_urban_innov_threshold = 30.0_r_kind
+    pm10_bg_threshold = 2.0_r_kind
+    pm10_obs_threshold = 140.0_r_kind ! Barry's manuscript
     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
+    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
     raod_radius_std_scale  = 1.0_r_kind  ! Tune standard deviation of particles when calculating AOD using CRTM with CMAQ LUTs.
     aod_qa_limit = 3

src/gsi/gsimod.F90

@@ -182,12 +182,15 @@ module gsimod
   use gsi_chemguess_mod, only: gsi_chemguess_init,gsi_chemguess_final
   use tcv_mod, only: init_tcps_errvals,tcp_refps,tcp_width,tcp_ermin,tcp_ermax
   use chemmod, only : init_chem,berror_chem,berror_fv3_cmaq_regional,oneobtest_chem,&
+       berror_fv3_sd_regional,&
        maginnov_chem,magoberr_chem,&
        oneob_type_chem,oblat_chem,&
+       anowbufr_ext,&
        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,pm2_5_innov_threshold,crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, &
+       laeroana_fv3cmaq,laeroana_fv3smoke,pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_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
@@ -1594,6 +1597,12 @@ module gsimod
 ! chem(options for gsi chem analysis) :
 !     berror_chem       - .true. when background  for chemical species that require
 !                          conversion to lower case and/or species names longer than 5 chars
+!     berror_fv3_cmaq_regional   - .true. use background error stat for online
+!                                         RRFS_CMAQ model. Control variable
+!                                         names extended up to 10 chars
+!     berror_fv3_sd_regional     - .true. use background error stat for online
+!                                         RRFS_SD model. Control variable
+!                                         names extended up to 10 chars
 !     oneobtest_chem    - one-ob trigger for chem constituent analysis
 !     maginnov_chem     - O-B make-believe residual for one-ob chem test
 !     magoberr_chem     - make-believe obs error for one-ob chem test
@@ -1615,13 +1624,15 @@ module gsimod
 !     luse_deepblue     - whether to use MODIS AOD from the deepblue   algorithm
 !     lread_ext_aerosol - if true, reads aerfNN file for aerosol arrays rather than sigfNN (NGAC NEMS IO)
 
-  namelist/chem/berror_chem,berror_fv3_cmaq_regional,oneobtest_chem,maginnov_chem,magoberr_chem,&
+  namelist/chem/berror_chem,berror_fv3_cmaq_regional,berror_fv3_sd_regional,& 
+       oneobtest_chem,anowbufr_ext,maginnov_chem,magoberr_chem,&
        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,laeroana_fv3smoke,l_aoderr_table, aod_qa_limit, &
        crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file, &
        icvt_cmaq_fv3,pm2_5_innov_threshold, &
+       pm2_5_innov_threshold,pm2_5_urban_innov_threshold,pm2_5_bg_threshold,&
        raod_radius_mean_scale,raod_radius_std_scale, luse_deepblue,&
        aero_ratios,wrf_pm2_5, lread_ext_aerosol
 

src/gsi/m_berror_stats_reg.f90

@@ -12,7 +12,7 @@ module m_berror_stats_reg
       use kinds,only : i_kind,r_kind
       use constants, only: zero,one,max_varname_length,half
       use gridmod, only: nsig
-      use chemmod, only : berror_chem,berror_fv3_cmaq_regional,upper2lower,lower2upper
+      use chemmod, only : berror_chem,berror_fv3_cmaq_regional,berror_fv3_sd_regional,upper2lower,lower2upper
       use m_berror_stats, only: usenewgfsberror,berror_stats
 
       implicit none
@@ -312,7 +312,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
       use constants, only: zero,one,ten,three
       use mpeu_util,only: getindex
       use radiance_mod, only: icloud_cv,n_clouds_fwd,cloud_names_fwd
-      use chemmod, only: berror_fv3_cmaq_regional
+      use chemmod, only: berror_fv3_cmaq_regional,berror_fv3_sd_regional
 
       implicit none
 
@@ -466,7 +466,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
         var=upper2lower(varshort)
         if (trim(var) == 'pm25') var = 'pm2_5'
      else 
-        if ( berror_fv3_cmaq_regional) then
+        if ( berror_fv3_cmaq_regional .or. berror_fv3_sd_regional) then
           read(inerr,iostat=istat) varlong, isig
           var=varlong
         else

src/gsi/read_anowbufr.f90

@@ -50,7 +50,9 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,&
         iconc,ierror,ilat,ilon,itime,iid,ielev,isite,iikx,ilate,ilone,&
         elev_missing,site_scale,tunable_error,&
         code_pm25_ncbufr,code_pm25_anowbufr,&
-        code_pm10_ncbufr,code_pm10_anowbufr
+        code_pm10_ncbufr,code_pm10_anowbufr,&
+        anowbufr_ext,pm2_5_teom_max,pm10_teom_max
+
   use mpimod, only: npe
 
   implicit none
@@ -71,10 +73,12 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,&
         nyob=3,ndhr=4,ntyp=5,ncopopm=6
 !see headr input format below
   integer(i_kind), parameter :: nfields=6
+  integer(i_kind), parameter :: nfields_b=12
 !output format parameters
   integer(i_kind), parameter:: nchanl=0,nreal=ilone
 
   real(r_kind),parameter :: r360 = 360.0_r_kind
+  real(r_kind),parameter :: r90 = 90.0_r_kind
   real(r_kind),parameter :: percent=1.e-2_r_kind
 
   real(r_kind), parameter :: anow_missing=1.0e11_r_kind,&
@@ -96,8 +100,10 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,&
   real(r_kind), dimension(5)  :: rinc
   character(len=8) :: subset
   character(len=80) :: headr
+  character(len=80) :: obstr
 
   real(r_double), dimension(nfields) :: indata
+  real(r_double), dimension(nfields_b) :: indata_a,indata_b
 
   real(r_kind) :: tdiff,obstime,t4dv
   real(r_kind) :: dlat,dlon,error_1,error_2,obserror,dlat_earth,dlon_earth
@@ -141,17 +147,23 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,&
 ! reading each report from bufr
 
   do while (ireadmg(lunin,subset,idate) == 0)
-
      if (trim(obstype)=='pm2_5') then
 
         if ( (subset == 'NC008031') .or. (subset == 'NC008032' ) ) then
            headr='PTID CLONH CLATH TPHR TYPO COPOPM'
            ncbufr=.true.
            write(6,*)'READ_PM2_5:  AIRNOW data type, subset=',subset
         else if (subset == 'ANOWPM') then
-           headr='SID XOB YOB DHR TYP COPOPM'
-           anowbufr=.true.
-           write(6,*)'READ_PM2_5:  AIRNOW data type, subset=',subset
+           if (anowbufr_ext) then
+              headr='SID XOB YOB DHR TYP T29 SQN PROCN RPT CAT TYPO TSIG'
+              obstr='TPHR QCIND COPOPM ELV COPOPM10 COPOCO'
+              anowbufr=.true.
+              write(6,*)'READ_PM2_5_BUFR_EXT:  AIRNOW data type, subset=',subset 
+           else ! default ANOWBUFR Table
+              headr='SID XOB YOB DHR TYP COPOPM'
+              anowbufr=.true.
+              write(6,*)'READ_PM2_5:  AIRNOW data type, subset=',subset
+           end if
         else
            cycle
         endif
@@ -162,6 +174,17 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,&
            headr='PTID CLONH CLATH TPHR TYPO COPOPM'
            ncbufr=.true.
            write(6,*)'READ_PM10:  AIRNOW data type, subset=',subset
+        else if (subset == 'ANOWPM') then
+           if (anowbufr_ext) then
+              headr='SID XOB YOB DHR TYP T29 SQN PROCN RPT CAT TYPO TSIG'
+              obstr='TPHR QCIND COPOPM ELV COPOPM10 COPOCO'
+              anowbufr=.true.
+              write(6,*)'READ_PM10_BUFR_EXT:  AIRNOW data type, subset=',subset
+           else
+              headr='SID XOB YOB DHR TYP COPOPM'
+              anowbufr=.true.
+              write(6,*)'READ_PM10:  AIRNOW data type, subset=',subset
+           end if
         else
            cycle
         endif
@@ -176,8 +199,17 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,&
      imin=0
 
      do while (ireadsb(lunin) == 0)
-        call ufbint(lunin,indata,nfields,1,iret,headr)
-
+        if (anowbufr_ext) then
+           call ufbint(lunin,indata_a,nfields_b,1,iret,headr)
+           indata(1:5) = indata_a(1:5)
+           call ufbint(lunin,indata_b,nfields_b,1,iret,obstr)
+           if (trim(obstype)=='pm2_5') indata(ncopopm)=indata_b(3)
+           if (trim(obstype)=='pm10')  indata(ncopopm)=indata_b(5)
+           site_elev = indata_b(4)
+        else
+           call ufbint(lunin,indata,nfields,1,iret,headr)
+        end if
+
         if (anowbufr) then
            kx=indata(ntyp)
            read(sid,'(Z8)')site_id
@@ -198,13 +230,15 @@ subroutine read_anowbufr(nread,ndata,nodata,gstime,&
 
         nread = nread + 1
         conc=indata(ncopopm)
-
         if ( iret > 0 .and. (conc < conc_missing ) .and. &
                (conc >= zero)) then
 
            if(indata(nxob) >= r360)  indata(nxob) = indata(nxob) - r360
            if(indata(nxob) <  zero)  indata(nxob) = indata(nxob) + r360
 
+           if(indata(nxob) > r360)cycle
+           if(indata(nyob) > r90)cycle
+
            dlon_earth_deg=indata(nxob)
            dlat_earth_deg=indata(nyob)
 

src/gsi/satthin.F90

@@ -134,6 +134,8 @@ module satthin
   use obsmod, only: time_window_max
   use constants, only: deg2rad,rearth_equator,zero,two,pi,half,one,&
        rad2deg,r1000
+  use chemmod, only: laeroana_fv3smoke
+
   implicit none
 
 ! set default to private
@@ -961,7 +963,10 @@ subroutine getsfc(mype,mype_io,use_sfc,use_sfc_any)
     end if
     if (.not.lobserver) then
        if(allocated(veg_frac)) deallocate(veg_frac)
-       if(allocated(veg_type)) deallocate(veg_type)
+!       veg_type will be used in setuppm2_5.f90 for rrfs_sd PM2.5 DA
+       if(.not. laeroana_fv3smoke )then
+         if(allocated(veg_type)) deallocate(veg_type)
+       endif
        if(allocated(soil_type)) deallocate(soil_type)
        if(allocated(soil_moi)) deallocate(soil_moi)
        if(allocated(sfc_rough)) deallocate(sfc_rough)