Skip to content

Commit ea621e1

Browse files
syhasyha
authored
Bug fixes for pm10 assimilation with chem_opt=108 (#1688)
Fixed bugs for the case where PM10 was treated as residuals (=PM10-PM2.5) when assimilated with PM2.5. TYPE: Bug fix KEYWORDS: WRFDA, PM10, chemicda_opt, WRF-Chem assimilation SOURCE: Soyoung Ha (MMM/NCAR) DESCRIPTION OF CHANGES: When assimilated with pm2.5 (e.g., chemicda_opt = 3 or 5) using chem_opt=108, pm10 is treated as (pm10 - pm2.5) residuals in both observations and priors. PM10 obs errors are also adjusted for such residuals. LIST OF MODIFIED FILES: var/da/da_chem_sfc/da_get_innov_vector_chem_sfc.inc var/da/da_obs_io/da_read_obs_chem_sfc.inc TEST CONDUCTED: WRF-Chem/WRFDA was built with gfortran (e.g., configure option 34) on cheyenne. The updated codes passed the cycling tests using chem_opt = 108 and chemicda_opt = 5, resolving the issue of overestimation in surface PM10 concentration.
1 parent c49e863 commit ea621e1

File tree

2 files changed

+99
-62
lines changed

2 files changed

+99
-62
lines changed

‎var/da/da_chem_sfc/da_get_innov_vector_chem_sfc.inc

+49-24
Original file line numberDiff line numberDiff line change
@@ -17,11 +17,13 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv)
1717
real, allocatable :: model_chemic(:,:)
1818
real, allocatable :: model_rho(:,:)
1919

20-
real, parameter :: scaleL = 3000. ! scaling factor for obs error [meters]
20+
! Raidus of influence of ground observations at urban sites (Table 3. in Elbernet al. 2007)
21+
! FIXME: It would be better to put this parameter in the namelist.
22+
real, parameter :: Lrepr = 3000. !2000. ! Scaling factor for obs error [meters]
2123

22-
! max thresholds of abs(o-b) in pm25, pm10, so2, no2, o3, co
24+
! Max thresholds of abs(o-b) in pm25, pm10, so2, no2, o3, co
2325
real,dimension(6):: maxomb = (/100., 100., 0.20, 0.20, 0.20, 20./)
24-
real :: err_o, err_r
26+
real :: err_o1, err_o2, err_r1, err_r2
2527

2628
if (trace_use) call da_trace_entry("da_get_innov_vector_chem_sfc")
2729

@@ -46,14 +48,14 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv)
4648
end do
4749

4850
! [1.1] Compute prior observations at obs sites
49-
if (chem_cv_options == 10) then
51+
if (chem_cv_options == 10) then ! gocart_test
5052
model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_p25)*(model_chemic(:,p_chem_ic_p25) + &
5153
1.375*(96.06/28.964*1000)*model_chemic(:,p_chem_ic_sulf) + &
5254
model_chemic(:,p_chem_ic_bc1)+model_chemic(:,p_chem_ic_bc2)+1.8*(model_chemic(:,p_chem_ic_oc1)+model_chemic(:,p_chem_ic_oc2)) + &
5355
model_chemic(:,p_chem_ic_dust_1)+0.286*model_chemic(:,p_chem_ic_dust_2)+model_chemic(:,p_chem_ic_seas_1) + &
5456
0.942*model_chemic(:,p_chem_ic_seas_2))
5557

56-
else if (chem_cv_options == 20) then
58+
else if (chem_cv_options == 20) then ! mosaic_test
5759

5860
if (chemicda_opt == 1 ) then
5961
model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + &
@@ -109,7 +111,7 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv)
109111

110112
end if
111113

112-
else if (chem_cv_options == 108) then ! racm_soa_vbs_da
114+
else if (chem_cv_options >= 108) then ! racm_soa_vbs_da or racm_soa_vbs_aqchem_da
113115

114116
if (chemicda_opt == 1 .or. chemicda_opt == 3 .or. chemicda_opt == 5) then ! pm2.5
115117

