Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Refine PM2.5 DA for the RRFS_SD model #609

Merged
merged 15 commits into from
Sep 29, 2023
Merged

Refine PM2.5 DA for the RRFS_SD model #609

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
c4d84a9
add anowbufr_ext into gsimod namelist
hongli-wang Aug 11, 2023
e66c837
cleanup print sentencce
hongli-wang Aug 11, 2023
f02b817
Cleanup and refine code for PR.
hongli-wang Sep 8, 2023
181b8b8
add anowbufr_ext into gsimod namelist
hongli-wang Aug 11, 2023
f274e81
cleanup print sentencce
hongli-wang Aug 11, 2023
addbb6e
Cleanup and refine code for PR.
hongli-wang Sep 8, 2023
21a61dc
Rebased to develop
hongli-wang Sep 8, 2023
10011ee
add anowbufr_ext into gsimod namelist
hongli-wang Aug 11, 2023
2346f2a
cleanup print sentencce
hongli-wang Aug 11, 2023
9714cc5
Cleanup and refine code for PR.
hongli-wang Sep 8, 2023
126ba88
add anowbufr_ext into gsimod namelist
hongli-wang Aug 11, 2023
3d4a21c
Cleanup and refine code for PR.
hongli-wang Sep 8, 2023
80fb69f
Merge branch 'GSI_SD' of https://github.com/hongli-wang/GSI into GSI_SD
hongli-wang Sep 19, 2023
a2a9de8
Merge remote-tracking branch 'upstream/develop' into GSI_SD
hongli-wang Sep 25, 2023
f4bb476
Remove pm10_da_scheme from namelist
hongli-wang Sep 25, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 29 additions & 11 deletions src/gsi/chemmod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand All @@ -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

Expand All @@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I searched this PR, it looks like we don't need separate berror_fv3_sd_regional and berror_fv3_cmaq_regional?

berror_fv3_sd_regional is only used at Line 460 of src/gsi/m_berror_stats_reg.f90 and there is no difference if we use berror_fv3_cmaq_regional instead. Is it possible to just use berror_fv3_cmaq_regional, or combine them to one logical variable, such as berror_fv3_aero_regional?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The current ones make it more clear to know which B file will be used, and the previous retros runs can repeat without any changes of GSI namelist parameter. This is good to users.
Your suggestion will make the changes be general and good in code development.
So what do you think?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you plan to use different BEC for cmaq and sd? Or will they be the same as they are right now?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@guoqing-noaa
Different B files are used. The one for CMAQ has 70 aerosol species whereas
rrfs-ad only has 3 tracers. I think a separate parameter is good for users. I have no problem in combining then if that is the way you prefer. Thanks.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Got it. Thanks for the explanation!

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,&
Expand All @@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This comment is not related to the code. But just curious whey bumping the max value from 150 to 3000?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This comment is not related to the code. But just curious whey bumping the max value from 150 to 3000?

@guoqing-noaa The default value may not consider the heavy fire event. The value of 3000 is a reasonable PM10 value when a heavy fire event or dust storm takes place.



real(r_kind),parameter :: elev_missing=-9999.0_r_kind
Expand Down Expand Up @@ -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']
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"coarsepm" is added here. However it is never used. And we have to modify Line 282 of src/gsi/setuppm2_5.f90
to remove the "coarsepm" contribution

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, you are correct. This was done by -1 in the do loop. For example,
do i=2,naero_smoke_fv3-1 ! remove contribution from coarsepm

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right, if we will remove coarsepm in src/gsi/setuppm2_5.f90 (line 281), why would we add it here?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@guoqing-noaa
This is because the dust PM10/AOD DA is not included in this PR. But we do consider EnVar to update the coarsepm, possibly in a next PR with FED DA.


! FV3CMAQ
integer(i_kind), parameter :: naero_cmaq_fv3=70 ! !number of cmaq aerosol species aero6
Expand Down Expand Up @@ -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'
Expand All @@ -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
Expand Down
15 changes: 13 additions & 2 deletions src/gsi/gsimod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -181,12 +181,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
Expand Down Expand Up @@ -1585,6 +1588,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
Expand All @@ -1606,13 +1615,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

Expand Down
6 changes: 3 additions & 3 deletions src/gsi/m_berror_stats_reg.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
50 changes: 42 additions & 8 deletions src/gsi/read_anowbufr.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,&
Expand All @@ -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
Expand Down Expand Up @@ -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
hongli-wang marked this conversation as resolved.
Show resolved Hide resolved
end if
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggest reverse the if and else block, i.e: if anowbufr_ext then; else ... (although they function the same, the latter reads more straightforward :) )

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Accepted. Thanks.

else
cycle
endif
Expand All @@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same comments as the above

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Accepted. Thanks.

else
cycle
endif
Expand All @@ -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
Expand All @@ -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)

Expand Down
7 changes: 6 additions & 1 deletion src/gsi/satthin.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
Loading