module_mosaic_addemiss.F
References to this file elsewhere.
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! MOSAIC module: see module_mosaic_driver.F for information and terms of use
8 !**********************************************************************************
9 MODULE module_mosaic_addemiss
10 !WRF:MODEL_LAYER:CHEMICS
11
12 ! rce 2005-feb-18 - one fix for indices of volumcen_sect, [now (isize,itype)]
13 ! rce 2005-jan-14 - added subr mosaic_seasalt_emiss (and a call to it)
14 ! rce 2004-dec-03 - many changes associated with the new aerosol "pointer"
15 ! variables in module_data_mosaic_asect
16
17
18
19 integer, parameter :: mosaic_addemiss_active = 1
20 ! only do emissions when this is positive
21 ! (when it is negative, emissions tendencies are zero)
22
23 integer, parameter :: mosaic_addemiss_masscheck = -1
24 ! only do emissions masscheck calcs when this is positive
25
26 integer, parameter :: mosaic_seasalt_emiss_active = 1
27 ! only do seasalt emissions when this is positive
28
29
30 CONTAINS
31
32
33
34 !----------------------------------------------------------------------
35 subroutine mosaic_addemiss( id, dtstep, u10, v10, alt, dz8w, xland, &
36 config_flags, chem, &
37 e_pm10,e_pm25,e_pm25i,e_pm25j,e_eci,e_ecj,e_orgi,e_orgj, &
38 e_so4j,e_so4c,e_no3j,e_no3c,e_orgc,e_ecc, &
39 ids,ide, jds,jde, kds,kde, &
40 ims,ime, jms,jme, kms,kme, &
41 its,ite, jts,jte, kts,kte )
42 !
43 ! adds emissions for mosaic aerosol species
44 ! (i.e., emissions tendencies over time dtstep are applied
45 ! to the aerosol concentrations)
46 !
47
48 USE module_configure, only: grid_config_rec_type
49 USE module_state_description, only: num_chem, param_first_scalar, &
50 emiss_inpt_default, emiss_inpt_pnnl_rs, emiss_inpt_pnnl_cm
51 USE module_data_mosaic_asect
52
53 IMPLICIT NONE
54
55 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
56
57 INTEGER, INTENT(IN ) :: id, &
58 ids,ide, jds,jde, kds,kde, &
59 ims,ime, jms,jme, kms,kme, &
60 its,ite, jts,jte, kts,kte
61
62 REAL, INTENT(IN ) :: dtstep
63
64 ! 10-m wind speed components (m/s)
65 REAL, DIMENSION( ims:ime , jms:jme ) , &
66 INTENT(IN ) :: u10, v10, xland
67
68 ! trace species mixing ratios (aerosol mass = ug/kg-air; number = #/kg-air)
69 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
70 INTENT(INOUT ) :: chem
71 !
72 ! aerosol emissions arrays ((ug/m3)*m/s)
73 !
74 ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
75 REAL, DIMENSION( ims:ime, kms:config_flags%kemit, jms:jme ), &
76 INTENT(IN ) :: &
77 e_pm10,e_pm25,e_pm25i,e_pm25j,e_eci,e_ecj,e_orgi,e_orgj, &
78 e_so4j,e_so4c,e_no3j,e_no3c,e_orgc,e_ecc
79
80 ! 1/(dry air density) and layer thickness (m)
81 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
82 INTENT(IN ) :: alt, dz8w
83
84 ! local variables
85 integer i, j, k, l, n
86 integer iphase, itype
87 integer p1st
88
89 real, parameter :: efact1 = 1.0
90 real :: aem_so4, aem_no3, aem_cl, aem_msa, aem_co3, aem_nh4, &
91 aem_na, aem_ca, aem_oin, aem_oc, aem_bc, aem_num
92 real dum, fact
93
94
95 ! fraction of sorgam i/aitken mode emissions that go to each
96 ! of the mosaic 8 "standard" sections
97 real, save :: fr8b_aem_sorgam_i(8) = &
98 (/ 0.965, 0.035, 0.000, 0.000, &
99 0.000, 0.000, 0.000, 0.000 /)
100
101 ! fraction of sorgam j/accum mode emissions that go to each
102 ! of the mosaic 8 "standard" sections
103 real, save :: fr8b_aem_sorgam_j(8) = &
104 (/ 0.026, 0.147, 0.350, 0.332, &
105 0.125, 0.019, 0.001, 0.000/)
106
107 ! fraction of sorgam coarse mode emissions that go to each
108 ! of the mosaic 8 "standard" sections
109 real, save :: fr8b_aem_sorgam_c(8) = &
110 (/ 0.000, 0.000, 0.000, 0.002, &
111 0.021, 0.110, 0.275, 0.592 /)
112
113 ! fraction of mosaic fine (< 2.5 um) emissions that go to each
114 ! of the mosaic 8 "standard" sections
115 !wig 1-Apr-2005, Updated fractional breakdown between bins. (See also
116 ! bdy_chem_value_mosaic and mosaic_init_wrf_mixrats_opt2
117 ! in module_mosaic_initmixrats.F.) Note that the values
118 ! here no longer match the other two subroutines.
119 !rce 10-may-2005, changed fr8b_aem_mosaic_f & fr8b_aem_mosaic_c
120 ! to values determined by jdf
121 real, save :: fr8b_aem_mosaic_f(8) = &
122 (/ 0.060, 0.045, 0.245, 0.400, 0.100, 0.150, 0., 0./) !10-may-2005
123 ! (/ 0.100, 0.045, 0.230, 0.375, 0.100, 0.150, 0., 0./) !1-Apr-2005 values
124 ! (/ 0.0275, 0.0426, 0.2303, 0.3885, 0.1100, 0.2011, 0., 0./) !15-Nov-2004 values
125 ! (/ 0.01, 0.05, 0.145, 0.60, 0.145, 0.05, 0.00, 0.00 /)
126 ! (/ 0.04, 0.10, 0.35, 0.29, 0.15, 0.07, 0.0, 0.0 /)
127
128 ! fraction of mosaic coarse (> 2.5 um) emissions that go to each
129 ! of the mosaic 8 "standard" sections
130 real, save :: fr8b_aem_mosaic_c(8) = &
131 (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.300, 0.700 /) ! 10-may-2005
132 ! (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.933, 0.067 /) ! as of apr-2005
133 ! (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.16, 0.84 /) ! "old"
134
135 ! following 5 arrays correspond to the above "fr8b_" arrays
136 ! but are set at run time base on input (namelist) parameters:
137 ! only the sorgam or mosaic arrays are non-zero, depending on
138 ! emiss_inpt_opt
139 ! when nsize_aer=4 (=number of sections), the values are
140 ! calculated by adding two of the 8-section values
141 real :: fr_aem_sorgam_i(8)
142 real :: fr_aem_sorgam_j(8)
143 real :: fr_aem_sorgam_c(8)
144 real :: fr_aem_mosaic_f(8)
145 real :: fr_aem_mosaic_c(8)
146 double precision :: chem_sum(num_chem)
147
148 character*80 msg
149
150
151 ! *** currently only works with ntype_aer = 1
152 itype = 1
153 iphase = ai_phase
154
155
156 !
157 ! compute factors used for apportioning either
158 ! the MADE-SORGAM emissions (i=aitken, j=accum, c=coarse modes) OR
159 ! the MOSAIC emission (f=fine (< 2.5 um), c=coarse (> 2.5 um))
160 ! to each size section
161 !
162 ! note: the fr8b_aer_xxxxxx_y values are specific to the mosaic 8 bin
163 ! structure with dlo_sect(1)=0.039 um and dhi_sect(8)=10.0 um,
164 ! also, the fr8b_aem_sorgam_y are specific for the assumed
165 ! dgvem_i/j/c used in the sorgam code
166 ! also, the fr8b_aem_mosaic_y values are specific for the assumed (by us)
167 ! size distribution for fine and coarse primary emissions
168 !
169 ! when there are 4 bins (nsize_aer=4), each of these "wider" bins
170 ! corresponds to 2 of the "narrower" bins of the 8 bin structure
171 !
172 ! note: if fr_aem_sorgam_y > 0, then fr_aem_mosaic_y = 0, and vice-versa
173 !
174 if ((nsize_aer(itype) .ne. 4) .and. (nsize_aer(itype) .ne. 8)) then
175 write(msg,'(a,i5)') &
176 'subr mosaic_addemiss - nsize_aer(itype) must be ' // &
177 '4 or 8 but = ', &
178 nsize_aer(itype)
179 call wrf_error_fatal( msg )
180 end if
181
182 fr_aem_sorgam_i(:) = 0.0
183 fr_aem_sorgam_j(:) = 0.0
184 fr_aem_sorgam_c(:) = 0.0
185 fr_aem_mosaic_f(:) = 0.0
186 fr_aem_mosaic_c(:) = 0.0
187
188 emiss_inpt_select_1: SELECT CASE( config_flags%emiss_inpt_opt )
189
190 CASE( emiss_inpt_default, emiss_inpt_pnnl_rs )
191 if (nsize_aer(itype) .eq. 8) then
192 fr_aem_sorgam_i(:) = fr8b_aem_sorgam_i(:)
193 fr_aem_sorgam_j(:) = fr8b_aem_sorgam_j(:)
194 fr_aem_sorgam_c(:) = fr8b_aem_sorgam_c(:)
195 else if (nsize_aer(itype) .eq. 4) then
196 do n = 1, nsize_aer(itype)
197 fr_aem_sorgam_i(n) = fr8b_aem_sorgam_i(2*n-1) &
198 + fr8b_aem_sorgam_i(2*n)
199 fr_aem_sorgam_j(n) = fr8b_aem_sorgam_j(2*n-1) &
200 + fr8b_aem_sorgam_j(2*n)
201 fr_aem_sorgam_c(n) = fr8b_aem_sorgam_c(2*n-1) &
202 + fr8b_aem_sorgam_c(2*n)
203 end do
204 end if
205
206 CASE( emiss_inpt_pnnl_cm )
207 if (nsize_aer(itype) .eq. 8) then
208 fr_aem_mosaic_f(:) = fr8b_aem_mosaic_f(:)
209 fr_aem_mosaic_c(:) = fr8b_aem_mosaic_c(:)
210 else if (nsize_aer(itype) .eq. 4) then
211 do n = 1, nsize_aer(itype)
212 fr_aem_mosaic_f(n) = fr8b_aem_mosaic_f(2*n-1) &
213 + fr8b_aem_mosaic_f(2*n)
214 fr_aem_mosaic_c(n) = fr8b_aem_mosaic_c(2*n-1) &
215 + fr8b_aem_mosaic_c(2*n)
216 end do
217 end if
218
219 CASE DEFAULT
220 return
221
222 END SELECT emiss_inpt_select_1
223
224 ! when mosaic_addemiss_active <= 0, set fr's to zero,
225 ! which causes the changes to chem(...) to be zero
226 if (mosaic_addemiss_active <= 0) then
227 fr_aem_sorgam_i(:) = 0.0
228 fr_aem_sorgam_j(:) = 0.0
229 fr_aem_sorgam_c(:) = 0.0
230 fr_aem_mosaic_f(:) = 0.0
231 fr_aem_mosaic_c(:) = 0.0
232 end if
233
234
235 ! do mass check initial calc
236 if (mosaic_addemiss_masscheck > 0) call addemiss_masscheck( &
237 id, config_flags, 1, 'mosaic_ademiss', &
238 dtstep, efact1, dz8w, chem, chem_sum, &
239 ids,ide, jds,jde, kds,kde, &
240 ims,ime, jms,jme, kms,kme, &
241 its,ite, jts,jte, kts,kte, &
242 14, &
243 e_pm10,e_pm25,e_pm25i,e_pm25j,e_eci,e_ecj,e_orgi,e_orgj, &
244 e_so4j,e_so4c,e_no3j,e_no3c,e_orgc,e_ecc, &
245 e_ecc,e_ecc,e_ecc,e_ecc,e_ecc,e_ecc,e_ecc )
246
247 p1st = param_first_scalar
248
249 !
250 ! apply emissions at each section and grid point
251 !
252 do 1900 n = 1, nsize_aer(itype)
253
254 do 1830 j = jts, jte
255 do 1820 k = kts, min(config_flags%kemit,kte-1)
256 do 1810 i = its, ite
257
258 ! compute mass emissions [(ug/m3)*m/s] for each species
259 ! using the apportioning fractions
260 aem_so4 = fr_aem_mosaic_f(n)*e_so4j(i,k,j) &
261 + fr_aem_mosaic_c(n)*e_so4c(i,k,j)
262
263 aem_no3 = fr_aem_mosaic_f(n)*e_no3j(i,k,j) &
264 + fr_aem_mosaic_c(n)*e_no3c(i,k,j)
265
266 aem_oc = fr_aem_mosaic_f(n)*e_orgj(i,k,j) &
267 + fr_aem_mosaic_c(n)*e_orgc(i,k,j) &
268 + fr_aem_sorgam_i(n)*e_orgi(i,k,j) &
269 + fr_aem_sorgam_j(n)*e_orgj(i,k,j)
270
271 aem_bc = fr_aem_mosaic_f(n)*e_ecj(i,k,j) &
272 + fr_aem_mosaic_c(n)*e_ecc(i,k,j) &
273 + fr_aem_sorgam_i(n)*e_eci(i,k,j) &
274 + fr_aem_sorgam_j(n)*e_ecj(i,k,j)
275
276 aem_oin = fr_aem_mosaic_f(n)*e_pm25j(i,k,j) &
277 + fr_aem_mosaic_c(n)*e_pm10(i,k,j) &
278 + fr_aem_sorgam_i(n)*e_pm25i(i,k,j) &
279 + fr_aem_sorgam_j(n)*e_pm25j(i,k,j) &
280 + fr_aem_sorgam_c(n)*e_pm10(i,k,j)
281
282 ! emissions for these species are currently zero
283 aem_nh4 = 0.0
284 aem_na = 0.0
285 aem_cl = 0.0
286 aem_ca = 0.0
287 aem_co3 = 0.0
288 aem_msa = 0.0
289
290 ! compute number emissions
291 ! first sum the mass-emissions/density
292 aem_num = &
293 (aem_so4/dens_so4_aer) + (aem_no3/dens_no3_aer) + &
294 (aem_cl /dens_cl_aer ) + (aem_msa/dens_msa_aer) + &
295 (aem_co3/dens_co3_aer) + (aem_nh4/dens_nh4_aer) + &
296 (aem_na /dens_na_aer ) + (aem_ca /dens_ca_aer ) + &
297 (aem_oin/dens_oin_aer) + (aem_oc /dens_oc_aer ) + &
298 (aem_bc /dens_bc_aer )
299
300 ! then multiply by 1.0e-6 to convert ug to g
301 ! and divide by particle volume at center of section (cm3)
302 aem_num = aem_num*1.0e-6/volumcen_sect(n,itype)
303
304 ! apply the emissions and convert from flux to mixing ratio
305 ! fact = (dtstep/dz8w(i,k,j))*(28.966/1000.)
306 fact = (dtstep/dz8w(i,k,j))*alt(i,k,j)
307
308 ! rce 22-nov-2004 - change to using the "..._aer" species pointers
309 l = lptr_so4_aer(n,itype,iphase)
310 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_so4*fact
311
312 l = lptr_no3_aer(n,itype,iphase)
313 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_no3*fact
314
315 l = lptr_cl_aer(n,itype,iphase)
316 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_cl*fact
317
318 l = lptr_msa_aer(n,itype,iphase)
319 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_msa*fact
320
321 l = lptr_co3_aer(n,itype,iphase)
322 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_co3*fact
323
324 l = lptr_nh4_aer(n,itype,iphase)
325 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_nh4*fact
326
327 l = lptr_na_aer(n,itype,iphase)
328 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_na*fact
329
330 l = lptr_ca_aer(n,itype,iphase)
331 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_ca*fact
332
333 l = lptr_oin_aer(n,itype,iphase)
334 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oin*fact
335
336 l = lptr_oc_aer(n,itype,iphase)
337 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oc*fact
338
339 l = lptr_bc_aer(n,itype,iphase)
340 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_bc*fact
341
342 l = numptr_aer(n,itype,iphase)
343 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_num*fact
344
345 1810 continue
346 1820 continue
347 1830 continue
348
349 1900 continue
350
351
352 ! do mass check final calc
353 if (mosaic_addemiss_masscheck > 0) call addemiss_masscheck( &
354 id, config_flags, 2, 'mosaic_ademiss', &
355 dtstep, efact1, dz8w, chem, chem_sum, &
356 ids,ide, jds,jde, kds,kde, &
357 ims,ime, jms,jme, kms,kme, &
358 its,ite, jts,jte, kts,kte, &
359 14, &
360 e_pm10,e_pm25,e_pm25i,e_pm25j,e_eci,e_ecj,e_orgi,e_orgj, &
361 e_so4j,e_so4c,e_no3j,e_no3c,e_orgc,e_ecc, &
362 e_ecc,e_ecc,e_ecc,e_ecc,e_ecc,e_ecc,e_ecc )
363
364
365 ! do seasalt emissions
366 if (mosaic_seasalt_emiss_active > 0) &
367 call mosaic_seasalt_emiss( &
368 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
369 ids,ide, jds,jde, kds,kde, &
370 ims,ime, jms,jme, kms,kme, &
371 its,ite, jts,jte, kts,kte )
372
373
374 return
375
376
377 END subroutine mosaic_addemiss
378
379
380
381 !----------------------------------------------------------------------
382 subroutine mosaic_seasalt_emiss( &
383 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
384 ids,ide, jds,jde, kds,kde, &
385 ims,ime, jms,jme, kms,kme, &
386 its,ite, jts,jte, kts,kte )
387 !
388 ! adds seasalt emissions for mosaic aerosol species
389 ! (i.e., seasalt emissions tendencies over time dtstep are applied
390 ! to the aerosol mixing ratios)
391 !
392
393 USE module_configure, only: grid_config_rec_type
394 USE module_state_description, only: num_chem, param_first_scalar
395 USE module_data_mosaic_asect
396
397 IMPLICIT NONE
398
399 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
400
401 INTEGER, INTENT(IN ) :: id, &
402 ids,ide, jds,jde, kds,kde, &
403 ims,ime, jms,jme, kms,kme, &
404 its,ite, jts,jte, kts,kte
405
406 REAL, INTENT(IN ) :: dtstep
407
408 ! 10-m wind speed components (m/s)
409 REAL, DIMENSION( ims:ime , jms:jme ), &
410 INTENT(IN ) :: u10, v10, xland
411
412 ! trace species mixing ratios (aerosol mass = ug/kg; number = #/kg)
413 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
414 INTENT(INOUT ) :: chem
415
416 ! alt = 1.0/(dry air density) in (m3/kg)
417 ! dz8w = layer thickness in (m)
418 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
419 INTENT(IN ) :: alt, dz8w
420
421 ! local variables
422 integer i, j, k, l, l_na, l_cl, n
423 integer iphase, itype
424 integer p1st
425
426 real dum, dumdlo, dumdhi, dumoceanfrac, dumspd10
427 real factaa, factbb, fracna, fraccl
428
429 real :: ssemfact_numb( maxd_asize, maxd_atype )
430 real :: ssemfact_mass( maxd_asize, maxd_atype )
431
432
433 ! for now just do itype=1
434 itype = 1
435 iphase = ai_phase
436
437 ! compute emissions factors for each size bin
438 ! (limit emissions to dp > 0.1 micrometer)
439 do n = 1, nsize_aer(itype)
440 dumdlo = max( dlo_sect(n,itype), 0.1e-4 )
441 dumdhi = max( dhi_sect(n,itype), 0.1e-4 )
442 call seasalt_emitfactors_1bin( 1, dumdlo, dumdhi, &
443 ssemfact_numb(n,itype), dum, ssemfact_mass(n,itype) )
444
445 ! convert mass emissions factor from (g/m2/s) to (ug/m2/s)
446 ssemfact_mass(n,itype) = ssemfact_mass(n,itype)*1.0e6
447 end do
448
449
450 ! loop over i,j and apply seasalt emissions
451 k = kts
452 do 1830 j = jts, jte
453 do 1820 i = its, ite
454
455 !Skip this point if over land. xland=1 for land and 2 for water.
456 !Also, there is no way to differentiate fresh from salt water.
457 !Currently, this assumes all water is salty.
458 if( xland(i,j) < 1.5 ) cycle
459
460 !wig: As far as I can tell, only real.exe knows the fractional breakdown
461 ! of land use. So, in wrf.exe, dumoceanfrac will always be 1.
462 dumoceanfrac = 1. !fraction of grid i,j that is salt water
463 dumspd10 = dumoceanfrac* &
464 ( (u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))**(0.5*3.41) )
465
466 ! factaa is (s*m2/kg-air)
467 ! factaa*ssemfact_mass(n) is (s*m2/kg-air)*(ug/m2/s) = ug/kg-air
468 ! factaa*ssemfact_numb(n) is (s*m2/kg-air)*( #/m2/s) = #/kg-air
469 factaa = (dtstep/dz8w(i,k,j))*alt(i,k,j)
470
471 factbb = factaa * dumspd10
472
473 ! apportion seasalt mass emissions assumming that seasalt is pure nacl
474 fracna = mw_na_aer / (mw_na_aer + mw_cl_aer)
475 fraccl = 1.0 - fracna
476
477 do 1810 n = 1, nsize_aer(itype)
478
479 ! only apply emissions if bin has both na and cl species
480 l_na = lptr_na_aer(n,itype,iphase)
481 l_cl = lptr_cl_aer(n,itype,iphase)
482 if ((l_na >= p1st) .and. (l_cl >= p1st)) then
483
484 chem(i,k,j,l_na) = chem(i,k,j,l_na) + &
485 factbb * ssemfact_mass(n,itype) * fracna
486
487 chem(i,k,j,l_cl) = chem(i,k,j,l_cl) + &
488 factbb * ssemfact_mass(n,itype) * fraccl
489
490 l = numptr_aer(n,itype,iphase)
491 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + &
492 factbb * ssemfact_numb(n,itype)
493
494 end if
495 1810 continue
496
497 1820 continue
498 1830 continue
499
500 return
501
502 END subroutine mosaic_seasalt_emiss
503
504
505
506 !c----------------------------------------------------------------------
507 !c following is from gong06b.f in
508 !c /net/cirrus/files1/home/rce/oldfiles1/box/seasaltg
509 !c----------------------------------------------------------------------
510 subroutine seasalt_emitfactors_1bin( ireduce_smallr_emit, &
511 dpdrylo_cm, dpdryhi_cm, &
512 emitfact_numb, emitfact_surf, emitfact_mass )
513 !c
514 !c computes seasalt emissions factors for a specifed
515 !c dry particle size range
516 !c dpdrylo_cm = lower dry diameter (cm)
517 !c dpdryhi_cm = upper dry diameter (cm)
518 !c
519 !c number and mass emissions are then computed as
520 !c number emissions (#/m2/s) == emitfact_numb * (spd10*3.41)
521 !c dry-sfc emissions (cm2/m2/s) == emitfact_surf * (spd10*3.41)
522 !c dry-mass emissions (g/m2/s) == emitfact_mass * (spd10*3.41)
523 !c
524 !c where spd10 = 10 m windspeed in m/s
525 !c
526 !c uses bubble emissions formula (eqn 5a) from
527 !c Gong et al. [JGR, 1997, p 3805-3818]
528 !c
529 !c *** for rdry < rdry_star, this formula overpredicts emissions.
530 !c A strictly ad hoc correction is applied to the formula,
531 !c based on sea-salt size measurements of
532 !c O'Dowd et al. [Atmos Environ, 1997, p 73-80]
533 !c
534 !c *** the correction is only applied when ireduce_smallr_emit > 0
535 !c
536 implicit none
537
538 !c subr arguments
539 integer ireduce_smallr_emit
540 real dpdrylo_cm, dpdryhi_cm, &
541 emitfact_numb, emitfact_surf, emitfact_mass
542
543 !c local variables
544 integer isub_bin, nsub_bin
545
546 real alnrdrylo
547 real drydens, drydens_43pi_em12, x_4pi_em8
548 real dum, dumadjust, dumb, dumexpb
549 real dumsum_na, dumsum_ma, dumsum_sa
550 real drwet, dlnrdry
551 real df0drwet, df0dlnrdry, df0dlnrdry_star
552 real relhum
553 real rdry, rdrylo, rdryhi, rdryaa, rdrybb
554 real rdrylowermost, rdryuppermost, rdry_star
555 real rwet, rwetaa, rwetbb
556 real rdry_cm, rwet_cm
557 real sigmag_star
558 real xmdry, xsdry
559
560 real pi
561 parameter (pi = 3.1415936536)
562
563 !c c1-c4 are constants for seasalt hygroscopic growth parameterization
564 !c in Eqn 3 and Table 2 of Gong et al. [1997]
565 real c1, c2, c3, c4, onethird
566 parameter (c1 = 0.7674)
567 parameter (c2 = 3.079)
568 parameter (c3 = 2.573e-11)
569 parameter (c4 = -1.424)
570 parameter (onethird = 1.0/3.0)
571
572
573 !c dry particle density (g/cm3)
574 drydens = 2.165
575 !c factor for radius (micrometers) to mass (g)
576 drydens_43pi_em12 = drydens*(4.0/3.0)*pi*1.0e-12
577 !c factor for radius (micrometers) to surface (cm2)
578 x_4pi_em8 = 4.0*pi*1.0e-8
579 !c bubble emissions formula assume 80% RH
580 relhum = 0.80
581
582 !c rdry_star = dry radius (micrometers) below which the
583 !c dF0/dr emission formula is adjusted downwards
584 rdry_star = 0.1
585 if (ireduce_smallr_emit .le. 0) rdry_star = -1.0e20
586 !c sigmag_star = geometric standard deviation used for
587 !c rdry < rdry_star
588 sigmag_star = 1.9
589
590 !c initialize sums
591 dumsum_na = 0.0
592 dumsum_sa = 0.0
593 dumsum_ma = 0.0
594
595 !c rdrylowermost, rdryuppermost = lower and upper
596 !c dry radii (micrometers) for overall integration
597 rdrylowermost = dpdrylo_cm*0.5e4
598 rdryuppermost = dpdryhi_cm*0.5e4
599
600 !c
601 !c "section 1"
602 !c integrate over rdry > rdry_star, where the dF0/dr emissions
603 !c formula is applicable
604 !c (when ireduce_smallr_emit <= 0, rdry_star = -1.0e20,
605 !c and the entire integration is done here)
606 !c
607 if (rdryuppermost .le. rdry_star) goto 2000
608
609 !c rdrylo, rdryhi = lower and upper dry radii (micrometers)
610 !c for this part of the integration
611 rdrylo = max( rdrylowermost, rdry_star )
612 rdryhi = rdryuppermost
613
614 nsub_bin = 1000
615
616 alnrdrylo = log( rdrylo )
617 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
618
619 !c compute rdry, rwet (micrometers) at lowest size
620 rdrybb = exp( alnrdrylo )
621 rdry_cm = rdrybb*1.0e-4
622 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
623 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
624 rwetbb = rwet_cm*1.0e4
625
626 do 1900 isub_bin = 1, nsub_bin
627
628 !c rdry, rwet at sub_bin lower boundary are those
629 !c at upper boundary of previous sub_bin
630 rdryaa = rdrybb
631 rwetaa = rwetbb
632
633 !c compute rdry, rwet (micrometers) at sub_bin upper boundary
634 dum = alnrdrylo + isub_bin*dlnrdry
635 rdrybb = exp( dum )
636
637 rdry_cm = rdrybb*1.0e-4
638 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
639 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
640 rwetbb = rwet_cm*1.0e4
641
642 !c geometric mean rdry, rwet (micrometers) for sub_bin
643 rdry = sqrt(rdryaa * rdrybb)
644 rwet = sqrt(rwetaa * rwetbb)
645 drwet = rwetbb - rwetaa
646
647 !c xmdry is dry mass in g
648 xmdry = drydens_43pi_em12 * (rdry**3.0)
649
650 !c xsdry is dry surface in cm2
651 xsdry = x_4pi_em8 * (rdry**2.0)
652
653 !c dumb is "B" in Gong's Eqn 5a
654 !c df0drwet is "dF0/dr" in Gong's Eqn 5a
655 dumb = ( 0.380 - log10(rwet) ) / 0.650
656 dumexpb = exp( -dumb*dumb)
657 df0drwet = 1.373 * (rwet**(-3.0)) * &
658 (1.0 + 0.057*(rwet**1.05)) * &
659 (10.0**(1.19*dumexpb))
660
661 dumsum_na = dumsum_na + drwet*df0drwet
662 dumsum_ma = dumsum_ma + drwet*df0drwet*xmdry
663 dumsum_sa = dumsum_sa + drwet*df0drwet*xsdry
664
665 1900 continue
666
667
668 !c
669 !c "section 2"
670 !c integrate over rdry < rdry_star, where the dF0/dr emissions
671 !c formula is just an extrapolation and predicts too many emissions
672 !c
673 !c 1. compute dF0/dln(rdry) = (dF0/drwet)*(drwet/dlnrdry)
674 !c at rdry_star
675 !c 2. for rdry < rdry_star, assume dF0/dln(rdry) is lognormal,
676 !c with the same lognormal parameters observed in
677 !c O'Dowd et al. [1997]
678 !c
679 2000 if (rdrylowermost .ge. rdry_star) goto 3000
680
681 !c compute dF0/dln(rdry) at rdry_star
682 rdryaa = 0.99*rdry_star
683 rdry_cm = rdryaa*1.0e-4
684 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
685 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
686 rwetaa = rwet_cm*1.0e4
687
688 rdrybb = 1.01*rdry_star
689 rdry_cm = rdrybb*1.0e-4
690 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
691 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
692 rwetbb = rwet_cm*1.0e4
693
694 rwet = 0.5*(rwetaa + rwetbb)
695 dumb = ( 0.380 - log10(rwet) ) / 0.650
696 dumexpb = exp( -dumb*dumb)
697 df0drwet = 1.373 * (rwet**(-3.0)) * &
698 (1.0 + 0.057*(rwet**1.05)) * &
699 (10.0**(1.19*dumexpb))
700
701 drwet = rwetbb - rwetaa
702 dlnrdry = log( rdrybb/rdryaa )
703 df0dlnrdry_star = df0drwet * (drwet/dlnrdry)
704
705
706 !c rdrylo, rdryhi = lower and upper dry radii (micrometers)
707 !c for this part of the integration
708 rdrylo = rdrylowermost
709 rdryhi = min( rdryuppermost, rdry_star )
710
711 nsub_bin = 1000
712
713 alnrdrylo = log( rdrylo )
714 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
715
716 do 2900 isub_bin = 1, nsub_bin
717
718 !c geometric mean rdry (micrometers) for sub_bin
719 dum = alnrdrylo + (isub_bin-0.5)*dlnrdry
720 rdry = exp( dum )
721
722 !c xmdry is dry mass in g
723 xmdry = drydens_43pi_em12 * (rdry**3.0)
724
725 !c xsdry is dry surface in cm2
726 xsdry = x_4pi_em8 * (rdry**2.0)
727
728 !c dumadjust is adjustment factor to reduce dF0/dr
729 dum = log( rdry/rdry_star ) / log( sigmag_star )
730 dumadjust = exp( -0.5*dum*dum )
731
732 df0dlnrdry = df0dlnrdry_star * dumadjust
733
734 dumsum_na = dumsum_na + dlnrdry*df0dlnrdry
735 dumsum_ma = dumsum_ma + dlnrdry*df0dlnrdry*xmdry
736 dumsum_sa = dumsum_sa + dlnrdry*df0dlnrdry*xsdry
737
738 2900 continue
739
740
741 !c
742 !c all done
743 !c
744 3000 emitfact_numb = dumsum_na
745 emitfact_mass = dumsum_ma
746 emitfact_surf = dumsum_sa
747
748 return
749 end subroutine seasalt_emitfactors_1bin
750
751
752
753 END MODULE module_mosaic_addemiss
754
755
756
757 !----------------------------------------------------------------------
758 subroutine addemiss_masscheck( id, config_flags, iflagaa, fromwhere, &
759 dtstep, efact1, dz8w, chem, chem_sum, &
760 ids,ide, jds,jde, kds,kde, &
761 ims,ime, jms,jme, kms,kme, &
762 its,ite, jts,jte, kts,kte, &
763 nemit, &
764 e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, &
765 e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21 )
766 !
767 ! produces test diagnostics for "addemiss" routines
768 !
769 ! 1. computes {sum over i,j,k ( chem * dz8w )} before and after
770 ! emissions tendencies are added to chem,
771 ! then prints (sum_after - sum_before)/(dtstep*efact1)
772 ! 2. computes {sum over i,j,k ( e_xxx )}, then prints them
773 ! the two should be equal
774 !
775
776 USE module_configure, only: grid_config_rec_type
777 USE module_state_description, only: num_chem
778
779 IMPLICIT NONE
780
781 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
782
783 INTEGER, INTENT(IN ) :: id, iflagaa, &
784 ids,ide, jds,jde, kds,kde, &
785 ims,ime, jms,jme, kms,kme, &
786 its,ite, jts,jte, kts,kte, &
787 nemit
788
789 REAL, INTENT(IN ) :: dtstep, efact1
790
791 ! trace species mixing ratios
792 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
793 INTENT(IN ) :: chem
794
795 ! trace species integrals
796 DOUBLE PRECISION, DIMENSION( num_chem ), &
797 INTENT(INOUT ) :: chem_sum
798
799 ! layer thickness (m)
800 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
801 INTENT(IN ) :: dz8w
802
803 ! emissions
804 ! REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
805 REAL, DIMENSION( ims:ime , kms:config_flags%kemit , jms:jme ), &
806 INTENT(IN ) :: &
807 e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, &
808 e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21
809
810 character(len=*), intent(in) :: fromwhere
811
812 ! local variables
813 integer, parameter :: nemit_maxd = 21
814 integer :: i, j, k, l
815 double precision :: chem_sum_prev
816 real :: fact
817 real :: emit_sum(nemit_maxd)
818
819
820 ! compute column integral, summed over i-j grids
821 ! compute {sum over i,j,k ( chem * dz8w ) }
822 ! on second pass (iflagaa==2), subtract the pass-one sum
823 do 1900 l = 1, num_chem
824
825 chem_sum_prev = chem_sum(l)
826 chem_sum(l) = 0.0
827
828 do j = jts, jte
829 do k = kts, kte-1
830 do i = its, ite
831 chem_sum(l) = chem_sum(l) + dble( chem(i,k,j,l)*dz8w(i,k,j) )
832 end do
833 end do
834 end do
835
836 if (iflagaa == 2) chem_sum(l) = (chem_sum(l) - chem_sum_prev)
837
838 1900 continue
839
840 if (iflagaa /= 2) return
841
842
843 ! compute {sum over i,j,k ( e_xxx ) }
844 emit_sum(:) = 0.0
845
846 do 2900 l = 1, min(nemit,nemit_maxd)
847 do j = jts, jte
848 do k = kts, min(config_flags%kemit,kte-1)
849 do i = its, ite
850 if (l== 1) emit_sum(l) = emit_sum(l) + e01(i,k,j)
851 if (l== 2) emit_sum(l) = emit_sum(l) + e02(i,k,j)
852 if (l== 3) emit_sum(l) = emit_sum(l) + e03(i,k,j)
853 if (l== 4) emit_sum(l) = emit_sum(l) + e04(i,k,j)
854 if (l== 5) emit_sum(l) = emit_sum(l) + e05(i,k,j)
855 if (l== 6) emit_sum(l) = emit_sum(l) + e06(i,k,j)
856 if (l== 7) emit_sum(l) = emit_sum(l) + e07(i,k,j)
857 if (l== 8) emit_sum(l) = emit_sum(l) + e08(i,k,j)
858 if (l== 9) emit_sum(l) = emit_sum(l) + e09(i,k,j)
859 if (l==10) emit_sum(l) = emit_sum(l) + e10(i,k,j)
860
861 if (l==11) emit_sum(l) = emit_sum(l) + e11(i,k,j)
862 if (l==12) emit_sum(l) = emit_sum(l) + e12(i,k,j)
863 if (l==13) emit_sum(l) = emit_sum(l) + e13(i,k,j)
864 if (l==14) emit_sum(l) = emit_sum(l) + e14(i,k,j)
865 if (l==15) emit_sum(l) = emit_sum(l) + e15(i,k,j)
866 if (l==16) emit_sum(l) = emit_sum(l) + e16(i,k,j)
867 if (l==17) emit_sum(l) = emit_sum(l) + e17(i,k,j)
868 if (l==18) emit_sum(l) = emit_sum(l) + e18(i,k,j)
869 if (l==19) emit_sum(l) = emit_sum(l) + e19(i,k,j)
870 if (l==20) emit_sum(l) = emit_sum(l) + e20(i,k,j)
871
872 if (l==21) emit_sum(l) = emit_sum(l) + e21(i,k,j)
873 end do
874 end do
875 end do
876 2900 continue
877
878 ! output the chem_sum and emit_sum
879 print 9110, fromwhere, its, ite, jts, jte
880 print 9100, 'chem_sum'
881 fact = 1.0/(dtstep*efact1)
882 print 9120, (l, fact*chem_sum(l), l=1,num_chem)
883 print 9100, 'emit_sum'
884 print 9120, (l, emit_sum(l), l=1,min(nemit,nemit_maxd))
885
886 9100 format( a )
887 9110 format( / 'addemiss_masscheck output, fromwhere = ', a / &
888 'its, ite, jts, jte =', 4i5 )
889 9120 format( 5( i5, 1pe11.3 ) )
890
891
892 return
893 END subroutine addemiss_masscheck
894