@@ -119,6 +121,13 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv)
119121
model_rho(:,ichem) * model_chemic(:,ichem)
120122
end do
121123

124+
if( (p_chem_ic_p25cwi .gt. p_chem_ic_p25i) .and. chem_cv_options == 109 ) then ! aqueous chemistry
125+
do ichem = p_chem_ic_so4cwj,p_chem_ic_p25cwi
126+
model_chemic_surf(:,p_chemsi_pm25) = model_chemic_surf(:,p_chemsi_pm25) + &
127+
model_rho(:,ichem) * model_chemic(:,ichem)
128+
enddo
129+
endif
130+
122131
end if
123132

124133
if (chemicda_opt == 2) then ! pm10
@@ -129,6 +138,13 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv)
129138
model_rho(:,ichem) * model_chemic(:,ichem)
130139
end do
131140

141+
if( (p_chem_ic_p25cwi .gt. p_chem_ic_p25i) .and. chem_cv_options == 109 ) then ! aqueous chemistry
142+
do ichem = p_chem_ic_so4cwj,p_chem_ic_soilcw
143+
model_chemic_surf(:,p_chemsi_pm10) = model_chemic_surf(:,p_chemsi_pm10) + &
144+
model_rho(:,ichem) * model_chemic(:,ichem)
145+
end do
146+
endif
147+
132148
end if !(chemicda_opt == 2) then
133149

134150
if (chemicda_opt == 3 .or. chemicda_opt == 5) then ! pm10 after pm2.5
@@ -140,6 +156,13 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv)
140156
model_rho(:,ichem) * model_chemic(:,ichem)
141157
end do
142158

159+
if( (p_chem_ic_p25cwi .gt. p_chem_ic_p25i) .and. chem_cv_options == 109 ) then ! aqueous chemistry
160+
do ichem = p_chem_ic_anthcw,p_chem_ic_soilcw
161+
model_chemic_surf(:,p_chemsi_pm10) = model_chemic_surf(:,p_chemsi_pm10) + &
162+
model_rho(:,ichem) * model_chemic(:,ichem)
163+
end do
164+
endif
165+
143166
end if !(chemicda_opt == 3 .or. chemicda_opt == 5)
144167

145168
if (chemicda_opt >= 4) then
@@ -151,9 +174,11 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv)
151174

152175
end if !(chemicda_opt >= 4 ) then
153176

154-
end if ! (chem_cv_options == 108)
177+
end if ! (chem_cv_options >= 108)
178+
179+
180+
! Innovations and obs errors for chemic_surf (e.g., surfce concentration)
155181

156-
! Loop over chemic_surf (e.g., surfce concentration)
157182
do ichem = PARAM_FIRST_SCALAR, num_chemic_surf
158183
ipm = ichem - PARAM_FIRST_SCALAR + 1
159184

@@ -170,22 +195,22 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv)
170195
iv % chemic_surf(n) % chem(ichem) % inv = &
171196
ob % chemic_surf(n) % chem(ichem) - model_chemic_surf(n,ichem)
172197

