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