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 p1st = PARAM_FIRST_SCALAR
433
434 ! for now just do itype=1
435 itype = 1
436 iphase = ai_phase
437
438 ! compute emissions factors for each size bin
439 ! (limit emissions to dp > 0.1 micrometer)
440 do n = 1, nsize_aer(itype)
441 dumdlo = max( dlo_sect(n,itype), 0.1e-4 )
442 dumdhi = max( dhi_sect(n,itype), 0.1e-4 )
443 call seasalt_emitfactors_1bin( 1, dumdlo, dumdhi, &
444 ssemfact_numb(n,itype), dum, ssemfact_mass(n,itype) )
445
446 ! convert mass emissions factor from (g/m2/s) to (ug/m2/s)
447 ssemfact_mass(n,itype) = ssemfact_mass(n,itype)*1.0e6
448 end do
449
450
451 ! loop over i,j and apply seasalt emissions
452 k = kts
453 do 1830 j = jts, jte
454 do 1820 i = its, ite
455
456 !Skip this point if over land. xland=1 for land and 2 for water.
457 !Also, there is no way to differentiate fresh from salt water.
458 !Currently, this assumes all water is salty.
459 if( xland(i,j) < 1.5 ) cycle
460
461 !wig: As far as I can tell, only real.exe knows the fractional breakdown
462 ! of land use. So, in wrf.exe, dumoceanfrac will always be 1.
463 dumoceanfrac = 1. !fraction of grid i,j that is salt water
464 dumspd10 = dumoceanfrac* &
465 ( (u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))**(0.5*3.41) )
466
467 ! factaa is (s*m2/kg-air)
468 ! factaa*ssemfact_mass(n) is (s*m2/kg-air)*(ug/m2/s) = ug/kg-air
469 ! factaa*ssemfact_numb(n) is (s*m2/kg-air)*( #/m2/s) = #/kg-air
470 factaa = (dtstep/dz8w(i,k,j))*alt(i,k,j)
471
472 factbb = factaa * dumspd10
473
474 ! apportion seasalt mass emissions assumming that seasalt is pure nacl
475 fracna = mw_na_aer / (mw_na_aer + mw_cl_aer)
476 fraccl = 1.0 - fracna
477
478 do 1810 n = 1, nsize_aer(itype)
479
480 ! only apply emissions if bin has both na and cl species
481 l_na = lptr_na_aer(n,itype,iphase)
482 l_cl = lptr_cl_aer(n,itype,iphase)
483 if ((l_na >= p1st) .and. (l_cl >= p1st)) then
484
485 chem(i,k,j,l_na) = chem(i,k,j,l_na) + &
486 factbb * ssemfact_mass(n,itype) * fracna
487
488 chem(i,k,j,l_cl) = chem(i,k,j,l_cl) + &
489 factbb * ssemfact_mass(n,itype) * fraccl
490
491 l = numptr_aer(n,itype,iphase)
492 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + &
493 factbb * ssemfact_numb(n,itype)
494
495 end if
496 1810 continue
497
498 1820 continue
499 1830 continue
500
501 return
502
503 END subroutine mosaic_seasalt_emiss
504
505
506
507 !c----------------------------------------------------------------------
508 !c following is from gong06b.f in
509 !c /net/cirrus/files1/home/rce/oldfiles1/box/seasaltg
510 !c----------------------------------------------------------------------
511 subroutine seasalt_emitfactors_1bin( ireduce_smallr_emit, &
512 dpdrylo_cm, dpdryhi_cm, &
513 emitfact_numb, emitfact_surf, emitfact_mass )
514 !c
515 !c computes seasalt emissions factors for a specifed
516 !c dry particle size range
517 !c dpdrylo_cm = lower dry diameter (cm)
518 !c dpdryhi_cm = upper dry diameter (cm)
519 !c
520 !c number and mass emissions are then computed as
521 !c number emissions (#/m2/s) == emitfact_numb * (spd10*3.41)
522 !c dry-sfc emissions (cm2/m2/s) == emitfact_surf * (spd10*3.41)
523 !c dry-mass emissions (g/m2/s) == emitfact_mass * (spd10*3.41)
524 !c
525 !c where spd10 = 10 m windspeed in m/s
526 !c
527 !c uses bubble emissions formula (eqn 5a) from
528 !c Gong et al. [JGR, 1997, p 3805-3818]
529 !c
530 !c *** for rdry < rdry_star, this formula overpredicts emissions.
531 !c A strictly ad hoc correction is applied to the formula,
532 !c based on sea-salt size measurements of
533 !c O'Dowd et al. [Atmos Environ, 1997, p 73-80]
534 !c
535 !c *** the correction is only applied when ireduce_smallr_emit > 0
536 !c
537 implicit none
538
539 !c subr arguments
540 integer ireduce_smallr_emit
541 real dpdrylo_cm, dpdryhi_cm, &
542 emitfact_numb, emitfact_surf, emitfact_mass
543
544 !c local variables
545 integer isub_bin, nsub_bin
546
547 real alnrdrylo
548 real drydens, drydens_43pi_em12, x_4pi_em8
549 real dum, dumadjust, dumb, dumexpb
550 real dumsum_na, dumsum_ma, dumsum_sa
551 real drwet, dlnrdry
552 real df0drwet, df0dlnrdry, df0dlnrdry_star
553 real relhum
554 real rdry, rdrylo, rdryhi, rdryaa, rdrybb
555 real rdrylowermost, rdryuppermost, rdry_star
556 real rwet, rwetaa, rwetbb
557 real rdry_cm, rwet_cm
558 real sigmag_star
559 real xmdry, xsdry
560
561 real pi
562 parameter (pi = 3.1415936536)
563
564 !c c1-c4 are constants for seasalt hygroscopic growth parameterization
565 !c in Eqn 3 and Table 2 of Gong et al. [1997]
566 real c1, c2, c3, c4, onethird
567 parameter (c1 = 0.7674)
568 parameter (c2 = 3.079)
569 parameter (c3 = 2.573e-11)
570 parameter (c4 = -1.424)
571 parameter (onethird = 1.0/3.0)
572
573
574 !c dry particle density (g/cm3)
575 drydens = 2.165
576 !c factor for radius (micrometers) to mass (g)
577 drydens_43pi_em12 = drydens*(4.0/3.0)*pi*1.0e-12
578 !c factor for radius (micrometers) to surface (cm2)
579 x_4pi_em8 = 4.0*pi*1.0e-8
580 !c bubble emissions formula assume 80% RH
581 relhum = 0.80
582
583 !c rdry_star = dry radius (micrometers) below which the
584 !c dF0/dr emission formula is adjusted downwards
585 rdry_star = 0.1
586 if (ireduce_smallr_emit .le. 0) rdry_star = -1.0e20
587 !c sigmag_star = geometric standard deviation used for
588 !c rdry < rdry_star
589 sigmag_star = 1.9
590
591 !c initialize sums
592 dumsum_na = 0.0
593 dumsum_sa = 0.0
594 dumsum_ma = 0.0
595
596 !c rdrylowermost, rdryuppermost = lower and upper
597 !c dry radii (micrometers) for overall integration
598 rdrylowermost = dpdrylo_cm*0.5e4
599 rdryuppermost = dpdryhi_cm*0.5e4
600
601 !c
602 !c "section 1"
603 !c integrate over rdry > rdry_star, where the dF0/dr emissions
604 !c formula is applicable
605 !c (when ireduce_smallr_emit <= 0, rdry_star = -1.0e20,
606 !c and the entire integration is done here)
607 !c
608 if (rdryuppermost .le. rdry_star) goto 2000
609
610 !c rdrylo, rdryhi = lower and upper dry radii (micrometers)
611 !c for this part of the integration
612 rdrylo = max( rdrylowermost, rdry_star )
613 rdryhi = rdryuppermost
614
615 nsub_bin = 1000
616
617 alnrdrylo = log( rdrylo )
618 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
619
620 !c compute rdry, rwet (micrometers) at lowest size
621 rdrybb = exp( alnrdrylo )
622 rdry_cm = rdrybb*1.0e-4
623 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
624 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
625 rwetbb = rwet_cm*1.0e4
626
627 do 1900 isub_bin = 1, nsub_bin
628
629 !c rdry, rwet at sub_bin lower boundary are those
630 !c at upper boundary of previous sub_bin
631 rdryaa = rdrybb
632 rwetaa = rwetbb
633
634 !c compute rdry, rwet (micrometers) at sub_bin upper boundary
635 dum = alnrdrylo + isub_bin*dlnrdry
636 rdrybb = exp( dum )
637
638 rdry_cm = rdrybb*1.0e-4
639 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
640 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
641 rwetbb = rwet_cm*1.0e4
642
643 !c geometric mean rdry, rwet (micrometers) for sub_bin
644 rdry = sqrt(rdryaa * rdrybb)
645 rwet = sqrt(rwetaa * rwetbb)
646 drwet = rwetbb - rwetaa
647
648 !c xmdry is dry mass in g
649 xmdry = drydens_43pi_em12 * (rdry**3.0)
650
651 !c xsdry is dry surface in cm2
652 xsdry = x_4pi_em8 * (rdry**2.0)
653
654 !c dumb is "B" in Gong's Eqn 5a
655 !c df0drwet is "dF0/dr" in Gong's Eqn 5a
656 dumb = ( 0.380 - log10(rwet) ) / 0.650
657 dumexpb = exp( -dumb*dumb)
658 df0drwet = 1.373 * (rwet**(-3.0)) * &
659 (1.0 + 0.057*(rwet**1.05)) * &
660 (10.0**(1.19*dumexpb))
661
662 dumsum_na = dumsum_na + drwet*df0drwet
663 dumsum_ma = dumsum_ma + drwet*df0drwet*xmdry
664 dumsum_sa = dumsum_sa + drwet*df0drwet*xsdry
665
666 1900 continue
667
668
669 !c
670 !c "section 2"
671 !c integrate over rdry < rdry_star, where the dF0/dr emissions
672 !c formula is just an extrapolation and predicts too many emissions
673 !c
674 !c 1. compute dF0/dln(rdry) = (dF0/drwet)*(drwet/dlnrdry)
675 !c at rdry_star
676 !c 2. for rdry < rdry_star, assume dF0/dln(rdry) is lognormal,
677 !c with the same lognormal parameters observed in
678 !c O'Dowd et al. [1997]
679 !c
680 2000 if (rdrylowermost .ge. rdry_star) goto 3000
681
682 !c compute dF0/dln(rdry) at rdry_star
683 rdryaa = 0.99*rdry_star
684 rdry_cm = rdryaa*1.0e-4
685 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
686 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
687 rwetaa = rwet_cm*1.0e4
688
689 rdrybb = 1.01*rdry_star
690 rdry_cm = rdrybb*1.0e-4
691 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
692 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
693 rwetbb = rwet_cm*1.0e4
694
695 rwet = 0.5*(rwetaa + rwetbb)
696 dumb = ( 0.380 - log10(rwet) ) / 0.650
697 dumexpb = exp( -dumb*dumb)
698 df0drwet = 1.373 * (rwet**(-3.0)) * &
699 (1.0 + 0.057*(rwet**1.05)) * &
700 (10.0**(1.19*dumexpb))
701
702 drwet = rwetbb - rwetaa
703 dlnrdry = log( rdrybb/rdryaa )
704 df0dlnrdry_star = df0drwet * (drwet/dlnrdry)
705
706
707 !c rdrylo, rdryhi = lower and upper dry radii (micrometers)
708 !c for this part of the integration
709 rdrylo = rdrylowermost
710 rdryhi = min( rdryuppermost, rdry_star )
711
712 nsub_bin = 1000
713
714 alnrdrylo = log( rdrylo )
715 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
716
717 do 2900 isub_bin = 1, nsub_bin
718
719 !c geometric mean rdry (micrometers) for sub_bin
720 dum = alnrdrylo + (isub_bin-0.5)*dlnrdry
721 rdry = exp( dum )
722
723 !c xmdry is dry mass in g
724 xmdry = drydens_43pi_em12 * (rdry**3.0)
725
726 !c xsdry is dry surface in cm2
727 xsdry = x_4pi_em8 * (rdry**2.0)
728
729 !c dumadjust is adjustment factor to reduce dF0/dr
730 dum = log( rdry/rdry_star ) / log( sigmag_star )
731 dumadjust = exp( -0.5*dum*dum )
732
733 df0dlnrdry = df0dlnrdry_star * dumadjust
734
735 dumsum_na = dumsum_na + dlnrdry*df0dlnrdry
736 dumsum_ma = dumsum_ma + dlnrdry*df0dlnrdry*xmdry
737 dumsum_sa = dumsum_sa + dlnrdry*df0dlnrdry*xsdry
738
739 2900 continue
740
741
742 !c
743 !c all done
744 !c
745 3000 emitfact_numb = dumsum_na
746 emitfact_mass = dumsum_ma
747 emitfact_surf = dumsum_sa
748
749 return
750 end subroutine seasalt_emitfactors_1bin
751
752
753
754 END MODULE module_mosaic_addemiss
755
756
757
758 !----------------------------------------------------------------------
759 subroutine addemiss_masscheck( id, config_flags, iflagaa, fromwhere, &
760 dtstep, efact1, dz8w, chem, chem_sum, &
761 ids,ide, jds,jde, kds,kde, &
762 ims,ime, jms,jme, kms,kme, &
763 its,ite, jts,jte, kts,kte, &
764 nemit, &
765 e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, &
766 e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21 )
767 !
768 ! produces test diagnostics for "addemiss" routines
769 !
770 ! 1. computes {sum over i,j,k ( chem * dz8w )} before and after
771 ! emissions tendencies are added to chem,
772 ! then prints (sum_after - sum_before)/(dtstep*efact1)
773 ! 2. computes {sum over i,j,k ( e_xxx )}, then prints them
774 ! the two should be equal
775 !
776
777 USE module_configure, only: grid_config_rec_type
778 USE module_state_description, only: num_chem
779
780 IMPLICIT NONE
781
782 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
783
784 INTEGER, INTENT(IN ) :: id, iflagaa, &
785 ids,ide, jds,jde, kds,kde, &
786 ims,ime, jms,jme, kms,kme, &
787 its,ite, jts,jte, kts,kte, &
788 nemit
789
790 REAL, INTENT(IN ) :: dtstep, efact1
791
792 ! trace species mixing ratios
793 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
794 INTENT(IN ) :: chem
795
796 ! trace species integrals
797 DOUBLE PRECISION, DIMENSION( num_chem ), &
798 INTENT(INOUT ) :: chem_sum
799
800 ! layer thickness (m)
801 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
802 INTENT(IN ) :: dz8w
803
804 ! emissions
805 ! REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
806 REAL, DIMENSION( ims:ime , kms:config_flags%kemit , jms:jme ), &
807 INTENT(IN ) :: &
808 e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, &
809 e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21
810
811 character(len=*), intent(in) :: fromwhere
812
813 ! local variables
814 integer, parameter :: nemit_maxd = 21
815 integer :: i, j, k, l
816 double precision :: chem_sum_prev
817 real :: fact
818 real :: emit_sum(nemit_maxd)
819
820
821 ! compute column integral, summed over i-j grids
822 ! compute {sum over i,j,k ( chem * dz8w ) }
823 ! on second pass (iflagaa==2), subtract the pass-one sum
824 do 1900 l = 1, num_chem
825
826 chem_sum_prev = chem_sum(l)
827 chem_sum(l) = 0.0
828
829 do j = jts, jte
830 do k = kts, kte-1
831 do i = its, ite
832 chem_sum(l) = chem_sum(l) + dble( chem(i,k,j,l)*dz8w(i,k,j) )
833 end do
834 end do
835 end do
836
837 if (iflagaa == 2) chem_sum(l) = (chem_sum(l) - chem_sum_prev)
838
839 1900 continue
840
841 if (iflagaa /= 2) return
842
843
844 ! compute {sum over i,j,k ( e_xxx ) }
845 emit_sum(:) = 0.0
846
847 do 2900 l = 1, min(nemit,nemit_maxd)
848 do j = jts, jte
849 do k = kts, min(config_flags%kemit,kte-1)
850 do i = its, ite
851 if (l== 1) emit_sum(l) = emit_sum(l) + e01(i,k,j)
852 if (l== 2) emit_sum(l) = emit_sum(l) + e02(i,k,j)
853 if (l== 3) emit_sum(l) = emit_sum(l) + e03(i,k,j)
854 if (l== 4) emit_sum(l) = emit_sum(l) + e04(i,k,j)
855 if (l== 5) emit_sum(l) = emit_sum(l) + e05(i,k,j)
856 if (l== 6) emit_sum(l) = emit_sum(l) + e06(i,k,j)
857 if (l== 7) emit_sum(l) = emit_sum(l) + e07(i,k,j)
858 if (l== 8) emit_sum(l) = emit_sum(l) + e08(i,k,j)
859 if (l== 9) emit_sum(l) = emit_sum(l) + e09(i,k,j)
860 if (l==10) emit_sum(l) = emit_sum(l) + e10(i,k,j)
861
862 if (l==11) emit_sum(l) = emit_sum(l) + e11(i,k,j)
863 if (l==12) emit_sum(l) = emit_sum(l) + e12(i,k,j)
864 if (l==13) emit_sum(l) = emit_sum(l) + e13(i,k,j)
865 if (l==14) emit_sum(l) = emit_sum(l) + e14(i,k,j)
866 if (l==15) emit_sum(l) = emit_sum(l) + e15(i,k,j)
867 if (l==16) emit_sum(l) = emit_sum(l) + e16(i,k,j)
868 if (l==17) emit_sum(l) = emit_sum(l) + e17(i,k,j)
869 if (l==18) emit_sum(l) = emit_sum(l) + e18(i,k,j)
870 if (l==19) emit_sum(l) = emit_sum(l) + e19(i,k,j)
871 if (l==20) emit_sum(l) = emit_sum(l) + e20(i,k,j)
872
873 if (l==21) emit_sum(l) = emit_sum(l) + e21(i,k,j)
874 end do
875 end do
876 end do
877 2900 continue
878
879 ! output the chem_sum and emit_sum
880 print 9110, fromwhere, its, ite, jts, jte
881 print 9100, 'chem_sum'
882 fact = 1.0/(dtstep*efact1)
883 print 9120, (l, fact*chem_sum(l), l=1,num_chem)
884 print 9100, 'emit_sum'
885 print 9120, (l, emit_sum(l), l=1,min(nemit,nemit_maxd))
886
887 9100 format( a )
888 9110 format( / 'addemiss_masscheck output, fromwhere = ', a / &
889 'its, ite, jts, jte =', 4i5 )
890 9120 format( 5( i5, 1pe11.3 ) )
891
892
893 return
894 END subroutine addemiss_masscheck
895