173-
if ( ichem <= PARAM_FIRST_SCALAR+1 ) then ! PM2.5 and PM10
174-
175-
err_o = 1.5 + 0.0075 * ob % chemic_surf(n) % chem(ichem)
176-
err_r = 0.5 * err_o * sqrt(grid % dx / scaleL)
177-
iv % chemic_surf(n) % chem(ichem) % error = sqrt(err_o**2 + err_r**2)
178-
179-
else ! gas species (so2, no2, o3, co)
180-
181-
err_o = 0.10 * ob % chemic_surf(n) % chem(ichem) ! obs error as 10% of the observation value
182-
iv % chemic_surf(n) % chem(ichem) % error = err_o
183-
184-
!err_o = 1.0 + 0.0075 * ob % chemic_surf(n) % chem(ichem) ! Wei Sun
185-
!err_r = 0.5 * err_o * sqrt(grid % dx / 2000.)
186-
!iv % chemic_surf(n) % chem(ichem) % error = 1.5*sqrt(err_o**2 + err_r**2)
187-
188-
end if !( ichem <= PARAM_FIRST_SCALAR+1 )
198+
! Observation error specification
199+
err_o1 = 1.0 + 0.0075 * ob % chemic_surf(n) % chem(ichem)
200+
err_o2 = 1.5 + 0.0075 * ob % chemic_surf(n) % chem(ichem)
201+
err_r1 = 0.5 * err_o1 * sqrt(grid%dx/Lrepr)
202+
err_r2 = 0.5 * err_o2 * sqrt(grid%dx/Lrepr)
203+
204+
if (chemicda_opt == 3 .or. chemicda_opt == 5 ) then ! PM10 is treated as (PM10 - PM2.5) residual.
205+
if (ichem == PARAM_FIRST_SCALAR+1 ) then
206+
iv % chemic_surf(n) % chem(ichem) % error = 1.0*sqrt(err_o1**2 + err_r1**2) ! reduced error for PM10
207+
else ! PM2.5 and gas species (so2, no2, o3, co)
208+
iv % chemic_surf(n) % chem(ichem) % error = 1.0*sqrt(err_o2**2 + err_r2**2) ! Ha (GMD2022)
209+
end if
210+
else ! PM10 is not a residual.
211+
iv % chemic_surf(n) % chem(ichem) % error = 1.0*sqrt(err_o2**2 + err_r2**2) ! Ha (GMD2022)
212+
!iv % chemic_surf(n) % chem(ichem) % error = 1.5*sqrt(err_o1**2 + err_r1**2) ! Wei's error
213+
end if
189214
190215
! Quality Control
191216
! ----------------

‎var/da/da_obs_io/da_read_obs_chem_sfc.inc

+50-38
Original file line numberDiff line numberDiff line change
@@ -235,65 +235,77 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid)
235235

236236
else if (chem_cv_options == 108) then
237237

