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 CONTAINS
27
28
29
30 !----------------------------------------------------------------------
31 subroutine mosaic_addemiss( id, dtstep, u10, v10, alt, dz8w, xland, &
32 config_flags, chem, slai, ust, smois, ivgtyp, isltyp, &
33 e_pm10,e_pm25,e_pm25i,e_pm25j,e_eci,e_ecj,e_orgi,e_orgj, &
34 e_so4i,e_so4j,e_so4c,e_no3i,e_no3j,e_no3c,e_orgc,e_ecc, &
35 dust_emiss_active, seasalt_emiss_active, &
36 ids,ide, jds,jde, kds,kde, &
37 ims,ime, jms,jme, kms,kme, &
38 its,ite, jts,jte, kts,kte )
39 !
40 ! adds emissions for mosaic aerosol species
41 ! (i.e., emissions tendencies over time dtstep are applied
42 ! to the aerosol concentrations)
43 !
44
45 USE module_configure, only: grid_config_rec_type
46 USE module_state_description, only: num_chem, param_first_scalar, &
47 emiss_inpt_default, emiss_inpt_pnnl_rs, emiss_inpt_pnnl_cm
48 USE module_data_mosaic_asect
49
50 IMPLICIT NONE
51
52 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
53
54 INTEGER, INTENT(IN ) :: id, &
55 ids,ide, jds,jde, kds,kde, &
56 ims,ime, jms,jme, kms,kme, &
57 its,ite, jts,jte, kts,kte
58
59 INTEGER, INTENT(IN) :: dust_emiss_active, seasalt_emiss_active
60
61 REAL, INTENT(IN ) :: dtstep
62
63 ! 10-m wind speed components (m/s)
64 REAL, DIMENSION( ims:ime , jms:jme ) , &
65 INTENT(IN ) :: u10, v10, xland, slai, ust
66 INTEGER, DIMENSION( ims:ime , jms:jme ) , &
67 INTENT(IN ) :: ivgtyp, isltyp
68
69 ! trace species mixing ratios (aerosol mass = ug/kg-air; number = #/kg-air)
70 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
71 INTENT(INOUT ) :: chem
72 !
73 ! aerosol emissions arrays ((ug/m3)*m/s)
74 !
75 ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
76 REAL, DIMENSION( ims:ime, kms:config_flags%kemit, jms:jme ), &
77 INTENT(IN ) :: &
78 e_pm10,e_pm25,e_pm25i,e_pm25j,e_eci,e_ecj,e_orgi,e_orgj, &
79 e_so4i,e_so4j,e_so4c,e_no3i,e_no3j,e_no3c,e_orgc,e_ecc
80
81 ! 1/(dry air density) and layer thickness (m)
82 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
83 INTENT(IN ) :: alt, dz8w
84
85 REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) , &
86 INTENT(INOUT) :: smois
87
88 ! local variables
89 integer i, j, k, l, n
90 integer iphase, itype
91 integer p1st
92
93 real, parameter :: efact1 = 1.0
94 real :: aem_so4, aem_no3, aem_cl, aem_msa, aem_co3, aem_nh4, &
95 aem_na, aem_ca, aem_oin, aem_oc, aem_bc, aem_num
96 real dum, fact
97
98
99 ! fraction of sorgam i/aitken mode emissions that go to each
100 ! of the mosaic 8 "standard" sections
101 real, save :: fr8b_aem_sorgam_i(8) = &
102 (/ 0.965, 0.035, 0.000, 0.000, &
103 0.000, 0.000, 0.000, 0.000 /)
104
105 ! fraction of sorgam j/accum mode emissions that go to each
106 ! of the mosaic 8 "standard" sections
107 real, save :: fr8b_aem_sorgam_j(8) = &
108 (/ 0.026, 0.147, 0.350, 0.332, &
109 0.125, 0.019, 0.001, 0.000/)
110
111 ! fraction of sorgam coarse mode emissions that go to each
112 ! of the mosaic 8 "standard" sections
113 real, save :: fr8b_aem_sorgam_c(8) = &
114 (/ 0.000, 0.000, 0.000, 0.002, &
115 0.021, 0.110, 0.275, 0.592 /)
116
117 ! fraction of mosaic fine (< 2.5 um) emissions that go to each
118 ! of the mosaic 8 "standard" sections
119 !wig 1-Apr-2005, Updated fractional breakdown between bins. (See also
120 ! bdy_chem_value_mosaic and mosaic_init_wrf_mixrats_opt2
121 ! in module_mosaic_initmixrats.F.) Note that the values
122 ! here no longer match the other two subroutines.
123 !rce 10-may-2005, changed fr8b_aem_mosaic_f & fr8b_aem_mosaic_c
124 ! to values determined by jdf
125 real, save :: fr8b_aem_mosaic_f(8) = &
126 (/ 0.060, 0.045, 0.245, 0.400, 0.100, 0.150, 0., 0./) !10-may-2005
127 ! (/ 0.100, 0.045, 0.230, 0.375, 0.100, 0.150, 0., 0./) !1-Apr-2005 values
128 ! (/ 0.0275, 0.0426, 0.2303, 0.3885, 0.1100, 0.2011, 0., 0./) !15-Nov-2004 values
129 ! (/ 0.01, 0.05, 0.145, 0.60, 0.145, 0.05, 0.00, 0.00 /)
130 ! (/ 0.04, 0.10, 0.35, 0.29, 0.15, 0.07, 0.0, 0.0 /)
131
132 ! fraction of mosaic coarse (> 2.5 um) emissions that go to each
133 ! of the mosaic 8 "standard" sections
134 real, save :: fr8b_aem_mosaic_c(8) = &
135 (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.300, 0.700 /) ! 10-may-2005
136 ! (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.933, 0.067 /) ! as of apr-2005
137 ! (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.16, 0.84 /) ! "old"
138
139 ! following 5 arrays correspond to the above "fr8b_" arrays
140 ! but are set at run time base on input (namelist) parameters:
141 ! only the sorgam or mosaic arrays are non-zero, depending on
142 ! emiss_inpt_opt
143 ! when nsize_aer=4 (=number of sections), the values are
144 ! calculated by adding two of the 8-section values
145 real :: fr_aem_sorgam_i(8)
146 real :: fr_aem_sorgam_j(8)
147 real :: fr_aem_sorgam_c(8)
148 real :: fr_aem_mosaic_f(8)
149 real :: fr_aem_mosaic_c(8)
150 double precision :: chem_sum(num_chem)
151
152 character*80 msg
153
154
155 ! *** currently only works with ntype_aer = 1
156 itype = 1
157 iphase = ai_phase
158
159
160 !
161 ! compute factors used for apportioning either
162 ! the MADE-SORGAM emissions (i=aitken, j=accum, c=coarse modes) OR
163 ! the MOSAIC emission (f=fine (< 2.5 um), c=coarse (> 2.5 um))
164 ! to each size section
165 !
166 ! note: the fr8b_aer_xxxxxx_y values are specific to the mosaic 8 bin
167 ! structure with dlo_sect(1)=0.039 um and dhi_sect(8)=10.0 um,
168 ! also, the fr8b_aem_sorgam_y are specific for the assumed
169 ! dgvem_i/j/c used in the sorgam code
170 ! also, the fr8b_aem_mosaic_y values are specific for the assumed (by us)
171 ! size distribution for fine and coarse primary emissions
172 !
173 ! when there are 4 bins (nsize_aer=4), each of these "wider" bins
174 ! corresponds to 2 of the "narrower" bins of the 8 bin structure
175 !
176 ! note: if fr_aem_sorgam_y > 0, then fr_aem_mosaic_y = 0, and vice-versa
177 !
178 if ((nsize_aer(itype) .ne. 4) .and. (nsize_aer(itype) .ne. 8)) then
179 write(msg,'(a,i5)') &
180 'subr mosaic_addemiss - nsize_aer(itype) must be ' // &
181 '4 or 8 but = ', &
182 nsize_aer(itype)
183 call wrf_error_fatal( msg )
184 end if
185
186 fr_aem_sorgam_i(:) = 0.0
187 fr_aem_sorgam_j(:) = 0.0
188 fr_aem_sorgam_c(:) = 0.0
189 fr_aem_mosaic_f(:) = 0.0
190 fr_aem_mosaic_c(:) = 0.0
191
192 emiss_inpt_select_1: SELECT CASE( config_flags%emiss_inpt_opt )
193
194 CASE( emiss_inpt_default, emiss_inpt_pnnl_rs )
195 if (nsize_aer(itype) .eq. 8) then
196 fr_aem_sorgam_i(:) = fr8b_aem_sorgam_i(:)
197 fr_aem_sorgam_j(:) = fr8b_aem_sorgam_j(:)
198 fr_aem_sorgam_c(:) = fr8b_aem_sorgam_c(:)
199 else if (nsize_aer(itype) .eq. 4) then
200 do n = 1, nsize_aer(itype)
201 fr_aem_sorgam_i(n) = fr8b_aem_sorgam_i(2*n-1) &
202 + fr8b_aem_sorgam_i(2*n)
203 fr_aem_sorgam_j(n) = fr8b_aem_sorgam_j(2*n-1) &
204 + fr8b_aem_sorgam_j(2*n)
205 fr_aem_sorgam_c(n) = fr8b_aem_sorgam_c(2*n-1) &
206 + fr8b_aem_sorgam_c(2*n)
207 end do
208 end if
209
210 CASE( emiss_inpt_pnnl_cm )
211 if (nsize_aer(itype) .eq. 8) then
212 fr_aem_mosaic_f(:) = fr8b_aem_mosaic_f(:)
213 fr_aem_mosaic_c(:) = fr8b_aem_mosaic_c(:)
214 else if (nsize_aer(itype) .eq. 4) then
215 do n = 1, nsize_aer(itype)
216 fr_aem_mosaic_f(n) = fr8b_aem_mosaic_f(2*n-1) &
217 + fr8b_aem_mosaic_f(2*n)
218 fr_aem_mosaic_c(n) = fr8b_aem_mosaic_c(2*n-1) &
219 + fr8b_aem_mosaic_c(2*n)
220 end do
221 end if
222
223 CASE DEFAULT
224 return
225
226 END SELECT emiss_inpt_select_1
227
228 ! when mosaic_addemiss_active <= 0, set fr's to zero,
229 ! which causes the changes to chem(...) to be zero
230 if (mosaic_addemiss_active <= 0) then
231 fr_aem_sorgam_i(:) = 0.0
232 fr_aem_sorgam_j(:) = 0.0
233 fr_aem_sorgam_c(:) = 0.0
234 fr_aem_mosaic_f(:) = 0.0
235 fr_aem_mosaic_c(:) = 0.0
236 end if
237
238
239 ! do mass check initial calc
240 if (mosaic_addemiss_masscheck > 0) call addemiss_masscheck( &
241 id, config_flags, 1, 'mosaic_ademiss', &
242 dtstep, efact1, dz8w, chem, chem_sum, &
243 ids,ide, jds,jde, kds,kde, &
244 ims,ime, jms,jme, kms,kme, &
245 its,ite, jts,jte, kts,kte, &
246 14, &
247 e_pm10,e_pm25,e_pm25i,e_pm25j,e_eci,e_ecj,e_orgi,e_orgj, &
248 e_so4j,e_so4c,e_no3j,e_no3c,e_orgc,e_ecc, &
249 e_ecc,e_ecc,e_ecc,e_ecc,e_ecc,e_ecc,e_ecc )
250
251 p1st = param_first_scalar
252
253 !
254 ! apply emissions at each section and grid point
255 !
256 do 1900 n = 1, nsize_aer(itype)
257
258 do 1830 j = jts, jte
259 do 1820 k = kts, min(config_flags%kemit,kte)
260 do 1810 i = its, ite
261
262 ! compute mass emissions [(ug/m3)*m/s] for each species
263 ! using the apportioning fractions
264 aem_so4 = fr_aem_mosaic_f(n)*e_so4j(i,k,j) &
265 + fr_aem_mosaic_c(n)*e_so4c(i,k,j) &
266 + fr_aem_sorgam_i(n)*e_so4i(i,k,j) &
267 + fr_aem_sorgam_j(n)*e_so4j(i,k,j)
268
269 aem_no3 = fr_aem_mosaic_f(n)*e_no3j(i,k,j) &
270 + fr_aem_mosaic_c(n)*e_no3c(i,k,j) &
271 + fr_aem_sorgam_i(n)*e_no3i(i,k,j) &
272 + fr_aem_sorgam_j(n)*e_no3j(i,k,j)
273
274 aem_oc = fr_aem_mosaic_f(n)*e_orgj(i,k,j) &
275 + fr_aem_mosaic_c(n)*e_orgc(i,k,j) &
276 + fr_aem_sorgam_i(n)*e_orgi(i,k,j) &
277 + fr_aem_sorgam_j(n)*e_orgj(i,k,j)
278
279 aem_bc = fr_aem_mosaic_f(n)*e_ecj(i,k,j) &
280 + fr_aem_mosaic_c(n)*e_ecc(i,k,j) &
281 + fr_aem_sorgam_i(n)*e_eci(i,k,j) &
282 + fr_aem_sorgam_j(n)*e_ecj(i,k,j)
283
284 aem_oin = fr_aem_mosaic_f(n)*e_pm25j(i,k,j) &
285 + fr_aem_mosaic_c(n)*e_pm10(i,k,j) &
286 + fr_aem_sorgam_i(n)*e_pm25i(i,k,j) &
287 + fr_aem_sorgam_j(n)*e_pm25j(i,k,j) &
288 + fr_aem_sorgam_c(n)*e_pm10(i,k,j)
289
290 ! emissions for these species are currently zero
291 aem_nh4 = 0.0
292 aem_na = 0.0
293 aem_cl = 0.0
294 aem_ca = 0.0
295 aem_co3 = 0.0
296 aem_msa = 0.0
297
298 ! compute number emissions
299 ! first sum the mass-emissions/density
300 aem_num = &
301 (aem_so4/dens_so4_aer) + (aem_no3/dens_no3_aer) + &
302 (aem_cl /dens_cl_aer ) + (aem_msa/dens_msa_aer) + &
303 (aem_co3/dens_co3_aer) + (aem_nh4/dens_nh4_aer) + &
304 (aem_na /dens_na_aer ) + (aem_ca /dens_ca_aer ) + &
305 (aem_oin/dens_oin_aer) + (aem_oc /dens_oc_aer ) + &
306 (aem_bc /dens_bc_aer )
307
308 ! then multiply by 1.0e-6 to convert ug to g
309 ! and divide by particle volume at center of section (cm3)
310 aem_num = aem_num*1.0e-6/volumcen_sect(n,itype)
311
312 ! apply the emissions and convert from flux to mixing ratio
313 ! fact = (dtstep/dz8w(i,k,j))*(28.966/1000.)
314 fact = (dtstep/dz8w(i,k,j))*alt(i,k,j)
315
316 ! rce 22-nov-2004 - change to using the "..._aer" species pointers
317 l = lptr_so4_aer(n,itype,iphase)
318 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_so4*fact
319
320 l = lptr_no3_aer(n,itype,iphase)
321 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_no3*fact
322
323 l = lptr_cl_aer(n,itype,iphase)
324 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_cl*fact
325
326 l = lptr_msa_aer(n,itype,iphase)
327 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_msa*fact
328
329 l = lptr_co3_aer(n,itype,iphase)
330 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_co3*fact
331
332 l = lptr_nh4_aer(n,itype,iphase)
333 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_nh4*fact
334
335 l = lptr_na_aer(n,itype,iphase)
336 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_na*fact
337
338 l = lptr_ca_aer(n,itype,iphase)
339 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_ca*fact
340
341 l = lptr_oin_aer(n,itype,iphase)
342 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oin*fact
343
344 l = lptr_oc_aer(n,itype,iphase)
345 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oc*fact
346
347 l = lptr_bc_aer(n,itype,iphase)
348 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_bc*fact
349
350 l = numptr_aer(n,itype,iphase)
351 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_num*fact
352
353 1810 continue
354 1820 continue
355 1830 continue
356
357 1900 continue
358
359
360 ! do mass check final calc
361 if (mosaic_addemiss_masscheck > 0) call addemiss_masscheck( &
362 id, config_flags, 2, 'mosaic_ademiss', &
363 dtstep, efact1, dz8w, chem, chem_sum, &
364 ids,ide, jds,jde, kds,kde, &
365 ims,ime, jms,jme, kms,kme, &
366 its,ite, jts,jte, kts,kte, &
367 14, &
368 e_pm10,e_pm25,e_pm25i,e_pm25j,e_eci,e_ecj,e_orgi,e_orgj, &
369 e_so4j,e_so4c,e_no3j,e_no3c,e_orgc,e_ecc, &
370 e_ecc,e_ecc,e_ecc,e_ecc,e_ecc,e_ecc,e_ecc )
371
372
373 ! do seasalt emissions
374 if (seasalt_emiss_active > 0) &
375 call mosaic_seasalt_emiss( &
376 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
377 ids,ide, jds,jde, kds,kde, &
378 ims,ime, jms,jme, kms,kme, &
379 its,ite, jts,jte, kts,kte )
380
381 ! do dust emissions
382 if (dust_emiss_active > 0) &
383 call mosaic_dust_emiss( slai, ust, smois, ivgtyp, isltyp, &
384 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
385 ids,ide, jds,jde, kds,kde, &
386 ims,ime, jms,jme, kms,kme, &
387 its,ite, jts,jte, kts,kte )
388
389 return
390
391
392 END subroutine mosaic_addemiss
393
394
395
396 !----------------------------------------------------------------------
397 subroutine mosaic_seasalt_emiss( &
398 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
399 ids,ide, jds,jde, kds,kde, &
400 ims,ime, jms,jme, kms,kme, &
401 its,ite, jts,jte, kts,kte )
402 !
403 ! adds seasalt emissions for mosaic aerosol species
404 ! (i.e., seasalt emissions tendencies over time dtstep are applied
405 ! to the aerosol mixing ratios)
406 !
407
408 USE module_configure, only: grid_config_rec_type
409 USE module_state_description, only: num_chem, param_first_scalar
410 USE module_data_mosaic_asect
411
412 IMPLICIT NONE
413
414 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
415
416 INTEGER, INTENT(IN ) :: id, &
417 ids,ide, jds,jde, kds,kde, &
418 ims,ime, jms,jme, kms,kme, &
419 its,ite, jts,jte, kts,kte
420
421 REAL, INTENT(IN ) :: dtstep
422
423 ! 10-m wind speed components (m/s)
424 REAL, DIMENSION( ims:ime , jms:jme ), &
425 INTENT(IN ) :: u10, v10, xland
426
427 ! trace species mixing ratios (aerosol mass = ug/kg; number = #/kg)
428 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
429 INTENT(INOUT ) :: chem
430
431 ! alt = 1.0/(dry air density) in (m3/kg)
432 ! dz8w = layer thickness in (m)
433 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
434 INTENT(IN ) :: alt, dz8w
435
436 ! local variables
437 integer i, j, k, l, l_na, l_cl, n
438 integer iphase, itype
439 integer p1st
440
441 real dum, dumdlo, dumdhi, dumoceanfrac, dumspd10
442 real factaa, factbb, fracna, fraccl
443
444 real :: ssemfact_numb( maxd_asize, maxd_atype )
445 real :: ssemfact_mass( maxd_asize, maxd_atype )
446
447 p1st = PARAM_FIRST_SCALAR
448
449 ! for now just do itype=1
450 itype = 1
451 iphase = ai_phase
452
453 ! compute emissions factors for each size bin
454 ! (limit emissions to dp > 0.1 micrometer)
455 do n = 1, nsize_aer(itype)
456 dumdlo = max( dlo_sect(n,itype), 0.1e-4 )
457 dumdhi = max( dhi_sect(n,itype), 0.1e-4 )
458 call seasalt_emitfactors_1bin( 1, dumdlo, dumdhi, &
459 ssemfact_numb(n,itype), dum, ssemfact_mass(n,itype) )
460
461 ! convert mass emissions factor from (g/m2/s) to (ug/m2/s)
462 ssemfact_mass(n,itype) = ssemfact_mass(n,itype)*1.0e6
463 end do
464
465
466 ! loop over i,j and apply seasalt emissions
467 k = kts
468 do 1830 j = jts, jte
469 do 1820 i = its, ite
470
471 !Skip this point if over land. xland=1 for land and 2 for water.
472 !Also, there is no way to differentiate fresh from salt water.
473 !Currently, this assumes all water is salty.
474 if( xland(i,j) < 1.5 ) cycle
475
476 !wig: As far as I can tell, only real.exe knows the fractional breakdown
477 ! of land use. So, in wrf.exe, dumoceanfrac will always be 1.
478 dumoceanfrac = 1. !fraction of grid i,j that is salt water
479 dumspd10 = dumoceanfrac* &
480 ( (u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))**(0.5*3.41) )
481
482 ! factaa is (s*m2/kg-air)
483 ! factaa*ssemfact_mass(n) is (s*m2/kg-air)*(ug/m2/s) = ug/kg-air
484 ! factaa*ssemfact_numb(n) is (s*m2/kg-air)*( #/m2/s) = #/kg-air
485 factaa = (dtstep/dz8w(i,k,j))*alt(i,k,j)
486
487 factbb = factaa * dumspd10
488
489 ! apportion seasalt mass emissions assumming that seasalt is pure nacl
490 fracna = mw_na_aer / (mw_na_aer + mw_cl_aer)
491 fraccl = 1.0 - fracna
492
493 do 1810 n = 1, nsize_aer(itype)
494
495 ! only apply emissions if bin has both na and cl species
496 l_na = lptr_na_aer(n,itype,iphase)
497 l_cl = lptr_cl_aer(n,itype,iphase)
498 if ((l_na >= p1st) .and. (l_cl >= p1st)) then
499
500 chem(i,k,j,l_na) = chem(i,k,j,l_na) + &
501 factbb * ssemfact_mass(n,itype) * fracna
502
503 chem(i,k,j,l_cl) = chem(i,k,j,l_cl) + &
504 factbb * ssemfact_mass(n,itype) * fraccl
505
506 l = numptr_aer(n,itype,iphase)
507 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + &
508 factbb * ssemfact_numb(n,itype)
509 end if
510 1810 continue
511
512 1820 continue
513 1830 continue
514
515 return
516
517 END subroutine mosaic_seasalt_emiss
518
519
520
521 !c----------------------------------------------------------------------
522 !c following is from gong06b.f in
523 !c /net/cirrus/files1/home/rce/oldfiles1/box/seasaltg
524 !c----------------------------------------------------------------------
525 subroutine seasalt_emitfactors_1bin( ireduce_smallr_emit, &
526 dpdrylo_cm, dpdryhi_cm, &
527 emitfact_numb, emitfact_surf, emitfact_mass )
528 !c
529 !c computes seasalt emissions factors for a specifed
530 !c dry particle size range
531 !c dpdrylo_cm = lower dry diameter (cm)
532 !c dpdryhi_cm = upper dry diameter (cm)
533 !c
534 !c number and mass emissions are then computed as
535 !c number emissions (#/m2/s) == emitfact_numb * (spd10*3.41)
536 !c dry-sfc emissions (cm2/m2/s) == emitfact_surf * (spd10*3.41)
537 !c dry-mass emissions (g/m2/s) == emitfact_mass * (spd10*3.41)
538 !c
539 !c where spd10 = 10 m windspeed in m/s
540 !c
541 !c uses bubble emissions formula (eqn 5a) from
542 !c Gong et al. [JGR, 1997, p 3805-3818]
543 !c
544 !c *** for rdry < rdry_star, this formula overpredicts emissions.
545 !c A strictly ad hoc correction is applied to the formula,
546 !c based on sea-salt size measurements of
547 !c O'Dowd et al. [Atmos Environ, 1997, p 73-80]
548 !c
549 !c *** the correction is only applied when ireduce_smallr_emit > 0
550 !c
551 implicit none
552
553 !c subr arguments
554 integer ireduce_smallr_emit
555 real dpdrylo_cm, dpdryhi_cm, &
556 emitfact_numb, emitfact_surf, emitfact_mass
557
558 !c local variables
559 integer isub_bin, nsub_bin
560
561 real alnrdrylo
562 real drydens, drydens_43pi_em12, x_4pi_em8
563 real dum, dumadjust, dumb, dumexpb
564 real dumsum_na, dumsum_ma, dumsum_sa
565 real drwet, dlnrdry
566 real df0drwet, df0dlnrdry, df0dlnrdry_star
567 real relhum
568 real rdry, rdrylo, rdryhi, rdryaa, rdrybb
569 real rdrylowermost, rdryuppermost, rdry_star
570 real rwet, rwetaa, rwetbb
571 real rdry_cm, rwet_cm
572 real sigmag_star
573 real xmdry, xsdry
574
575 real pi
576 parameter (pi = 3.1415936536)
577
578 !c c1-c4 are constants for seasalt hygroscopic growth parameterization
579 !c in Eqn 3 and Table 2 of Gong et al. [1997]
580 real c1, c2, c3, c4, onethird
581 parameter (c1 = 0.7674)
582 parameter (c2 = 3.079)
583 parameter (c3 = 2.573e-11)
584 parameter (c4 = -1.424)
585 parameter (onethird = 1.0/3.0)
586
587
588 !c dry particle density (g/cm3)
589 drydens = 2.165
590 !c factor for radius (micrometers) to mass (g)
591 drydens_43pi_em12 = drydens*(4.0/3.0)*pi*1.0e-12
592 !c factor for radius (micrometers) to surface (cm2)
593 x_4pi_em8 = 4.0*pi*1.0e-8
594 !c bubble emissions formula assume 80% RH
595 relhum = 0.80
596
597 !c rdry_star = dry radius (micrometers) below which the
598 !c dF0/dr emission formula is adjusted downwards
599 rdry_star = 0.1
600 if (ireduce_smallr_emit .le. 0) rdry_star = -1.0e20
601 !c sigmag_star = geometric standard deviation used for
602 !c rdry < rdry_star
603 sigmag_star = 1.9
604
605 !c initialize sums
606 dumsum_na = 0.0
607 dumsum_sa = 0.0
608 dumsum_ma = 0.0
609
610 !c rdrylowermost, rdryuppermost = lower and upper
611 !c dry radii (micrometers) for overall integration
612 rdrylowermost = dpdrylo_cm*0.5e4
613 rdryuppermost = dpdryhi_cm*0.5e4
614
615 !c
616 !c "section 1"
617 !c integrate over rdry > rdry_star, where the dF0/dr emissions
618 !c formula is applicable
619 !c (when ireduce_smallr_emit <= 0, rdry_star = -1.0e20,
620 !c and the entire integration is done here)
621 !c
622 if (rdryuppermost .le. rdry_star) goto 2000
623
624 !c rdrylo, rdryhi = lower and upper dry radii (micrometers)
625 !c for this part of the integration
626 rdrylo = max( rdrylowermost, rdry_star )
627 rdryhi = rdryuppermost
628
629 nsub_bin = 1000
630
631 alnrdrylo = log( rdrylo )
632 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
633
634 !c compute rdry, rwet (micrometers) at lowest size
635 rdrybb = exp( alnrdrylo )
636 rdry_cm = rdrybb*1.0e-4
637 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
638 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
639 rwetbb = rwet_cm*1.0e4
640
641 do 1900 isub_bin = 1, nsub_bin
642
643 !c rdry, rwet at sub_bin lower boundary are those
644 !c at upper boundary of previous sub_bin
645 rdryaa = rdrybb
646 rwetaa = rwetbb
647
648 !c compute rdry, rwet (micrometers) at sub_bin upper boundary
649 dum = alnrdrylo + isub_bin*dlnrdry
650 rdrybb = exp( dum )
651
652 rdry_cm = rdrybb*1.0e-4
653 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
654 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
655 rwetbb = rwet_cm*1.0e4
656
657 !c geometric mean rdry, rwet (micrometers) for sub_bin
658 rdry = sqrt(rdryaa * rdrybb)
659 rwet = sqrt(rwetaa * rwetbb)
660 drwet = rwetbb - rwetaa
661
662 !c xmdry is dry mass in g
663 xmdry = drydens_43pi_em12 * (rdry**3.0)
664
665 !c xsdry is dry surface in cm2
666 xsdry = x_4pi_em8 * (rdry**2.0)
667
668 !c dumb is "B" in Gong's Eqn 5a
669 !c df0drwet is "dF0/dr" in Gong's Eqn 5a
670 dumb = ( 0.380 - log10(rwet) ) / 0.650
671 dumexpb = exp( -dumb*dumb)
672 df0drwet = 1.373 * (rwet**(-3.0)) * &
673 (1.0 + 0.057*(rwet**1.05)) * &
674 (10.0**(1.19*dumexpb))
675
676 dumsum_na = dumsum_na + drwet*df0drwet
677 dumsum_ma = dumsum_ma + drwet*df0drwet*xmdry
678 dumsum_sa = dumsum_sa + drwet*df0drwet*xsdry
679
680 1900 continue
681
682
683 !c
684 !c "section 2"
685 !c integrate over rdry < rdry_star, where the dF0/dr emissions
686 !c formula is just an extrapolation and predicts too many emissions
687 !c
688 !c 1. compute dF0/dln(rdry) = (dF0/drwet)*(drwet/dlnrdry)
689 !c at rdry_star
690 !c 2. for rdry < rdry_star, assume dF0/dln(rdry) is lognormal,
691 !c with the same lognormal parameters observed in
692 !c O'Dowd et al. [1997]
693 !c
694 2000 if (rdrylowermost .ge. rdry_star) goto 3000
695
696 !c compute dF0/dln(rdry) at rdry_star
697 rdryaa = 0.99*rdry_star
698 rdry_cm = rdryaa*1.0e-4
699 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
700 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
701 rwetaa = rwet_cm*1.0e4
702
703 rdrybb = 1.01*rdry_star
704 rdry_cm = rdrybb*1.0e-4
705 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
706 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
707 rwetbb = rwet_cm*1.0e4
708
709 rwet = 0.5*(rwetaa + rwetbb)
710 dumb = ( 0.380 - log10(rwet) ) / 0.650
711 dumexpb = exp( -dumb*dumb)
712 df0drwet = 1.373 * (rwet**(-3.0)) * &
713 (1.0 + 0.057*(rwet**1.05)) * &
714 (10.0**(1.19*dumexpb))
715
716 drwet = rwetbb - rwetaa
717 dlnrdry = log( rdrybb/rdryaa )
718 df0dlnrdry_star = df0drwet * (drwet/dlnrdry)
719
720
721 !c rdrylo, rdryhi = lower and upper dry radii (micrometers)
722 !c for this part of the integration
723 rdrylo = rdrylowermost
724 rdryhi = min( rdryuppermost, rdry_star )
725
726 nsub_bin = 1000
727
728 alnrdrylo = log( rdrylo )
729 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
730
731 do 2900 isub_bin = 1, nsub_bin
732
733 !c geometric mean rdry (micrometers) for sub_bin
734 dum = alnrdrylo + (isub_bin-0.5)*dlnrdry
735 rdry = exp( dum )
736
737 !c xmdry is dry mass in g
738 xmdry = drydens_43pi_em12 * (rdry**3.0)
739
740 !c xsdry is dry surface in cm2
741 xsdry = x_4pi_em8 * (rdry**2.0)
742
743 !c dumadjust is adjustment factor to reduce dF0/dr
744 dum = log( rdry/rdry_star ) / log( sigmag_star )
745 dumadjust = exp( -0.5*dum*dum )
746
747 df0dlnrdry = df0dlnrdry_star * dumadjust
748
749 dumsum_na = dumsum_na + dlnrdry*df0dlnrdry
750 dumsum_ma = dumsum_ma + dlnrdry*df0dlnrdry*xmdry
751 dumsum_sa = dumsum_sa + dlnrdry*df0dlnrdry*xsdry
752
753 2900 continue
754
755
756 !c
757 !c all done
758 !c
759 3000 emitfact_numb = dumsum_na
760 emitfact_mass = dumsum_ma
761 emitfact_surf = dumsum_sa
762
763 return
764 end subroutine seasalt_emitfactors_1bin
765
766
767
768 END MODULE module_mosaic_addemiss
769
770 !----------------------------------------------------------------------
771 subroutine mosaic_dust_emiss( slai,ust, smois, ivgtyp, isltyp, &
772 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
773 ids,ide, jds,jde, kds,kde, &
774 ims,ime, jms,jme, kms,kme, &
775 its,ite, jts,jte, kts,kte )
776 !
777 ! adds dust emissions for mosaic aerosol species
778 ! This is a simple dust scheme at this point, based more on Mahowald et al., 1999
779 ! Many things are hard coded in to this first version (e.g. 8 aerosol bins)
780 !
781 ! (i.e., emissions tendencies over time dtstep are applied
782 ! to the aerosol mixing ratios)
783 !
784
785 USE module_configure, only: grid_config_rec_type
786 USE module_state_description, only: num_chem, param_first_scalar
787 USE module_data_mosaic_asect
788
789 IMPLICIT NONE
790
791 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
792
793 INTEGER, INTENT(IN ) :: id, &
794 ids,ide, jds,jde, kds,kde, &
795 ims,ime, jms,jme, kms,kme, &
796 its,ite, jts,jte, kts,kte
797
798 REAL, INTENT(IN ) :: dtstep
799
800 ! 10-m wind speed components (m/s)
801 REAL, DIMENSION( ims:ime , jms:jme ), &
802 INTENT(IN ) :: u10, v10, xland, slai, ust
803 INTEGER, DIMENSION( ims:ime , jms:jme ), &
804 INTENT(IN ) :: ivgtyp, isltyp
805
806 ! trace species mixing ratios (aerosol mass = ug/kg; number = #/kg)
807 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
808 INTENT(INOUT ) :: chem
809
810 ! alt = 1.0/(dry air density) in (m3/kg)
811 ! dz8w = layer thickness in (m)
812 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
813 INTENT(IN ) :: alt, dz8w
814
815 REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) , &
816 INTENT(INOUT) :: smois
817
818 ! local variables
819 integer i, j, k, l, l_oin, l_ca, l_co3, n, ii
820 integer iphase, itype, izob
821 integer p1st
822
823 real dum, dumdlo, dumdhi, dumlandfrac, dumspd10
824 real factaa, factbb, fracoin, fracca, fracco3, fractot
825 real ustart, ustar1, ustart0
826 real alphamask, f8, f50, f51, f52, wetfactor, sumdelta, ftot
827 real smois_grav, wp, pclay
828 real :: beta(4,7)
829 real :: gamma(4), delta(4)
830 real :: sz(8)
831 real :: dustflux, densdust, mass1part
832 !
833 ! from Nickovic et al., JGR, 2001 and Shaw et al. 2007
834 ! beta: fraction of clay, small silt, large silt, and sand correcsponding to Zobler class (7)
835 ! beta (1,*) for 0.5-1 um
836 ! beta (2,*) for 1-10 um
837 ! beta (3,*) for 10-25 um
838 ! beta (4,*) for 25-50 um
839 !
840 beta(1,1)=0.12
841 beta(2,1)=0.04
842 beta(3,1)=0.04
843 beta(4,1)=0.80
844 beta(1,2)=0.34
845 beta(2,2)=0.28
846 beta(3,2)=0.28
847 beta(4,2)=0.10
848 beta(1,3)=0.45
849 beta(2,3)=0.15
850 beta(3,3)=0.15
851 beta(4,3)=0.25
852 beta(1,4)=0.12
853 beta(2,4)=0.09
854 beta(3,4)=0.09
855 beta(4,4)=0.70
856 beta(1,5)=0.40
857 beta(2,5)=0.05
858 beta(3,5)=0.05
859 beta(4,5)=0.50
860 beta(1,6)=0.34
861 beta(2,6)=0.18
862 beta(3,6)=0.18
863 beta(4,6)=0.30
864 beta(1,7)=0.22
865 beta(2,7)=0.09
866 beta(3,7)=0.09
867 beta(4,7)=0.60
868 gamma(1)=0.08
869 gamma(2)=1.00
870 gamma(3)=1.00
871 gamma(4)=0.12
872 !
873 ! * Mass fractions for each size bin. These values were recommended by
874 ! Natalie Mahowold, with bins 7 and 8 the same as bins 3 and 4 from CAM.
875 ! * Changed slightly since Natelie's estimates do not add up to 1.0
876 ! * This would need to be made more generic for other bin sizes.
877 ! sz(1)=0
878 ! sz(2)=1.78751e-06
879 ! sz(3)=0.000273786
880 ! sz(4)=0.00847978
881 ! sz(5)=0.056055
882 ! sz(6)=0.0951896
883 ! sz(7)=0.17
884 ! sz(8)=0.67
885 sz(1)=0.0
886 sz(2)=0.0
887 sz(3)=0.0005
888 sz(4)=0.0095
889 sz(5)=0.03
890 sz(6)=0.10
891 sz(7)=0.18
892 sz(8)=0.68
893
894 ! for now just do itype=1
895 itype = 1
896 iphase = ai_phase
897
898 ! loop over i,j and apply dust emissions
899 k = kts
900 do 1830 j = jts, jte
901 do 1820 i = its, ite
902
903 if( xland(i,j) > 1.5 ) cycle
904
905 ! compute wind speed anyway, even though ustar is used below
906
907 dumlandfrac = 1.
908 dumspd10=(u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))**(0.5)
909 if(dumspd10 >= 5.0) then
910 dumspd10 = dumlandfrac* &
911 ( dumspd10*dumspd10*(dumspd10-5.0))
912 else
913 dumspd10=0.
914 endif
915
916 ! The parameters are based on Shaw et al. (2007) to be submitted to Atmos. Environ.
917
918 ! part1 - compute vegetation mask
919 !
920 ! * f8, f50, f51, f52 refer to vegetation classes from the Olsen categories
921 ! for desert, sand desert, grass aemi-desert, and shrub semi-desert
922 ! * in WRF, vegetation types 7, 8 and 10 are grassland, shrubland, and savanna
923 ! that are dominate types in Mexico and probably have some erodable surface
924 ! during the dry season
925 ! * currently modified these values so that only a small fraction of cell
926 ! area is erodable
927 ! * these values are highly tuneable!
928
929 alphamask=0.001
930 if (ivgtyp(i,j) .eq. 7) then
931 ! f8=0.05
932 ! f50=0.05
933 ! f51=0.70
934 ! f52=0.00
935 f8=0.005
936 f50=0.00
937 f51=0.10
938 f52=0.00
939 alphamask=(f8+f50)*1.0+(f51+f52)*0.5
940 endif
941 if (ivgtyp(i,j) .eq. 8) then
942 ! f8=0.10
943 ! f50=0.05
944 ! f51=0.00
945 ! f52=0.85
946 f8=0.010
947 f50=0.00
948 f51=0.00
949 f52=0.15
950 alphamask=(f8+f50)*1.0+(f51+f52)*0.5
951 endif
952 if (ivgtyp(i,j) .eq. 10) then
953 ! f8=0.05
954 ! f50=0.05
955 ! f51=0.50
956 ! f52=0.00
957 f8=0.00
958 f50=0.00
959 f51=0.01
960 f52=0.00
961 alphamask=(f8+f50)*1.0+(f51+f52)*0.5
962 endif
963
964 ! part2 - zobler
965 !
966 ! * in Shaw's paper, dust is computed for 4 size ranges:
967 ! 0.5-1 um
968 ! 1-10 um
969 ! 10-25 um
970 ! 25-50 um
971 ! * Shaw's paper also accounts for sub-grid variability in soil
972 ! texture, but here we just assume the same soil texture for each
973 ! grid cell
974 ! * since MOSAIC is currently has a maximum size range up to 10 um,
975 ! neglect upper 2 size ranges and lowest size range (assume small)
976 ! * map WRF soil classes arbitrarily to Zolber soil textural classes
977 ! * skip dust computations for WRF soil classes greater than 13, i.e.
978 ! do not compute dust over water, bedrock, and other surfaces
979 ! * should be skipping for water surface at this point anyway
980 !
981 izob=0
982 if(isltyp(i,j).eq.1) izob=1
983 if(isltyp(i,j).eq.2) izob=1
984 if(isltyp(i,j).eq.3) izob=4
985 if(isltyp(i,j).eq.4) izob=2
986 if(isltyp(i,j).eq.5) izob=2
987 if(isltyp(i,j).eq.6) izob=2
988 if(isltyp(i,j).eq.7) izob=7
989 if(isltyp(i,j).eq.8) izob=2
990 if(isltyp(i,j).eq.9) izob=6
991 if(isltyp(i,j).eq.10) izob=5
992 if(isltyp(i,j).eq.11) izob=2
993 if(isltyp(i,j).eq.12) izob=3
994 if(isltyp(i,j).ge.13) izob=0
995 if(izob.eq.0) goto 1840
996 beta = 0.28
997 !
998 ! part3 - dustprod
999 !
1000 do ii=1,4
1001 delta(ii)=0.0
1002 enddo
1003 sumdelta=0.0
1004 do ii=1,4
1005 delta(ii)=beta(ii,izob)*gamma(ii)
1006 if(ii.lt.4) then
1007 sumdelta=sumdelta+delta(ii)
1008 endif
1009 enddo
1010 do ii=1,4
1011 delta(ii)=delta(ii)/sumdelta
1012 enddo
1013
1014 ! part4 - wetness
1015 !
1016 ! * assume dry for now, have passed in soil moisture to this routine
1017 ! but needs to be included here
1018 ! * wetfactor less than 1 would reduce dustflux
1019 ! * convert model soil moisture (m3/m3) to gravimetric soil moisture
1020 ! (mass of water / mass of soil in %) assuming a constant density
1021 ! for soil
1022 pclay=beta(1,izob)*100.
1023 wp=0.0014*pclay*pclay+0.17*pclay
1024 smois_grav=(smois(i,1,j)/2.6)*100.
1025 if(smois_grav.gt.wp) then
1026 wetfactor=sqrt(1.0+1.21*(smois_grav-wp)**0.68)
1027 else
1028 wetfactor=1.0
1029 endif
1030 ! wetfactor=1.0
1031
1032 ! part5 - dustflux
1033
1034 ustar1=ust(i,j)*100.0
1035 ustart0=20.0
1036 ustart=ustart0*wetfactor
1037 if(ustar1.le.ustart) then
1038 dustflux=0.0
1039 else
1040 dustflux=1.0e-14*(ustar1**4)*(1.0-(ustart/ustar1))
1041 endif
1042 dustflux=dustflux*10.0
1043 ! units kg m-2 s-1
1044 ftot=0.0
1045 do ii=1,2
1046 ftot=ftot+dustflux*alphamask*delta(ii)
1047 enddo
1048 ! convert to ug m-2 s-1
1049 ftot=ftot*1.0e+09
1050
1051 ! apportion other inorganics only
1052 factaa = (dtstep/dz8w(i,k,j))*alt(i,k,j)
1053 factbb = factaa * ftot
1054 fracoin = 0.97
1055 fracca = 0.03*0.4
1056 fracco3 = 0.03*0.6
1057 fractot = fracoin + fracca + fracco3
1058 ! if (ivgtyp(i,j) .eq. 8) print*,'jdf',i,j,ustar1,ustart0,factaa,ftot
1059 do 1810 n = 1, nsize_aer(itype)
1060 l_oin = lptr_oin_aer(n,itype,iphase)
1061 l_ca = lptr_ca_aer(n,itype,iphase)
1062 l_co3 = lptr_co3_aer(n,itype,iphase)
1063 if (l_oin >= p1st) then
1064 chem(i,k,j,l_oin) = chem(i,k,j,l_oin) + &
1065 factbb * sz(n) * fracoin
1066 if (l_ca >= p1st) &
1067 chem(i,k,j,l_ca) = chem(i,k,j,l_ca) + &
1068 factbb * sz(n) * fracca
1069 if (l_co3 >= p1st) &
1070 chem(i,k,j,l_co3) = chem(i,k,j,l_co3) + &
1071 factbb * sz(n) * fracco3
1072 ! mass1part is mass of a single particle in ug, density of dust ~2.5 g cm-3
1073 densdust=2.5
1074 mass1part=0.523598*(dcen_sect(n,itype)**3)*densdust*1.0e06
1075 l = numptr_aer(n,itype,iphase)
1076 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + &
1077 factbb * sz(n) * fractot / mass1part
1078 end if
1079 1810 continue
1080
1081 1840 continue
1082
1083 1820 continue
1084 1830 continue
1085
1086 return
1087
1088 END subroutine mosaic_dust_emiss
1089
1090
1091 !----------------------------------------------------------------------
1092 subroutine addemiss_masscheck( id, config_flags, iflagaa, fromwhere, &
1093 dtstep, efact1, dz8w, chem, chem_sum, &
1094 ids,ide, jds,jde, kds,kde, &
1095 ims,ime, jms,jme, kms,kme, &
1096 its,ite, jts,jte, kts,kte, &
1097 nemit, &
1098 e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, &
1099 e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21 )
1100 !
1101 ! produces test diagnostics for "addemiss" routines
1102 !
1103 ! 1. computes {sum over i,j,k ( chem * dz8w )} before and after
1104 ! emissions tendencies are added to chem,
1105 ! then prints (sum_after - sum_before)/(dtstep*efact1)
1106 ! 2. computes {sum over i,j,k ( e_xxx )}, then prints them
1107 ! the two should be equal
1108 !
1109
1110 USE module_configure, only: grid_config_rec_type
1111 USE module_state_description, only: num_chem
1112
1113 IMPLICIT NONE
1114
1115 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1116
1117 INTEGER, INTENT(IN ) :: id, iflagaa, &
1118 ids,ide, jds,jde, kds,kde, &
1119 ims,ime, jms,jme, kms,kme, &
1120 its,ite, jts,jte, kts,kte, &
1121 nemit
1122
1123 REAL, INTENT(IN ) :: dtstep, efact1
1124
1125 ! trace species mixing ratios
1126 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
1127 INTENT(IN ) :: chem
1128
1129 ! trace species integrals
1130 DOUBLE PRECISION, DIMENSION( num_chem ), &
1131 INTENT(INOUT ) :: chem_sum
1132
1133 ! layer thickness (m)
1134 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
1135 INTENT(IN ) :: dz8w
1136
1137 ! emissions
1138 ! REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
1139 REAL, DIMENSION( ims:ime , kms:config_flags%kemit , jms:jme ), &
1140 INTENT(IN ) :: &
1141 e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, &
1142 e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21
1143
1144 character(len=*), intent(in) :: fromwhere
1145
1146 ! local variables
1147 integer, parameter :: nemit_maxd = 21
1148 integer :: i, j, k, l
1149 double precision :: chem_sum_prev
1150 real :: fact
1151 real :: emit_sum(nemit_maxd)
1152
1153
1154 ! compute column integral, summed over i-j grids
1155 ! compute {sum over i,j,k ( chem * dz8w ) }
1156 ! on second pass (iflagaa==2), subtract the pass-one sum
1157 do 1900 l = 1, num_chem
1158
1159 chem_sum_prev = chem_sum(l)
1160 chem_sum(l) = 0.0
1161
1162 do j = jts, jte
1163 do k = kts, kte
1164 do i = its, ite
1165 chem_sum(l) = chem_sum(l) + dble( chem(i,k,j,l)*dz8w(i,k,j) )
1166 end do
1167 end do
1168 end do
1169
1170 if (iflagaa == 2) chem_sum(l) = (chem_sum(l) - chem_sum_prev)
1171
1172 1900 continue
1173
1174 if (iflagaa /= 2) return
1175
1176
1177 ! compute {sum over i,j,k ( e_xxx ) }
1178 emit_sum(:) = 0.0
1179
1180 do 2900 l = 1, min(nemit,nemit_maxd)
1181 do j = jts, jte
1182 do k = kts, min(config_flags%kemit,kte)
1183 do i = its, ite
1184 if (l== 1) emit_sum(l) = emit_sum(l) + e01(i,k,j)
1185 if (l== 2) emit_sum(l) = emit_sum(l) + e02(i,k,j)
1186 if (l== 3) emit_sum(l) = emit_sum(l) + e03(i,k,j)
1187 if (l== 4) emit_sum(l) = emit_sum(l) + e04(i,k,j)
1188 if (l== 5) emit_sum(l) = emit_sum(l) + e05(i,k,j)
1189 if (l== 6) emit_sum(l) = emit_sum(l) + e06(i,k,j)
1190 if (l== 7) emit_sum(l) = emit_sum(l) + e07(i,k,j)
1191 if (l== 8) emit_sum(l) = emit_sum(l) + e08(i,k,j)
1192 if (l== 9) emit_sum(l) = emit_sum(l) + e09(i,k,j)
1193 if (l==10) emit_sum(l) = emit_sum(l) + e10(i,k,j)
1194
1195 if (l==11) emit_sum(l) = emit_sum(l) + e11(i,k,j)
1196 if (l==12) emit_sum(l) = emit_sum(l) + e12(i,k,j)
1197 if (l==13) emit_sum(l) = emit_sum(l) + e13(i,k,j)
1198 if (l==14) emit_sum(l) = emit_sum(l) + e14(i,k,j)
1199 if (l==15) emit_sum(l) = emit_sum(l) + e15(i,k,j)
1200 if (l==16) emit_sum(l) = emit_sum(l) + e16(i,k,j)
1201 if (l==17) emit_sum(l) = emit_sum(l) + e17(i,k,j)
1202 if (l==18) emit_sum(l) = emit_sum(l) + e18(i,k,j)
1203 if (l==19) emit_sum(l) = emit_sum(l) + e19(i,k,j)
1204 if (l==20) emit_sum(l) = emit_sum(l) + e20(i,k,j)
1205
1206 if (l==21) emit_sum(l) = emit_sum(l) + e21(i,k,j)
1207 end do
1208 end do
1209 end do
1210 2900 continue
1211
1212 ! output the chem_sum and emit_sum
1213 print 9110, fromwhere, its, ite, jts, jte
1214 print 9100, 'chem_sum'
1215 fact = 1.0/(dtstep*efact1)
1216 print 9120, (l, fact*chem_sum(l), l=1,num_chem)
1217 print 9100, 'emit_sum'
1218 print 9120, (l, emit_sum(l), l=1,min(nemit,nemit_maxd))
1219
1220 9100 format( a )
1221 9110 format( / 'addemiss_masscheck output, fromwhere = ', a / &
1222 'its, ite, jts, jte =', 4i5 )
1223 9120 format( 5( i5, 1pe11.3 ) )
1224
1225
1226 return
1227 END subroutine addemiss_masscheck
1228