238-
read (iunit, fmt='(A6,F10.6,1X,F10.6,1X,I4,3I2,2F9.1,4F9.3)', iostat=iost) &
239-
SiteID, &
240-
platform%info%lat, &
241-
platform%info%lon, &
242-
iyear, imonth, iday, ihour, &
243-
pm25, pm10, so2, no2, o3, co
244-
245-
platform%info%id = trim(SiteID)
246-
read(SiteID,'(I6)') isite
247-
write (platform%info%date_char,'(i4,3(i2.2))') iyear, imonth, iday, ihour
248-
249-
!!!!!!!!!!!!!!!!!!!
250-
! Caution: Observations are assumed to have the same unit of ppmv.
251-
! Unit conversion (to ppmv) must have been done in the preprocessing.
252-
! (ex. The original chinese data had ug/m3 for gas species.)
253-
!!!!!!!!!!!!!!!!!!!
254-
255-
select case ( chemicda_opt )
256-
case ( 1 )
238+
read (iunit, fmt='(A6,F10.6,1X,F10.6,1X,I4,3I2,2F9.1,4F9.3)', iostat=iost) &
239+
SiteID, &
240+
platform%info%lat, &
241+
platform%info%lon, &
242+
iyear, imonth, iday, ihour, &
243+
pm25, pm10, so2, no2, o3, co
244+
245+
platform%info%id = trim(SiteID)
246+
read(SiteID,'(I6)') isite
247+
write (platform%info%date_char,'(i4,3(i2.2))') iyear, imonth, iday, ihour
248+
249+
!!!!!!!!!!!!!!!!!!!
250+
! Caution: Observations are assumed to have the same unit of ppmv.
251+
! Unit conversion (to ppmv) must have been done in the preprocessing.
252+
! (ex. The original chinese data had ug/m3 for gas species.)
253+
!!!!!!!!!!!!!!!!!!!
254+
255+
select case ( chemicda_opt )
256+
case ( 1 )
257257
platform%chem(PARAM_FIRST_SCALAR )%inv = pm25
258258
platform%info%name = "pm25"
259-
case ( 2 )
259+
case ( 2 )
260260
platform%chem(PARAM_FIRST_SCALAR )%inv = pm10
261261
platform%info%name = "pm10"
262-
case ( 3 )
262+
case ( 3 )
263263
platform%chem(PARAM_FIRST_SCALAR )%inv = pm25
264264
platform%chem(PARAM_FIRST_SCALAR+1)%inv = pm10
265265
platform%info%name = "all_pm"
266-
case ( 4 )
266+
case ( 4 )
267267
platform%chem(PARAM_FIRST_SCALAR )%inv = so2
268268
platform%chem(PARAM_FIRST_SCALAR+1)%inv = no2
269269
platform%chem(PARAM_FIRST_SCALAR+2)%inv = o3
270270
platform%chem(PARAM_FIRST_SCALAR+3)%inv = co
271271
platform%info%name = "gas"
272-
case ( 5 )
272+
case ( 5 )
273273
platform%chem(PARAM_FIRST_SCALAR )%inv = pm25
274274
platform%chem(PARAM_FIRST_SCALAR+1)%inv = pm10
275275
platform%chem(PARAM_FIRST_SCALAR+2)%inv = so2
276276
platform%chem(PARAM_FIRST_SCALAR+3)%inv = no2
277277
platform%chem(PARAM_FIRST_SCALAR+4)%inv = o3
278278
platform%chem(PARAM_FIRST_SCALAR+5)%inv = co
279279
platform%info%name = "all"
280-
case default
280+
case default
281281
platform%chem(PARAM_FIRST_SCALAR )%inv = pm25
282282
platform%info%name = ""
283-
end select
284-
285-
! Sanity check
286-
!!!!!!!!!!!!!!!!!!!
287-
do ichem = PARAM_FIRST_SCALAR, num_chemic_surf
288-
ipm = ichem - PARAM_FIRST_SCALAR + 1
289-
if (platform%chem(ichem)%inv<0.or.abs(platform%chem(ichem)%inv)>=max_pm(ipm)) then
290-
platform%chem(ichem)%inv=missing_r
291-
else
292-
platform%chem(ichem)%qc = 0
293-
endif
294-
end do
295-
296-
end if ! (chem_cv_options == 10) then
283+
end select
284+
285+
! Sanity check
286+
!!!!!!!!!!!!!!!!!!!
287+
do ichem = PARAM_FIRST_SCALAR, num_chemic_surf
288+
ipm = ichem - PARAM_FIRST_SCALAR + 1
289+
if (platform%chem(ichem)%inv<0.or.abs(platform%chem(ichem)%inv)>=max_pm(ipm)) then
290+
platform%chem(ichem)%inv=missing_r
291+
else
292+
platform%chem(ichem)%qc = 0
293+
endif
294+
end do
295+
296+
if ((chemicda_opt == 3) .or. (chemicda_opt == 5)) then
297+
if ((platform%chem(PARAM_FIRST_SCALAR+1)%inv.lt.platform%chem(PARAM_FIRST_SCALAR)%inv) .or. &
298+
(platform%chem(PARAM_FIRST_SCALAR+1)%inv<0).or.(platform%chem(PARAM_FIRST_SCALAR)%inv<0)) then ! Feb-27-2022
299+
platform%chem(PARAM_FIRST_SCALAR+1)%inv = missing_r
300+
platform%chem(PARAM_FIRST_SCALAR+1)%qc = missing
301+
else
302+
! pm10 <= pm10 - pm25 residual
303+
platform%chem(PARAM_FIRST_SCALAR+1)%inv = platform%chem(PARAM_FIRST_SCALAR+1)%inv - &
304+
platform%chem(PARAM_FIRST_SCALAR)%inv
305+
endif
306+
end if ! ((chemicda_opt == 3) .or. (chemicda_opt == 5)) then
307+
308+
end if ! if (chem_cv_options == 108)
297309

298310
if (iost /= 0) then
299311
! FIX? This is expected, but its unclear how we handle failure

0 commit comments

Comments
 (0)