| blob | 0681ac36f6d7972c0a1d6cdcad178de8ad269ed1 |
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 references and terms of use
8 !**********************************************************************************
9 MODULE module_mosaic_addemiss
10 !WRF:MODEL_LAYER:CHEMICS
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
19 integer, parameter :: mosaic_addemiss_active = 1
20 ! only do emissions when this is positive
21 ! (when it is negative, emissions tendencies are zero)
23 integer, parameter :: mosaic_addemiss_masscheck = -1
24 ! only do emissions masscheck calcs when this is positive
26 CONTAINS
30 !----------------------------------------------------------------------
31 subroutine mosaic_addemiss( id, dtstep, u10, v10, alt, dz8w, xland, &
32 config_flags, chem, slai, ust, smois, ivgtyp, isltyp, &
33 emis_ant,ebu,biom_active,dust_opt, &
34 !czhao
35 ktau, u_phy,v_phy,rho_phy,g,dx,erod, &
36 dust_emiss_active, seasalt_emiss_active, &
37 seas_flux,emis_dust, &
38 ids,ide, jds,jde, kds,kde, &
39 ims,ime, jms,jme, kms,kme, &
40 its,ite, jts,jte, kts,kte )
41 !
42 ! adds emissions for mosaic aerosol species
43 ! (i.e., emissions tendencies over time dtstep are applied
44 ! to the aerosol concentrations)
45 !
47 USE module_configure, only: grid_config_rec_type
48 USE module_state_description
49 ! USE module_state_description, only: num_chem, param_first_scalar, &
50 ! emiss_inpt_default, emiss_inpt_pnnl_rs, emiss_inpt_pnnl_cm,num_emis_ant
51 USE module_data_mosaic_asect
53 IMPLICIT NONE
55 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
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
62 INTEGER, INTENT(IN) :: dust_emiss_active, seasalt_emiss_active, biom_active, dust_opt
64 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: seas_flux
66 REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL, &
67 INTENT(INOUT ) :: &
68 emis_dust
70 !czhao
71 INTEGER, INTENT(IN) :: ktau
72 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
73 INTENT(IN) :: u_phy,v_phy,rho_phy
74 REAL, INTENT(IN ) :: dx, g
75 REAL, DIMENSION( ims:ime, jms:jme,3),&
76 INTENT(IN ) :: erod
79 REAL, INTENT(IN ) :: dtstep
81 ! 10-m wind speed components (m/s)
82 REAL, DIMENSION( ims:ime , jms:jme ) , &
83 INTENT(IN ) :: u10, v10, xland, slai, ust
84 INTEGER, DIMENSION( ims:ime , jms:jme ) , &
85 INTENT(IN ) :: ivgtyp, isltyp
87 ! trace species mixing ratios (aerosol mass = ug/kg-air; number = #/kg-air)
88 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
89 INTENT(INOUT ) :: chem
90 !
91 ! aerosol emissions arrays ((ug/m3)*m/s)
92 !
93 ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
94 REAL, DIMENSION( ims:ime, 1:config_flags%kemit, jms:jme,num_emis_ant ), &
95 INTENT(IN ) :: &
96 emis_ant
98 REAL, DIMENSION( ims:ime, kms:kme, jms:jme,num_ebu ), &
99 INTENT(IN ) :: &
100 ebu
102 ! 1/(dry air density) and layer thickness (m)
103 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
104 INTENT(IN ) :: alt, dz8w
106 REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) , &
107 INTENT(INOUT) :: smois
109 ! local variables
110 integer i, j, k, l, n
111 integer iphase, itype
112 integer p1st
114 real, parameter :: efact1 = 1.0
115 real :: aem_so4, aem_no3, aem_cl, aem_msa, aem_co3, aem_nh4, &
116 aem_na, aem_ca, aem_oin, aem_oc, aem_bc, aem_num
117 real dum, fact
120 ! fraction of sorgam i/aitken mode emissions that go to each
121 ! of the mosaic 8 "standard" sections
122 !now the anthropogenic and biomass burning emissiosn use same distribution
123 !it can be improve later --czhao
124 ! not 100%, because most mass falling down to the size less than the lowest boundary
125 real, save :: fr8b_aem_sorgam_i(8) = &
126 (/ 0.965, 0.035, 0.000, 0.000, &
127 0.000, 0.000, 0.000, 0.000 /)
129 ! fraction of sorgam j/accum mode emissions that go to each
130 ! of the mosaic 8 "standard" sections
131 real, save :: fr8b_aem_sorgam_j(8) = &
132 (/ 0.026, 0.147, 0.350, 0.332, &
133 0.125, 0.019, 0.001, 0.000/)
135 ! fraction of sorgam coarse mode emissions that go to each
136 ! of the mosaic 8 "standard" sections
137 ! not 100%, because most mass falling up to the size larger than the highest boundary
138 real, save :: fr8b_aem_sorgam_c(8) = &
139 (/ 0.000, 0.000, 0.000, 0.002, &
140 0.021, 0.110, 0.275, 0.592 /)
142 ! fraction of mosaic fine (< 2.5 um) emissions that go to each
143 ! of the mosaic 8 "standard" sections
144 !wig 1-Apr-2005, Updated fractional breakdown between bins. (See also
145 ! bdy_chem_value_mosaic and mosaic_init_wrf_mixrats_opt2
146 ! in module_mosaic_initmixrats.F.) Note that the values
147 ! here no longer match the other two subroutines.
148 !rce 10-may-2005, changed fr8b_aem_mosaic_f & fr8b_aem_mosaic_c
149 ! to values determined by jdf
150 real, save :: fr8b_aem_mosaic_f(8) = &
151 (/ 0.060, 0.045, 0.245, 0.400, 0.100, 0.150, 0., 0./) !10-may-2005
152 ! (/ 0.100, 0.045, 0.230, 0.375, 0.100, 0.150, 0., 0./) !1-Apr-2005 values
153 ! (/ 0.0275, 0.0426, 0.2303, 0.3885, 0.1100, 0.2011, 0., 0./) !15-Nov-2004 values
154 ! (/ 0.01, 0.05, 0.145, 0.60, 0.145, 0.05, 0.00, 0.00 /)
155 ! (/ 0.04, 0.10, 0.35, 0.29, 0.15, 0.07, 0.0, 0.0 /)
157 ! fraction of mosaic coarse (> 2.5 um) emissions that go to each
158 ! of the mosaic 8 "standard" sections
159 real, save :: fr8b_aem_mosaic_c(8) = &
160 (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.300, 0.700 /) ! 10-may-2005
161 ! (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.933, 0.067 /) ! as of apr-2005
162 ! (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.16, 0.84 /) ! "old"
164 !++ SAN Fractionation for PREP-Chem Fire emissions ++!
165 ! New fractions for biomass burning emissions. Based on Janhall et. al. (2010)
166 ! Based on average fresh smoke emissions from 20 data points
167 ! mean_Dp = .117um, geometric_std = 1.7
168 real, save :: fr8b_bburn_mosaic(8) = &
169 (/ 0.0092, 0.1385, 0.4548, 0.3388, 0.0567, 0.002, 0.0, 0.0/)
170 !--SAN
173 ! fraction of TNO black carbon <1um emissions that go into each of
174 ! the mosaic 8 "standard" sections - Doug (4/3/2011)
175 real, save :: fr8b_tno_bc1(8) = &
176 (/ 0.0494, 0.3795, 0.4714, 0.0967, 0.003, 0.0, 0.0, 0.0 /)
178 ! fraction of TNO elemental carbon 1um-2.5um emissions that go into each of
179 ! the mosaic 8 "standard" sections - Doug (4/3/2011)
180 real, save :: fr8b_tno_ec25(8) = &
181 (/ 0.0, 0.0, 0.0, 0.0, 0.40, 0.60, 0.0, 0.0 /)
183 ! fraction of TNO organic carbon <2.5um domestic combustion emissions that go into each of
184 ! the mosaic 8 "standard" sections - Doug (4/3/2011)
185 real, save :: fr8b_tno_oc_dom(8) = &
186 (/ 0.0358, 0.1325, 0.2704, 0.3502, 0.1904, 0.0657, 0.0, 0.0 /)
188 ! fraction of TNO organic carbon <2.5um traffic (and other) emissions that go into each of
189 ! the mosaic 8 "standard" sections - Doug (4/3/2011)
190 real, save :: fr8b_tno_oc_tra(8) = &
191 (/ 0.0063, 0.0877, 0.3496, 0.4054, 0.1376, 0.0134, 0.0, 0.0 /)
193 ! fraction of TNO "OIN" [PM2.5 minus sum(carbon<2.5)] emissions that go into each of
194 ! the mosaic 8 "standard" sections - Doug (4/3/2011) --- based on mosaic fine mode
195 real, save :: fr8b_tno_mosaic_f(8) = &
196 (/ 0.060, 0.045, 0.245, 0.400, 0.100, 0.150, 0.0, 0.0/)
198 ! fraction of TNO 2.5-10um emissions that go into each of
199 ! the mosaic 8 "standard" sections - Doug (4/3/2011) --- based on mosaic coarse mode
200 real, save :: fr8b_tno_mosaic_c(8) = &
201 (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.300, 0.700 /)
205 ! following 5 arrays correspond to the above "fr8b_" arrays
206 ! but are set at run time base on input (namelist) parameters:
207 ! only the sorgam or mosaic arrays are non-zero, depending on
208 ! emiss_inpt_opt
209 ! when nsize_aer=4 (=number of sections), the values are
210 ! calculated by adding two of the 8-section values
211 real :: fr_aem_sorgam_i(8)
212 real :: fr_aem_sorgam_j(8)
213 real :: fr_aem_sorgam_c(8)
214 real :: fr_aem_mosaic_f(8)
215 real :: fr_aem_mosaic_c(8)
216 real :: fr_tno_bc1(8)
217 real :: fr_tno_ec25(8)
218 real :: fr_tno_oc_dom(8)
219 real :: fr_tno_oc_tra(8)
220 real :: fr_tno_mosaic_f(8)
221 real :: fr_tno_mosaic_c(8)
222 !++ SAN
223 real :: fr_aem_gc2mosaic_f(8) ! extra arrays for GOCART->MOSAIC mapping
224 real :: fr_aem_gc2mosaic_c(8)
225 real :: bburn_mosaic_f(8) ! Arrays to hold biomass-burning size
226 real :: bburn_mosaic_c(8) ! Arrays to hold biomass-burning size
227 !-- SAN
228 double precision :: chem_sum(num_chem)
230 character(len=80) :: msg
233 ! *** currently only works with ntype_aer = 1
234 itype = 1
235 iphase = ai_phase
238 ! if (num_ebu.le.0 ) then
239 ! call wrf_error_fatal( 'mosaic_addemiss: no biomass burning species' )
240 ! print*,'no biomass burning species',num_ebu
241 ! stop
242 ! endif
244 ! compute factors used for apportioning either
245 ! the MADE-SORGAM emissions (i=aitken, j=accum, c=coarse modes) OR
246 ! the MOSAIC emission (f=fine (< 2.5 um), c=coarse (> 2.5 um))
247 ! to each size section
248 !
249 ! note: the fr8b_aer_xxxxxx_y values are specific to the mosaic 8 bin
250 ! structure with dlo_sect(1)=0.039 um and dhi_sect(8)=10.0 um,
251 ! also, the fr8b_aem_sorgam_y are specific for the assumed
252 ! dgvem_i/j/c used in the module_data_sorgam.F code
253 ! also, the fr8b_aem_mosaic_y values are specific for the assumed (by us)
254 ! size distribution for fine and coarse primary emissions
255 !
256 ! when there are 4 bins (nsize_aer=4), each of these "wider" bins
257 ! corresponds to 2 of the "narrower" bins of the 8 bin structure
258 !
259 ! note: if fr_aem_sorgam_y > 0, then fr_aem_mosaic_y = 0, and vice-versa
260 !
261 if ((nsize_aer(itype) .ne. 4) .and. (nsize_aer(itype) .ne. 8)) then
262 write(msg,'(a,i5)') &
263 'subr mosaic_addemiss - nsize_aer(itype) must be ' // &
264 '4 or 8 but = ', nsize_aer(itype)
265 call wrf_error_fatal( msg )
266 end if
268 fr_aem_sorgam_i(:) = 0.0
269 fr_aem_sorgam_j(:) = 0.0
270 fr_aem_sorgam_c(:) = 0.0
271 fr_aem_mosaic_f(:) = 0.0
272 fr_aem_mosaic_c(:) = 0.0
273 fr_tno_bc1(:) = 0.0
274 fr_tno_ec25(:) = 0.0
275 fr_tno_oc_dom(:) = 0.0
276 fr_tno_oc_tra(:) = 0.0
277 fr_tno_mosaic_f(:) = 0.0
278 fr_tno_mosaic_c(:) = 0.0
279 fr_aem_gc2mosaic_f(:) = 0.0
280 fr_aem_gc2mosaic_c(:) = 0.0
282 emiss_inpt_select_1: SELECT CASE( config_flags%emiss_inpt_opt )
284 CASE( emiss_inpt_default, emiss_inpt_pnnl_rs )
285 if (nsize_aer(itype) .eq. 8) then
286 fr_aem_sorgam_i(:) = fr8b_aem_sorgam_i(:)
287 fr_aem_sorgam_j(:) = fr8b_aem_sorgam_j(:)
288 fr_aem_sorgam_c(:) = fr8b_aem_sorgam_c(:)
289 else if (nsize_aer(itype) .eq. 4) then
290 do n = 1, nsize_aer(itype)
291 fr_aem_sorgam_i(n) = fr8b_aem_sorgam_i(2*n-1) &
292 + fr8b_aem_sorgam_i(2*n)
293 fr_aem_sorgam_j(n) = fr8b_aem_sorgam_j(2*n-1) &
294 + fr8b_aem_sorgam_j(2*n)
295 fr_aem_sorgam_c(n) = fr8b_aem_sorgam_c(2*n-1) &
296 + fr8b_aem_sorgam_c(2*n)
297 end do
298 end if
300 !++ SAN, 2015-04-09
301 ! Added mapping routine from GOCART aerosol emission variables to MOSAIC
302 ! This will be the case if emissions were prepared using PREP-Chem in RAMD2-GOCART mode.
303 ! (Use with emiss_opt = 5/6, ECPTEC/GOCART_ECPTEC, emiss_inpt_opt == emiss_inpt_default)
304 !!! Maybe more consistant to add a new emiss_inpt_opt, but this requires more extensive changes
305 !!! and testing in other parts of WRF-Chem, may cause confusion... Thoughts?
307 if( config_flags%emiss_opt == 5 .or. config_flags%emiss_opt == 6 ) then
308 CALL wrf_debug(15,'mosaic_addemiss: emiss_opt = eptec')
309 CALL wrf_debug(15,'mosaic_addemiss: gocart speciation being mapped to mosaic')
311 fr_aem_sorgam_i(:) = 0.0
312 fr_aem_sorgam_j(:) = 0.0
313 fr_aem_sorgam_c(:) = 0.0
315 ! Use FM MOSAIC mapping for SO4, OC, BC and PM2.5, CM mosaic for PM10
316 if (nsize_aer(itype) .eq. 8) then
317 fr_aem_gc2mosaic_f(:) = fr8b_aem_mosaic_f(:)
318 fr_aem_gc2mosaic_c(:) = fr8b_aem_mosaic_c(:)
319 elseif (nsize_aer(itype) .eq. 4) then
320 do n = 1, nsize_aer(itype)
321 fr_aem_gc2mosaic_f(n) = fr8b_aem_mosaic_f(2*n-1) + fr8b_aem_mosaic_f(2*n)
322 fr_aem_gc2mosaic_c(n) = fr8b_aem_mosaic_c(2*n-1) + fr8b_aem_mosaic_c(2*n)
323 end do
324 endif
325 endif
326 !-- SAN
327 CASE( emiss_inpt_pnnl_cm )
328 if (nsize_aer(itype) .eq. 8) then
329 fr_aem_mosaic_f(:) = fr8b_aem_mosaic_f(:)
330 fr_aem_mosaic_c(:) = fr8b_aem_mosaic_c(:)
331 else if (nsize_aer(itype) .eq. 4) then
332 do n = 1, nsize_aer(itype)
333 fr_aem_mosaic_f(n) = fr8b_aem_mosaic_f(2*n-1) &
334 + fr8b_aem_mosaic_f(2*n)
335 fr_aem_mosaic_c(n) = fr8b_aem_mosaic_c(2*n-1) &
336 + fr8b_aem_mosaic_c(2*n)
337 end do
338 end if
340 CASE( emiss_inpt_tno ) ! Doug - added 4/3/2011
341 if (nsize_aer(itype) .eq. 8) then
342 fr_tno_bc1(:) = fr8b_tno_bc1(:)
343 fr_tno_ec25(:) = fr8b_tno_ec25(:)
344 fr_tno_oc_dom(:) = fr8b_tno_oc_dom(:)
345 fr_tno_oc_tra(:) = fr8b_tno_oc_tra(:)
346 fr_tno_mosaic_c(:) = fr8b_tno_mosaic_c(:)
347 fr_tno_mosaic_f(:) = fr8b_tno_mosaic_f(:)
348 else if (nsize_aer(itype) .eq. 4) then
349 do n = 1, nsize_aer(itype)
350 fr_tno_bc1(n) = fr8b_tno_bc1(2*n-1) &
351 + fr8b_tno_bc1(2*n)
352 fr_tno_ec25(n) = fr8b_tno_ec25(2*n-1) &
353 + fr8b_tno_ec25(2*n)
354 fr_tno_oc_dom(n) = fr8b_tno_oc_dom(2*n-1) &
355 + fr8b_tno_oc_dom(2*n)
356 fr_tno_oc_tra(n) = fr8b_tno_oc_tra(2*n-1) &
357 + fr8b_tno_oc_tra(2*n)
358 fr_tno_mosaic_c(n) = fr8b_tno_mosaic_c(2*n-1) &
359 + fr8b_tno_mosaic_c(2*n)
360 fr_tno_mosaic_f(n) = fr8b_tno_mosaic_f(2*n-1) &
361 + fr8b_tno_mosaic_f(2*n)
362 end do
363 end if
365 CASE DEFAULT
366 return
368 END SELECT emiss_inpt_select_1
370 !++ SAN 2015-04-09: Mapping fire emissions to MOSAIC
371 ! get arrays for apportioning mass into MOSAIC bins for fire emissions
372 fire_inpt_select: SELECT CASE (config_flags%biomass_burn_opt)
373 CASE (BIOMASSB,BIOMASSB_MOZC)
374 if (nsize_aer(itype) .eq. 8) then
375 bburn_mosaic_f(:) = fr8b_bburn_mosaic(:)
376 !++SAN
377 bburn_mosaic_c(:) = fr8b_aem_mosaic_c(:)
378 !--SAN
379 else if (nsize_aer(itype) .eq. 4) then
380 do n = 1, nsize_aer(itype)
381 bburn_mosaic_f(n) = fr8b_bburn_mosaic(2*n-1) + fr8b_bburn_mosaic(2*n)
382 !++SAN
383 bburn_mosaic_c(n) = fr8b_aem_mosaic_c(2*n-1) + fr8b_aem_mosaic_c(2*n)
384 !--SAN
385 end do
386 end if
387 END SELECT fire_inpt_select
389 ! when mosaic_addemiss_active <= 0, set fr's to zero,
390 ! which causes the changes to chem(...) to be zero
391 if (mosaic_addemiss_active <= 0 .and. biom_active <= 0) then
392 fr_aem_sorgam_i(:) = 0.0
393 fr_aem_sorgam_j(:) = 0.0
394 fr_aem_sorgam_c(:) = 0.0
395 fr_aem_mosaic_f(:) = 0.0
396 fr_aem_mosaic_c(:) = 0.0
397 fr_tno_bc1(:) = 0.0
398 fr_tno_ec25(:) = 0.0
399 fr_tno_oc_dom(:) = 0.0
400 fr_tno_oc_tra(:) = 0.0
401 fr_tno_mosaic_f(:) = 0.0
402 fr_tno_mosaic_c(:) = 0.0
403 end if
406 ! do mass check initial calc
407 if (mosaic_addemiss_masscheck > 0) call addemiss_masscheck( &
408 id, config_flags, 1, 'mosaic_ademiss', &
409 dtstep, efact1, dz8w, chem, chem_sum, &
410 ids,ide, jds,jde, kds,kde, &
411 ims,ime, jms,jme, kms,kme, &
412 its,ite, jts,jte, kts,kte, &
413 14, &
414 emis_ant(ims, 1,jms,p_e_pm_10),emis_ant(ims, 1,jms,p_e_pm_25), &
415 emis_ant(ims, 1,jms,p_e_pm25i),emis_ant(ims, 1,jms,p_e_pm25j), &
416 emis_ant(ims, 1,jms,p_e_eci),emis_ant(ims, 1,jms,p_e_ecj), &
417 emis_ant(ims, 1,jms,p_e_orgi),emis_ant(ims, 1,jms,p_e_orgj), &
418 emis_ant(ims, 1,jms,p_e_so4j),emis_ant(ims, 1,jms,p_e_so4c), &
419 emis_ant(ims, 1,jms,p_e_no3j),emis_ant(ims, 1,jms,p_e_no3c), &
420 emis_ant(ims, 1,jms,p_e_orgc),emis_ant(ims, 1,jms,p_e_ecc), &
421 emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), &
422 emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), &
423 emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), &
424 emis_ant(ims, 1,jms,p_e_ecc))
427 p1st = param_first_scalar
429 !
430 ! apply emissions at each section and grid point
431 !
432 do 1900 n = 1, nsize_aer(itype)
434 do 1830 j = jts, jte
435 do 1820 k = 1, min(config_flags%kemit,kte)
436 do 1810 i = its, ite
438 ! compute anthropogenic mass emissions [(ug/m3)*m/s] for each species
439 ! using the apportioning fractions
440 aem_so4 = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_so4j) &
441 + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_so4c) &
442 + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_so4i) &
443 + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_so4j)
445 aem_no3 = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_no3j) &
446 + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_no3c) &
447 + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_no3i) &
448 + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_no3j)
450 aem_oc = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_orgj) &
451 + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_orgc) &
452 + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_orgi) &
453 + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_orgj) &
454 + fr_tno_oc_dom(n)*emis_ant(i,k,j,p_e_oc_dom) &
455 + fr_tno_oc_tra(n)*emis_ant(i,k,j,p_e_oc_tra) &
456 + fr_tno_mosaic_c(n)*emis_ant(i,k,j,p_e_oc_25_10) &
457 !++SAN
458 + fr_aem_gc2mosaic_f(n)*emis_ant(i,k,j,p_e_oc)
459 !--SAN
461 chem_select_1 : SELECT CASE( config_flags%chem_opt )
462 CASE(SAPRC99_MOSAIC_4BIN_VBS2_KPP, SAPRC99_MOSAIC_8BIN_VBS2_AQ_KPP, &
463 SAPRC99_MOSAIC_8BIN_VBS2_KPP)!BSINGH(12/04/2013): Added SAPRC 8 bin and non-aq on 04/03/2014! Set the oc to zero for VBS
464 aem_oc = 0.0
465 END SELECT chem_select_1
467 aem_bc = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_ecj) &
468 + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_ecc) &
469 + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_eci) &
470 + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_ecj) &
471 + fr_tno_bc1(n)*emis_ant(i,k,j,p_e_bc_1) &
472 + fr_tno_ec25(n)*emis_ant(i,k,j,p_e_ec_1_25) &
473 + fr_tno_mosaic_c(n)*emis_ant(i,k,j,p_e_ec_25_10) &
474 !++SAN
475 + fr_aem_gc2mosaic_f(n)*emis_ant(i,k,j,p_e_bc)
476 !--SAN
478 aem_oin = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_pm25j) &
479 + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_pm_10) &
480 + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_pm25i) &
481 + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_pm25j) &
482 !++SAN
483 + fr_aem_sorgam_c(n)*emis_ant(i,k,j,p_e_pm_10) &
484 !--SAN
485 + fr_tno_mosaic_f(n)*emis_ant(i,k,j,p_e_oin_25) &
486 + fr_tno_mosaic_c(n)*emis_ant(i,k,j,p_e_oin_10) &
487 !++SAN
488 + fr_aem_gc2mosaic_f(n)*emis_ant(i,k,j,p_e_pm25) &
489 + fr_aem_gc2mosaic_f(n)*emis_ant(i,k,j,p_e_pm10)
490 !--SAN
493 ! emissions for these species are currently zero
494 aem_nh4 = 0.0
495 aem_na = 0.0
496 aem_cl = 0.0
497 aem_ca = 0.0
498 aem_co3 = 0.0
499 aem_msa = 0.0
501 chem_select_2 : SELECT CASE( config_flags%chem_opt )
502 CASE(MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP)
503 aem_nh4 = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_nh4j) &
504 + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_nh4c) &
505 + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_nh4i) &
506 + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_nh4j)
508 aem_na = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_naj) &
509 + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_nac) &
510 + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_nai) &
511 + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_naj)
513 aem_cl = fr_aem_mosaic_f(n)*emis_ant(i,k,j,p_e_clj) &
514 + fr_aem_mosaic_c(n)*emis_ant(i,k,j,p_e_clc) &
515 + fr_aem_sorgam_i(n)*emis_ant(i,k,j,p_e_cli) &
516 + fr_aem_sorgam_j(n)*emis_ant(i,k,j,p_e_clj)
517 END SELECT chem_select_2
519 ! compute number emissions
520 ! first sum the mass-emissions/density
521 aem_num = &
522 (aem_so4/dens_so4_aer) + (aem_no3/dens_no3_aer) + &
523 (aem_cl /dens_cl_aer ) + (aem_msa/dens_msa_aer) + &
524 (aem_co3/dens_co3_aer) + (aem_nh4/dens_nh4_aer) + &
525 (aem_na /dens_na_aer ) + (aem_ca /dens_ca_aer ) + &
526 (aem_oin/dens_oin_aer) + (aem_oc /dens_oc_aer ) + &
527 (aem_bc /dens_bc_aer )
529 ! then multiply by 1.0e-6 to convert ug to g
530 ! and divide by particle volume at center of section (cm3)
531 aem_num = aem_num*1.0e-6/volumcen_sect(n,itype)
533 ! apply the emissions and convert from flux to mixing ratio
534 ! fact = (dtstep/dz8w(i,k,j))*(28.966/1000.)
535 fact = (dtstep/dz8w(i,k,j))*alt(i,k,j)
537 ! rce 22-nov-2004 - change to using the "..._aer" species pointers
538 l = lptr_so4_aer(n,itype,iphase)
539 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_so4*fact
541 l = lptr_no3_aer(n,itype,iphase)
542 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_no3*fact
544 l = lptr_cl_aer(n,itype,iphase)
545 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_cl*fact
547 l = lptr_msa_aer(n,itype,iphase)
548 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_msa*fact
550 l = lptr_co3_aer(n,itype,iphase)
551 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_co3*fact
553 l = lptr_nh4_aer(n,itype,iphase)
554 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_nh4*fact
556 l = lptr_na_aer(n,itype,iphase)
557 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_na*fact
559 l = lptr_ca_aer(n,itype,iphase)
560 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_ca*fact
562 l = lptr_oin_aer(n,itype,iphase)
563 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oin*fact
565 l = lptr_oc_aer(n,itype,iphase)
566 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oc*fact
568 l = lptr_bc_aer(n,itype,iphase)
569 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_bc*fact
571 l = numptr_aer(n,itype,iphase)
572 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_num*fact
574 1810 continue
575 1820 continue
576 1830 continue
578 1900 continue
580 if (num_ebu <= 0 ) then
581 return
582 endif
584 !++ SAN, 2015-04-09.
585 ! ---------------------- ADD BIOMASS BURNING EMISSIONS ----------------------------!
586 ! - new case select for adding biomass burning emissions.
587 ! First reset all emissions to make sure we don't double count anthropogenic emissions
588 aem_so4 = 0.
589 aem_no3 = 0.
590 aem_oc = 0.
591 aem_bc = 0.
592 aem_oin = 0.
593 aem_num = 0.
595 fire_select: SELECT CASE(config_flags%biomass_burn_opt)
596 CASE (BIOMASSB,BIOMASSB_MOZC)
597 CALL wrf_debug(15,'mosaic_addemiss: adding fire emissions to MOSAIC')
598 ! CALL wrf_debug(15,'e_oin_fm = PM2.5 - OC - BC')
599 ! CALL wrf_debug(15,'e_oin_cm = PM10 - PM2.5')
601 size_loop: &
602 do n = 1, nsize_aer(itype)
603 do j = jts, jte
604 do k = kts, kte ! Loop up to kte for fire emissions (with plumerise)
605 do i = its, ite
606 ! compute mass biomass burning emissions [(ug/m3)*m/s] for each species
607 ! using the apportioning fractions
608 IF ( p_ebu_sulf .gt. 1) aem_so4 = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_sulf)
609 aem_oc = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_oc)
610 aem_bc = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_bc)
611 ! Option to calculate OIN fraction of total PM for fire emissions
612 ! Assume all OC and BC is in FM.
613 aem_oin = bburn_mosaic_f(n)*ebu(i,k,j,p_ebu_pm25) &
614 + bburn_mosaic_c(n)*ebu(i,k,j,p_ebu_pm10)
616 chem_select_3 : SELECT CASE( config_flags%chem_opt )
617 CASE(SAPRC99_MOSAIC_4BIN_VBS2_KPP) ! Set the oc to zero for VBS
618 aem_oc = 0.0
619 END SELECT chem_select_3
621 ! emissions for these species are currently zero
622 aem_nh4 = 0.0
623 aem_na = 0.0
624 aem_cl = 0.0
625 aem_ca = 0.0
626 aem_co3 = 0.0
627 aem_msa = 0.0
629 ! compute number emissions
630 ! first sum the mass-emissions/density
631 aem_num = &
632 (aem_so4/dens_so4_aer) + (aem_no3/dens_no3_aer) + &
633 (aem_cl /dens_cl_aer ) + (aem_msa/dens_msa_aer) + &
634 (aem_co3/dens_co3_aer) + (aem_nh4/dens_nh4_aer) + &
635 (aem_na /dens_na_aer ) + (aem_ca /dens_ca_aer ) + &
636 (aem_oin/dens_oin_aer) + (aem_oc /dens_oc_aer ) + &
637 (aem_bc /dens_bc_aer )
639 ! then multiply by 1.0e-6 to convert ug to g
640 ! and divide by particle volume at center of section (cm3)
641 aem_num = aem_num*1.0e-6/volumcen_sect(n,itype)
643 ! apply the emissions and convert from flux to mixing ratio
644 ! fact = (dtstep/dz8w(i,k,j))*(28.966/1000.)
645 fact = (dtstep/dz8w(i,k,j))*alt(i,k,j)
647 ! rce 22-nov-2004 - change to using the "..._aer" species pointers
648 l = lptr_so4_aer(n,itype,iphase)
649 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_so4*fact
651 l = lptr_no3_aer(n,itype,iphase)
652 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_no3*fact
654 l = lptr_cl_aer(n,itype,iphase)
655 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_cl*fact
657 l = lptr_msa_aer(n,itype,iphase)
658 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_msa*fact
660 l = lptr_co3_aer(n,itype,iphase)
661 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_co3*fact
663 l = lptr_nh4_aer(n,itype,iphase)
664 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_nh4*fact
666 l = lptr_na_aer(n,itype,iphase)
667 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_na*fact
669 l = lptr_ca_aer(n,itype,iphase)
670 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_ca*fact
672 l = lptr_oin_aer(n,itype,iphase)
673 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oin*fact
675 l = lptr_oc_aer(n,itype,iphase)
676 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_oc*fact
678 l = lptr_bc_aer(n,itype,iphase)
679 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_bc*fact
681 l = numptr_aer(n,itype,iphase)
682 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + aem_num*fact
684 end do
685 end do
686 end do
687 end do size_loop
688 END SELECT fire_select
689 !-- SAN - end adding fire emissions.
692 ! do mass check final calc
693 if (mosaic_addemiss_masscheck > 0) call addemiss_masscheck( &
694 id, config_flags, 2, 'mosaic_ademiss', &
695 dtstep, efact1, dz8w, chem, chem_sum, &
696 ids,ide, jds,jde, kds,kde, &
697 ims,ime, jms,jme, kms,kme, &
698 its,ite, jts,jte, kts,kte, &
699 14, &
700 emis_ant(ims, 1,jms,p_e_pm_10),emis_ant(ims, 1,jms,p_e_pm_25), &
701 emis_ant(ims, 1,jms,p_e_pm25i),emis_ant(ims, 1,jms,p_e_pm25j), &
702 emis_ant(ims, 1,jms,p_e_eci),emis_ant(ims, 1,jms,p_e_ecj), &
703 emis_ant(ims, 1,jms,p_e_orgi),emis_ant(ims, 1,jms,p_e_orgj), &
704 emis_ant(ims, 1,jms,p_e_so4j),emis_ant(ims, 1,jms,p_e_so4c), &
705 emis_ant(ims, 1,jms,p_e_no3j),emis_ant(ims, 1,jms,p_e_no3c), &
706 emis_ant(ims, 1,jms,p_e_orgc),emis_ant(ims, 1,jms,p_e_ecc), &
707 emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), &
708 emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), &
709 emis_ant(ims, 1,jms,p_e_ecc),emis_ant(ims, 1,jms,p_e_ecc), &
710 emis_ant(ims, 1,jms,p_e_ecc))
713 ! do seasalt emissions
714 if (seasalt_emiss_active > 0) &
715 call mosaic_seasalt_emiss( &
716 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
717 seas_flux, &
718 ids,ide, jds,jde, kds,kde, &
719 ims,ime, jms,jme, kms,kme, &
720 its,ite, jts,jte, kts,kte, seasalt_emiss_active )
722 ! if (seasalt_emiss_active == 2) &
723 ! call mosaic_seasalt_emiss( &
724 ! id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
725 ! ids,ide, jds,jde, kds,kde, &
726 ! ims,ime, jms,jme, kms,kme, &
727 ! its,ite, jts,jte, kts,kte )
728 !
729 ! do dust emissions
730 ! jdf: preliminary version that has not been made generic for situation
731 if (dust_opt == 2) then
732 !czhao+++++++++++++++++++++++++++
733 call wrf_message("WARNING: You are calling DUSTRAN dust emission scheme with MOSAIC, which is highly experimental and not recommended for use. Please use dust_opt==13")
734 !czhao---------------------------
735 call mosaic_dust_emiss( slai, ust, smois, ivgtyp, isltyp, &
736 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
737 ids,ide, jds,jde, kds,kde, &
738 ims,ime, jms,jme, kms,kme, &
739 its,ite, jts,jte, kts,kte )
740 endif
741 !czhao+++++++++++++++++++++++++++
742 if (dust_opt == 13) then
743 !czhao---------------------------
744 call mosaic_dust_gocartemis (dtstep,config_flags%start_month, &
745 config_flags%num_soil_layers,alt,u_phy, &
746 v_phy,chem,rho_phy,dz8w,smois,u10,v10,erod, &
747 isltyp,xland,g, &
748 emis_dust, &
749 ids,ide, jds,jde, kds,kde, &
750 ims,ime, jms,jme, kms,kme, &
751 its,ite, jts,jte, kts,kte)
752 endif
755 return
758 END subroutine mosaic_addemiss
762 !----------------------------------------------------------------------
763 subroutine mosaic_seasalt_emiss( &
764 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
765 seas_flux, &
766 ids,ide, jds,jde, kds,kde, &
767 ims,ime, jms,jme, kms,kme, &
768 its,ite, jts,jte, kts,kte, seasalt_emiss_active )
769 !
770 ! adds seasalt emissions for mosaic aerosol species
771 ! (i.e., seasalt emissions tendencies over time dtstep are applied
772 ! to the aerosol mixing ratios)
773 !
775 USE module_configure, only: grid_config_rec_type
776 USE module_state_description, only: num_chem, param_first_scalar
777 USE module_data_mosaic_asect
779 IMPLICIT NONE
781 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
783 INTEGER, INTENT(IN ) :: id, &
784 ids,ide, jds,jde, kds,kde, &
785 ims,ime, jms,jme, kms,kme, &
786 its,ite, jts,jte, kts,kte
788 REAL, INTENT(IN ) :: dtstep
789 INTEGER, INTENT(IN) :: seasalt_emiss_active
791 ! 10-m wind speed components (m/s)
792 REAL, DIMENSION( ims:ime , jms:jme ), &
793 INTENT(IN ) :: u10, v10, xland
795 ! trace species mixing ratios (aerosol mass = ug/kg; number = #/kg)
796 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
797 INTENT(INOUT ) :: chem
799 ! alt = 1.0/(dry air density) in (m3/kg)
800 ! dz8w = layer thickness in (m)
801 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
802 INTENT(IN ) :: alt, dz8w
804 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: seas_flux
806 ! local variables
807 integer i, j, k, l, l_na, l_cl, n, l_oc
808 integer iphase, itype
809 integer p1st
811 real dum, dumdlo, dumdhi, dumoceanfrac, dumspd10
812 real factaa, factbb, fracna, fraccl, fracorg
814 real :: ssemfact_numb( maxd_asize, maxd_atype )
815 real :: ssemfact_mass( maxd_asize, maxd_atype )
818 ! Factors for Fuentes et al (ACP, 2011, doi:10.5194/acp-11-2585-2011)
819 ! seasalt emission scheme.
820 ! These are for calculating the seawater_OC content dependent factors.
821 real, save :: alpha_f1(4) = &
822 (/ 12.328, 38.077, 102.31, 281.65 /)
823 real, save :: alpha_f2(4) = &
824 (/ 2.2958, 8.0935, 25.251, 46.80 /)
825 real, save :: alpha_f3(4) = &
826 (/ 0.00452, 0.00705, 0.00080, 0.000761 /)
827 real, save :: beta_f1(4) = &
828 (/ 0.0311, -0.031633, 0.013154, -0.0017762 /)
829 real, save :: beta_f2(4) = &
830 (/ -13.916, 35.73, -9.7651, 1.1665 /)
831 real, save :: beta_f3(4) = &
832 (/ 4747.8, 12920.0, 7313.4, 6610.0 /)
833 real :: nti(4), dp0gi(4)
834 ! seawater OC<0.2um concentration (in uM) - first is for low activity, the 2nd for high activity
835 ! High activity conc is from the average for RHaMBLe cruise in high activity regions
836 real, save :: oc02um(2) = (/ 0.0, 280.0 /)
837 ! Organic fraction for seasalt emissions (by size bin).
838 ! These are estimated from Figure 5 of Fuentes et al (ACP, 2011), assuming a value of
839 ! 2 ug/l for Chl-a for the high activity
840 real, save :: org_frac_low_activity(8) = &
841 (/ 0.05, 0.05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
842 real, save :: org_frac_high_activity(8) = &
843 (/ 0.2, 0.2, 0.1, 0.01, 0.01, 0.01, 0.01, 0.01 /)
846 p1st = PARAM_FIRST_SCALAR
848 ! for now just do itype=1
849 itype = 1
850 iphase = ai_phase
852 ! if using Feuntes et al (2011) then calculate the seawater OC dependent factors
853 if(seasalt_emiss_active .eq. 3 .or. seasalt_emiss_active .eq. 4)then
854 do i=1,4
855 nti(i) = beta_f1(i) * oc02um(seasalt_emiss_active-2)**2.0 + beta_f2(i) &
856 * oc02um(seasalt_emiss_active-2) + beta_f3(i)
857 dp0gi(i) = alpha_f1(i) + alpha_f2(i) &
858 * exp(-alpha_f3(i)*oc02um(seasalt_emiss_active-2))
859 end do
860 end if
863 ! compute emissions factors for each size bin
864 do n = 1, nsize_aer(itype)
865 ! changed the lower bound of the emission size limit to match lowest section edge, Qing.Yang@pnl.gov
866 ! DL - 30/3/2012 - select emission scheme (1=original, 3&4=Feuntes et al, 2011)
867 if(seasalt_emiss_active == 1)then
868 ! changed the lower bound of the emission size limit to match lowest section edge, Qing.Yang@pnl.gov
869 ! dumdlo = max( dlo_sect(n,itype), 0.1e-4 )
870 ! dumdhi = max( dhi_sect(n,itype), 0.1e-4 )
871 dumdlo = max( dlo_sect(n,itype), 0.02e-4 )
872 dumdhi = max( dhi_sect(n,itype), 0.02e-4 )
873 call seasalt_emitfactors_1bin( 1, dumdlo, dumdhi, &
874 ssemfact_numb(n,itype), dum, ssemfact_mass(n,itype) )
875 elseif(seasalt_emiss_active .eq. 3 .or. seasalt_emiss_active .eq. 4)then
876 call seasalt_emit_feuntes_1bin( dlo_sect(n,itype), dhi_sect(n,itype), &
877 ssemfact_numb(n,itype), dum, ssemfact_mass(n,itype), nti, dp0gi, oc02um(seasalt_emiss_active-2) )
878 endif
880 ! convert mass emissions factor from (g/m2/s) to (ug/m2/s)
881 ssemfact_mass(n,itype) = ssemfact_mass(n,itype)*1.0e6
882 end do
884 seas_flux(:,:) = 0.0
886 ! loop over i,j and apply seasalt emissions
887 k = kts
888 do j = jts, jte
889 do i = its, ite
891 ! Skip this point if over land. xland=1 for land and 2 for water.
892 ! Also, there is no way to differentiate fresh from salt water.
893 ! Currently, this assumes all water is salty.
894 if( xland(i,j) < 1.5 ) cycle
896 !wig: As far as I can tell, only real.exe knows the fractional breakdown
897 ! of land use. So, in wrf.exe, dumoceanfrac will always be 1.
898 dumoceanfrac = 1. !fraction of grid i,j that is salt water
899 dumspd10 = dumoceanfrac* &
900 ( (u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))**(0.5*3.41) )
902 ! factaa is (s*m2/kg-air)
903 ! factaa*ssemfact_mass(n) is (s*m2/kg-air)*(ug/m2/s) = ug/kg-air
904 ! factaa*ssemfact_numb(n) is (s*m2/kg-air)*( #/m2/s) = #/kg-air
905 factaa = (dtstep/dz8w(i,k,j))*alt(i,k,j)
907 factbb = factaa * dumspd10
909 if(seasalt_emiss_active == 1)then
911 seas_flux(i,j) = dumspd10*SUM(ssemfact_mass(1:nsize_aer(itype),itype))
913 ! apportion seasalt mass emissions assumming that seasalt is pure nacl
914 fracna = mw_na_aer / (mw_na_aer + mw_cl_aer)
915 fraccl = 1.0 - fracna
917 do n = 1, nsize_aer(itype)
919 ! only apply emissions if bin has both na and cl species
920 l_na = lptr_na_aer(n,itype,iphase)
921 l_cl = lptr_cl_aer(n,itype,iphase)
922 if ((l_na >= p1st) .and. (l_cl >= p1st)) then
924 chem(i,k,j,l_na) = chem(i,k,j,l_na) + &
925 factbb * ssemfact_mass(n,itype) * fracna
927 chem(i,k,j,l_cl) = chem(i,k,j,l_cl) + &
928 factbb * ssemfact_mass(n,itype) * fraccl
930 l = numptr_aer(n,itype,iphase)
931 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + &
932 factbb * ssemfact_numb(n,itype)
934 end if
936 end do ! n = 1, nsize_aer(itype)
939 elseif(seasalt_emiss_active.eq.3 .or. seasalt_emiss_active.eq.4)then
940 do n = 1, nsize_aer(itype)
942 ! apportion seasalt mass emissions assumming that seasalt is a
943 ! simple mix of pure nacl and a single primary organic
944 if(seasalt_emiss_active.eq.3)then
945 fracorg = org_frac_low_activity(n)
946 elseif(seasalt_emiss_active.eq.4)then
947 fracorg = org_frac_high_activity(n)
948 endif
949 fracna = mw_na_aer / (mw_na_aer + mw_cl_aer)
950 fraccl = 1.0 - fracna
951 fracna = (1.0-fracorg)*fracna
952 fraccl = (1.0-fracorg)*fraccl
954 ! only apply emissions if bin has both na and cl species
955 l_na = lptr_na_aer(n,itype,iphase)
956 l_cl = lptr_cl_aer(n,itype,iphase)
957 l_oc = lptr_oc_aer(n,itype,iphase)
958 if ((l_na >= p1st) .and. (l_cl >= p1st) .and. (l_oc >= p1st)) then
960 chem(i,k,j,l_na) = chem(i,k,j,l_na) + &
961 factbb * ssemfact_mass(n,itype) * fracna
963 chem(i,k,j,l_cl) = chem(i,k,j,l_cl) + &
964 factbb * ssemfact_mass(n,itype) * fraccl
966 chem(i,k,j,l_oc) = chem(i,k,j,l_oc) + &
967 factbb * ssemfact_mass(n,itype) * fracorg
969 l = numptr_aer(n,itype,iphase)
970 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + &
971 factbb * ssemfact_numb(n,itype)
973 end if
975 end do ! n = 1, nsize_aer(itype)
976 endif
978 end do ! i = its, ite
979 end do ! j = jts, jte
981 return
983 END subroutine mosaic_seasalt_emiss
987 !c----------------------------------------------------------------------
988 !c following is from gong06b.f in
989 !c /net/cirrus/files1/home/rce/oldfiles1/box/seasaltg
990 !c----------------------------------------------------------------------
991 subroutine seasalt_emitfactors_1bin( ireduce_smallr_emit, &
992 dpdrylo_cm, dpdryhi_cm, &
993 emitfact_numb, emitfact_surf, emitfact_mass )
994 !c
995 !c computes seasalt emissions factors for a specifed
996 !c dry particle size range
997 !c dpdrylo_cm = lower dry diameter (cm)
998 !c dpdryhi_cm = upper dry diameter (cm)
999 !c
1000 !c number and mass emissions are then computed as
1001 !c number emissions (#/m2/s) == emitfact_numb * (spd10*3.41)
1002 !c dry-sfc emissions (cm2/m2/s) == emitfact_surf * (spd10*3.41)
1003 !c dry-mass emissions (g/m2/s) == emitfact_mass * (spd10*3.41)
1004 !c
1005 !c where spd10 = 10 m windspeed in m/s
1006 !c
1007 !c uses bubble emissions formula (eqn 5a) from
1008 !c Gong et al. [JGR, 1997, p 3805-3818]
1009 !c
1010 !c *** for rdry < rdry_star, this formula overpredicts emissions.
1011 !c A strictly ad hoc correction is applied to the formula,
1012 !c based on sea-salt size measurements of
1013 !c O'Dowd et al. [Atmos Environ, 1997, p 73-80]
1014 !c
1015 !c *** the correction is only applied when ireduce_smallr_emit > 0
1016 !c
1017 implicit none
1019 !c subr arguments
1020 integer ireduce_smallr_emit
1021 real dpdrylo_cm, dpdryhi_cm, &
1022 emitfact_numb, emitfact_surf, emitfact_mass
1024 !c local variables
1025 integer isub_bin, nsub_bin
1027 real alnrdrylo
1028 real drydens, drydens_43pi_em12, x_4pi_em8
1029 real dum, dumadjust, dumb, dumexpb
1030 real dumsum_na, dumsum_ma, dumsum_sa
1031 real drwet, dlnrdry
1032 real df0drwet, df0dlnrdry, df0dlnrdry_star
1033 real relhum
1034 real rdry, rdrylo, rdryhi, rdryaa, rdrybb
1035 real rdrylowermost, rdryuppermost, rdry_star
1036 real rwet, rwetaa, rwetbb
1037 real rdry_cm, rwet_cm
1038 real sigmag_star
1039 real xmdry, xsdry
1041 real pi
1042 parameter (pi = 3.1415936536)
1044 !c c1-c4 are constants for seasalt hygroscopic growth parameterization
1045 !c in Eqn 3 and Table 2 of Gong et al. [1997]
1046 real c1, c2, c3, c4, onethird
1047 parameter (c1 = 0.7674)
1048 parameter (c2 = 3.079)
1049 parameter (c3 = 2.573e-11)
1050 parameter (c4 = -1.424)
1051 parameter (onethird = 1.0/3.0)
1054 !c dry particle density (g/cm3)
1055 drydens = 2.165
1056 !c factor for radius (micrometers) to mass (g)
1057 drydens_43pi_em12 = drydens*(4.0/3.0)*pi*1.0e-12
1058 !c factor for radius (micrometers) to surface (cm2)
1059 x_4pi_em8 = 4.0*pi*1.0e-8
1060 !c bubble emissions formula assume 80% RH
1061 relhum = 0.80
1063 !c rdry_star = dry radius (micrometers) below which the
1064 !c dF0/dr emission formula is adjusted downwards
1065 rdry_star = 0.1
1066 if (ireduce_smallr_emit .le. 0) rdry_star = -1.0e20
1067 !c sigmag_star = geometric standard deviation used for
1068 !c rdry < rdry_star
1069 sigmag_star = 1.9
1071 !c initialize sums
1072 dumsum_na = 0.0
1073 dumsum_sa = 0.0
1074 dumsum_ma = 0.0
1076 !c rdrylowermost, rdryuppermost = lower and upper
1077 !c dry radii (micrometers) for overall integration
1078 rdrylowermost = dpdrylo_cm*0.5e4
1079 rdryuppermost = dpdryhi_cm*0.5e4
1081 !c
1082 !c "section 1"
1083 !c integrate over rdry > rdry_star, where the dF0/dr emissions
1084 !c formula is applicable
1085 !c (when ireduce_smallr_emit <= 0, rdry_star = -1.0e20,
1086 !c and the entire integration is done here)
1087 !c
1088 if (rdryuppermost .le. rdry_star) goto 2000
1090 !c rdrylo, rdryhi = lower and upper dry radii (micrometers)
1091 !c for this part of the integration
1092 rdrylo = max( rdrylowermost, rdry_star )
1093 rdryhi = rdryuppermost
1095 nsub_bin = 1000
1097 alnrdrylo = log( rdrylo )
1098 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
1100 !c compute rdry, rwet (micrometers) at lowest size
1101 rdrybb = exp( alnrdrylo )
1102 rdry_cm = rdrybb*1.0e-4
1103 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
1104 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
1105 rwetbb = rwet_cm*1.0e4
1107 do 1900 isub_bin = 1, nsub_bin
1109 !c rdry, rwet at sub_bin lower boundary are those
1110 !c at upper boundary of previous sub_bin
1111 rdryaa = rdrybb
1112 rwetaa = rwetbb
1114 !c compute rdry, rwet (micrometers) at sub_bin upper boundary
1115 dum = alnrdrylo + isub_bin*dlnrdry
1116 rdrybb = exp( dum )
1118 rdry_cm = rdrybb*1.0e-4
1119 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
1120 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
1121 rwetbb = rwet_cm*1.0e4
1123 !c geometric mean rdry, rwet (micrometers) for sub_bin
1124 rdry = sqrt(rdryaa * rdrybb)
1125 rwet = sqrt(rwetaa * rwetbb)
1126 drwet = rwetbb - rwetaa
1128 !c xmdry is dry mass in g
1129 xmdry = drydens_43pi_em12 * (rdry**3.0)
1131 !c xsdry is dry surface in cm2
1132 xsdry = x_4pi_em8 * (rdry**2.0)
1134 !c dumb is "B" in Gong's Eqn 5a
1135 !c df0drwet is "dF0/dr" in Gong's Eqn 5a
1136 dumb = ( 0.380 - log10(rwet) ) / 0.650
1137 dumexpb = exp( -dumb*dumb)
1138 df0drwet = 1.373 * (rwet**(-3.0)) * &
1139 (1.0 + 0.057*(rwet**1.05)) * &
1140 (10.0**(1.19*dumexpb))
1142 dumsum_na = dumsum_na + drwet*df0drwet
1143 dumsum_ma = dumsum_ma + drwet*df0drwet*xmdry
1144 dumsum_sa = dumsum_sa + drwet*df0drwet*xsdry
1146 1900 continue
1149 !c
1150 !c "section 2"
1151 !c integrate over rdry < rdry_star, where the dF0/dr emissions
1152 !c formula is just an extrapolation and predicts too many emissions
1153 !c
1154 !c 1. compute dF0/dln(rdry) = (dF0/drwet)*(drwet/dlnrdry)
1155 !c at rdry_star
1156 !c 2. for rdry < rdry_star, assume dF0/dln(rdry) is lognormal,
1157 !c with the same lognormal parameters observed in
1158 !c O'Dowd et al. [1997]
1159 !c
1160 2000 if (rdrylowermost .ge. rdry_star) goto 3000
1162 !c compute dF0/dln(rdry) at rdry_star
1163 rdryaa = 0.99*rdry_star
1164 rdry_cm = rdryaa*1.0e-4
1165 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
1166 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
1167 rwetaa = rwet_cm*1.0e4
1169 rdrybb = 1.01*rdry_star
1170 rdry_cm = rdrybb*1.0e-4
1171 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
1172 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
1173 rwetbb = rwet_cm*1.0e4
1175 rwet = 0.5*(rwetaa + rwetbb)
1176 dumb = ( 0.380 - log10(rwet) ) / 0.650
1177 dumexpb = exp( -dumb*dumb)
1178 df0drwet = 1.373 * (rwet**(-3.0)) * &
1179 (1.0 + 0.057*(rwet**1.05)) * &
1180 (10.0**(1.19*dumexpb))
1182 drwet = rwetbb - rwetaa
1183 dlnrdry = log( rdrybb/rdryaa )
1184 df0dlnrdry_star = df0drwet * (drwet/dlnrdry)
1187 !c rdrylo, rdryhi = lower and upper dry radii (micrometers)
1188 !c for this part of the integration
1189 rdrylo = rdrylowermost
1190 rdryhi = min( rdryuppermost, rdry_star )
1192 nsub_bin = 1000
1194 alnrdrylo = log( rdrylo )
1195 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
1197 do 2900 isub_bin = 1, nsub_bin
1199 !c geometric mean rdry (micrometers) for sub_bin
1200 dum = alnrdrylo + (isub_bin-0.5)*dlnrdry
1201 rdry = exp( dum )
1203 !c xmdry is dry mass in g
1204 xmdry = drydens_43pi_em12 * (rdry**3.0)
1206 !c xsdry is dry surface in cm2
1207 xsdry = x_4pi_em8 * (rdry**2.0)
1209 !c dumadjust is adjustment factor to reduce dF0/dr
1210 dum = log( rdry/rdry_star ) / log( sigmag_star )
1211 dumadjust = exp( -0.5*dum*dum )
1213 df0dlnrdry = df0dlnrdry_star * dumadjust
1215 dumsum_na = dumsum_na + dlnrdry*df0dlnrdry
1216 dumsum_ma = dumsum_ma + dlnrdry*df0dlnrdry*xmdry
1217 dumsum_sa = dumsum_sa + dlnrdry*df0dlnrdry*xsdry
1219 2900 continue
1222 !c
1223 !c all done
1224 !c
1225 3000 emitfact_numb = dumsum_na
1226 emitfact_mass = dumsum_ma
1227 emitfact_surf = dumsum_sa
1229 return
1230 end subroutine seasalt_emitfactors_1bin
1233 !c----------------------------------------------------------------------
1234 !c Following is an adaption of the subroutine above.
1235 !c It has been modified to include flux from Fuentes et al, ACP, 2010
1236 !c for dealing with smaller particle emissions.
1237 !c----------------------------------------------------------------------
1238 subroutine seasalt_emit_feuntes_1bin( &
1239 dpdrylo_cm, dpdryhi_cm, &
1240 emitfact_numb, emitfact_surf, emitfact_mass, nti, dp0gi, oc02um )
1241 !c
1242 !c computes seasalt emissions factors for a specifed
1243 !c dry particle size range
1244 !c dpdrylo_cm = lower dry diameter (cm)
1245 !c dpdryhi_cm = upper dry diameter (cm)
1246 !c
1247 !c number and mass emissions are then computed as
1248 !c number emissions (#/m2/s) == emitfact_numb * (spd10*3.41)
1249 !c dry-sfc emissions (cm2/m2/s) == emitfact_surf * (spd10*3.41)
1250 !c dry-mass emissions (g/m2/s) == emitfact_mass * (spd10*3.41)
1251 !c
1252 !c where spd10 = 10 m windspeed in m/s
1253 !c
1254 !c Uses bubble emissions formula (eqn 5a) from
1255 !c Gong et al. [JGR, 1997, p 3805-3818] for particles
1256 !c with dry diameters of 0.45 um or more
1257 !c
1258 !c For smaller particles we use the parameterisation of
1259 !c Fuentes et al, ACP, 2010.
1260 !c
1263 implicit none
1265 ! subr arguments
1266 !integer ireduce_smallr_emit
1267 real dpdrylo_cm, dpdryhi_cm, &
1268 emitfact_numb, emitfact_surf, emitfact_mass
1269 real, intent(in) :: nti(4), dp0gi(4), oc02um
1271 ! local variables
1272 integer isub_bin, nsub_bin, jd, nsub_lower_bin, nsub_upper_bin
1274 real alnrdrylo
1275 real drydens, drydens_43pi_em12, x_4pi_em8, drydens_16pi_em12, x_pi_em8
1276 real dum, dumadjust, dumb, dumexpb
1277 real dumsum_na, dumsum_ma, dumsum_sa
1278 real drwet, dlnrdry, ddwet, ddry, dwet
1279 real df0drwet, df0dlnrdry, df0dlnrdry_star
1280 real relhum
1281 real rdry, rdrylo, rdryhi, rdryaa, rdrybb
1282 real rdrylowermost, rdryuppermost, rdry_star
1283 real rwet, rwetaa, rwetbb
1284 real rdry_cm, rwet_cm
1285 real sigmag_star
1286 real xmdry, xsdry
1288 real ddrylo, ddryhi, alogddrylo
1289 real ddrybb, ddry_cm, dwet_cm, dwetbb
1290 real ddryaa, dwetaa, dlogddry, df0dlogddry
1292 real pi
1293 parameter (pi = 3.1415936536)
1295 ! c1-c4 are constants for seasalt hygroscopic growth parameterization
1296 ! in Eqn 3 and Table 2 of Gong et al. [1997]
1297 real c1, c2, c3, c4, onethird
1298 parameter (c1 = 0.7674)
1299 parameter (c2 = 3.079)
1300 parameter (c3 = 2.573e-11)
1301 parameter (c4 = -1.424)
1302 parameter (onethird = 1.0/3.0)
1305 ! constants for seasalt/organic aerosol distribution parameterisation
1306 ! from Fuentes et al, ACP, 2010
1307 real, save :: width_ssinorg_f(4) = &
1308 (/ 1.55, 1.7, 1.5, 1.7 /)
1309 real, save :: width_ssorg_f(4) = &
1310 (/ 1.55, 1.9, 1.5, 1.7 /)
1311 real :: width_ss_f(4), frac
1312 ! scale_factor is a combination of:
1313 ! 1) (Q) sweep air flow (58.3 cm3/s)
1314 ! 2) (Ab) total surface area covered by bubbles (0.0146 m2)
1315 ! 4) scaling factor of 3.84e-6 from the whitecap coverage formulation of
1316 ! Monahan and O'Muircheartaigh (1980) - as the Gong et al formulation
1317 ! includes a whitecap coverage factor too
1318 ! this is for converting from dNt to dFp for Fuentes et al (ACP, 2010)
1319 real, parameter :: scale_factor = 58.3/(0.0146)*3.84e-6 !0.01533
1322 ! select which distribution width to use for Fuentes et al,
1323 if(oc02um.gt.0e0)then
1324 width_ss_f = width_ssorg_f
1325 else
1326 width_ss_f = width_ssinorg_f
1327 end if
1330 ! dry particle density (g/cm3)
1331 drydens = 2.165
1332 ! factor for radius (micrometers) to mass (g)
1333 drydens_43pi_em12 = drydens*(4.0/3.0)*pi*1.0e-12
1334 ! factor for diameter (micrometers) to mass (g)
1335 drydens_16pi_em12 = drydens*(1.0/6.0)*pi*1.0e-12
1336 ! factor for radius (micrometers) to surface (cm2)
1337 x_4pi_em8 = 4.0*pi*1.0e-8
1338 ! factor for diameter (micrometers) to surface (cm2)
1339 x_pi_em8 = pi*1.0e-8
1340 ! bubble emissions formula assume 80% RH
1341 relhum = 0.80
1343 ! rdry_star = dry radius (micrometers) below which the
1344 ! dF0/dr emission formula from Fuentes et al, ACP, 2010
1345 ! will be used.
1346 rdry_star = 0.45 / 2.0
1347 !if (ireduce_smallr_emit .le. 0) rdry_star = -1.0e20
1348 ! sigmag_star = geometric standard deviation used for
1349 ! rdry < rdry_star
1350 sigmag_star = 1.9
1352 ! initialize sums
1353 dumsum_na = 0.0
1354 dumsum_sa = 0.0
1355 dumsum_ma = 0.0
1357 ! rdrylowermost, rdryuppermost = lower and upper
1358 ! dry radii (micrometers) for overall integration
1359 rdrylowermost = dpdrylo_cm*0.5e4
1360 rdryuppermost = dpdryhi_cm*0.5e4
1363 ! "section 1"
1364 ! integrate over rdry > rdry_star, where the dF0/dr emissions
1365 ! formula from the original subroutine is applicable
1366 ! (so using Gong et al for all emission profile)
1367 if (rdrylowermost .ge. rdry_star) then
1369 ! rdrylo, rdryhi = lower and upper dry radii (micrometers)
1370 ! for this part of the integration
1371 rdrylo = rdrylowermost
1372 rdryhi = rdryuppermost
1374 nsub_bin = 1000
1376 alnrdrylo = log( rdrylo )
1377 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
1379 ! compute rdry, rwet (micrometers) at lowest size
1380 rdrybb = exp( alnrdrylo )
1381 rdry_cm = rdrybb*1.0e-4
1382 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
1383 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
1384 rwetbb = rwet_cm*1.0e4
1386 do isub_bin = 1, nsub_bin
1388 ! rdry, rwet at sub_bin lower boundary are those
1389 ! at upper boundary of previous sub_bin
1390 rdryaa = rdrybb
1391 rwetaa = rwetbb
1393 ! compute rdry, rwet (micrometers) at sub_bin upper boundary
1394 dum = alnrdrylo + isub_bin*dlnrdry
1395 rdrybb = exp( dum )
1397 rdry_cm = rdrybb*1.0e-4
1398 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
1399 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
1400 rwetbb = rwet_cm*1.0e4
1402 ! geometric mean rdry, rwet (micrometers) for sub_bin
1403 rdry = sqrt(rdryaa * rdrybb)
1404 rwet = sqrt(rwetaa * rwetbb)
1405 drwet = rwetbb - rwetaa
1407 ! xmdry is dry mass in g
1408 xmdry = drydens_43pi_em12 * (rdry**3.0)
1410 ! xsdry is dry surface in cm2
1411 xsdry = x_4pi_em8 * (rdry**2.0)
1413 ! dumb is "B" in Gong's Eqn 5a
1414 ! df0drwet is "dF0/dr" in Gong's Eqn 5a
1415 dumb = ( 0.380 - log10(rwet) ) / 0.650
1416 dumexpb = exp( -dumb*dumb)
1417 df0drwet = 1.373 * (rwet**(-3.0)) * &
1418 (1.0 + 0.057*(rwet**1.05)) * &
1419 (10.0**(1.19*dumexpb))
1421 dumsum_na = dumsum_na + drwet*df0drwet
1422 dumsum_ma = dumsum_ma + drwet*df0drwet*xmdry
1423 dumsum_sa = dumsum_sa + drwet*df0drwet*xsdry
1425 end do
1428 ! "section 2"
1429 ! integrate over rdry < rdry_star for old parameterisation, and use
1430 ! Fuentes et al below that
1431 !
1432 ! 1. for rdry < rdry_star, use Fuentes et al parameterisation
1433 ! 2. for rdry > rdry_star, use Gong et al parameterisation, as above
1434 elseif (rdryuppermost .gt. rdry_star) then
1436 ! determine the fraction of size bin below rdry_star
1437 frac = (log(rdry_star)-log(rdrylowermost)) / (log(rdryuppermost)-log(rdrylowermost))
1438 nsub_lower_bin = floor(frac*1000.0) ! calc number of size bins to use for section below rdry_star
1439 nsub_upper_bin = 1000-nsub_lower_bin ! use remaining size bins for section above rdry_star
1442 ! 1.................
1443 ! ddrylo, ddryhi = lower and upper dry diameter (micrometers)
1444 ! for this part of the integration
1445 ddrylo = rdrylowermost*2.0
1446 ddryhi = rdry_star*2.0
1448 nsub_bin = nsub_lower_bin
1450 alogddrylo = log10( ddrylo )
1451 dlogddry = (log10( ddryhi ) - alogddrylo)/nsub_bin
1453 ! compute ddry, dwet (micrometers) at lowest size
1454 ddrybb = 10.0**( alogddrylo )
1455 ddry_cm = ddrybb*1.0e-4
1456 dwet_cm = ( ddry_cm**3 + (c1*(ddry_cm**c2))/ &
1457 ( (c3*(ddry_cm**c4)) - log10(relhum) ) )**onethird
1458 dwetbb = dwet_cm*1.0e4
1460 do isub_bin = 1, nsub_bin
1462 ! ddry, dwet at sub_bin lower boundary are those
1463 ! at upper boundary of previous sub_bin
1464 ddryaa = ddrybb
1465 !dwetaa = dwetbb
1467 ! compute ddry, dwet (micrometers) at sub_bin upper boundary
1468 dum = alogddrylo + isub_bin*dlogddry
1469 ddrybb = 10.0**( dum )
1471 ddry_cm = ddrybb*1.0e-4
1472 !dwet_cm = ( ddry_cm**3 + (c1*(ddry_cm**c2))/ &
1473 ! ( (c3*(ddry_cm**c4)) - log10(relhum) ) )**onethird
1474 !dwetbb = dwet_cm*1.0e4
1476 ! geometric mean rdry, rwet (micrometers) for sub_bin
1477 ddry = sqrt(ddryaa * ddrybb)
1478 !dwet = sqrt(dwetaa * dwetbb)
1480 ! xmdry is dry mass in g
1481 xmdry = drydens_16pi_em12 * (ddry**3.0)
1483 ! xsdry is dry surface in cm2
1484 xsdry = x_pi_em8 * (ddry**2.0)
1486 ! Equation 2 from Fuentes et al (ACP, 2011). Wet diameter is scaled by a factor
1487 ! of 1e3 to convert from um to nm for this calculation
1488 df0dlogddry = 0.0
1489 do jd = 1, 4
1490 df0dlogddry = df0dlogddry + nti(jd)/( (2.0*pi)**0.5 * log10(width_ss_f(jd)) ) &
1491 * exp(-1.0/2.0 * (log10(ddry*1e3/dp0gi(jd))/log10(width_ss_f(jd)))**2.0 )
1492 end do
1494 dumsum_na = dumsum_na + dlogddry*df0dlogddry*scale_factor
1495 dumsum_ma = dumsum_ma + dlogddry*df0dlogddry*scale_factor*xmdry
1496 dumsum_sa = dumsum_sa + dlogddry*df0dlogddry*scale_factor*xsdry
1498 end do
1501 ! 2................
1502 ! rdrylo, rdryhi = lower and upper dry radii (micrometers)
1503 ! for this part of the integration
1504 rdrylo = rdry_star
1505 rdryhi = rdryuppermost
1507 nsub_bin = nsub_upper_bin
1509 alnrdrylo = log( rdrylo )
1510 dlnrdry = (log( rdryhi ) - alnrdrylo)/nsub_bin
1512 ! compute rdry, rwet (micrometers) at lowest size
1513 rdrybb = exp( alnrdrylo )
1514 rdry_cm = rdrybb*1.0e-4
1515 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
1516 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
1517 rwetbb = rwet_cm*1.0e4
1519 do isub_bin = 1, nsub_bin
1521 ! rdry, rwet at sub_bin lower boundary are those
1522 ! at upper boundary of previous sub_bin
1523 rdryaa = rdrybb
1524 rwetaa = rwetbb
1526 ! compute rdry, rwet (micrometers) at sub_bin upper boundary
1527 dum = alnrdrylo + isub_bin*dlnrdry
1528 rdrybb = exp( dum )
1530 rdry_cm = rdrybb*1.0e-4
1531 rwet_cm = ( rdry_cm**3 + (c1*(rdry_cm**c2))/ &
1532 ( (c3*(rdry_cm**c4)) - log10(relhum) ) )**onethird
1533 rwetbb = rwet_cm*1.0e4
1535 ! geometric mean rdry, rwet (micrometers) for sub_bin
1536 rdry = sqrt(rdryaa * rdrybb)
1537 rwet = sqrt(rwetaa * rwetbb)
1538 drwet = rwetbb - rwetaa
1540 ! xmdry is dry mass in g
1541 xmdry = drydens_43pi_em12 * (rdry**3.0)
1543 ! xsdry is dry surface in cm2
1544 xsdry = x_4pi_em8 * (rdry**2.0)
1546 ! dumb is "B" in Gong's Eqn 5a
1547 ! df0drwet is "dF0/dr" in Gong's Eqn 5a
1548 dumb = ( 0.380 - log10(rwet) ) / 0.650
1549 dumexpb = exp( -dumb*dumb)
1550 df0drwet = 1.373 * (rwet**(-3.0)) * &
1551 (1.0 + 0.057*(rwet**1.05)) * &
1552 (10.0**(1.19*dumexpb))
1554 dumsum_na = dumsum_na + drwet*df0drwet
1555 dumsum_ma = dumsum_ma + drwet*df0drwet*xmdry
1556 dumsum_sa = dumsum_sa + drwet*df0drwet*xsdry
1558 end do
1562 ! "section 3"
1563 ! where rdry < rdry_star for the whole bin we use Fuentes et al
1564 ! for all emissions
1565 !
1566 else
1568 ! ddrylo, ddryhi = lower and upper dry diameter (micrometers)
1569 ! for this part of the integration
1570 ddrylo = rdrylowermost*2.0
1571 ddryhi = rdryuppermost*2.0
1573 nsub_bin = 1000
1575 alogddrylo = log10( ddrylo )
1576 dlogddry = (log10( ddryhi ) - alogddrylo)/nsub_bin
1578 ! compute ddry, dwet (micrometers) at lowest size
1579 ddrybb = 10.0**( alogddrylo )
1580 ddry_cm = ddrybb*1.0e-4
1581 dwet_cm = ( ddry_cm**3 + (c1*(ddry_cm**c2))/ &
1582 ( (c3*(ddry_cm**c4)) - log10(relhum) ) )**onethird
1583 dwetbb = dwet_cm*1.0e4
1585 do isub_bin = 1, nsub_bin
1587 ! ddry, dwet at sub_bin lower boundary are those
1588 ! at upper boundary of previous sub_bin
1589 ddryaa = ddrybb
1590 !dwetaa = dwetbb
1592 ! compute ddry, dwet (micrometers) at sub_bin upper boundary
1593 dum = alogddrylo + isub_bin*dlogddry
1594 ddrybb = 10.0**( dum )
1596 ddry_cm = ddrybb*1.0e-4
1597 !dwet_cm = ( ddry_cm**3 + (c1*(ddry_cm**c2))/ &
1598 ! ( (c3*(ddry_cm**c4)) - log10(relhum) ) )**onethird
1599 !dwetbb = dwet_cm*1.0e4
1601 ! geometric mean rdry, rwet (micrometers) for sub_bin
1602 ddry = sqrt(ddryaa * ddrybb)
1603 !dwet = sqrt(dwetaa * dwetbb)
1605 ! xmdry is dry mass in g
1606 xmdry = drydens_16pi_em12 * (ddry**3.0)
1608 ! xsdry is dry surface in cm2
1609 xsdry = x_pi_em8 * (ddry**2.0)
1611 ! Equation 2 from Fuentes et al (ACP, 2011). Wet diameter is scaled by a factor
1612 ! of 1e3 to convert from um to nm for this calculation
1613 df0dlogddry = 0.0
1614 do jd = 1, 4
1615 df0dlogddry = df0dlogddry + nti(jd)/( (2.0*pi)**0.5 * log10(width_ss_f(jd)) ) &
1616 * exp(-1.0/2.0 * (log10(ddry*1e3/dp0gi(jd))/log10(width_ss_f(jd)))**2.0 )
1617 end do
1619 dumsum_na = dumsum_na + dlogddry*df0dlogddry*scale_factor
1620 dumsum_ma = dumsum_ma + dlogddry*df0dlogddry*scale_factor*xmdry
1621 dumsum_sa = dumsum_sa + dlogddry*df0dlogddry*scale_factor*xsdry
1623 end do
1627 end if
1628 !c
1629 !c all done
1630 !c
1631 emitfact_numb = dumsum_na
1632 emitfact_mass = dumsum_ma
1633 emitfact_surf = dumsum_sa
1636 return
1637 end subroutine seasalt_emit_feuntes_1bin
1644 END MODULE module_mosaic_addemiss
1646 !----------------------------------------------------------------------
1648 subroutine mosaic_dust_emiss( slai,ust, smois, ivgtyp, isltyp, &
1649 id, dtstep, u10, v10, alt, dz8w, xland, config_flags, chem, &
1650 ids,ide, jds,jde, kds,kde, &
1651 ims,ime, jms,jme, kms,kme, &
1652 its,ite, jts,jte, kts,kte )
1653 !
1654 ! adds dust emissions for mosaic aerosol species (i.e. emission tendencies
1655 ! over time dtstep are applied to the aerosol mixing ratios)
1656 !
1657 ! This is a simple dust scheme based on Shaw et al. (2008) to appear in
1658 ! Atmospheric Environment, recoded by Jerome Fast
1659 !
1660 ! NOTE:
1661 ! 1) This version only works with the 8-bin version of MOSAIC.
1662 ! 2) Dust added to MOSAIC's other inorganic specie, OIN. If Ca and CO3 are
1663 ! activated in the Registry, a small fraction also added to Ca and CO3.
1664 ! 3) The main departure from Shaw et al., is now alphamask is computed since
1665 ! the land-use categories in that paper and in WRF differ. WRF currently
1666 ! does not have that many land-use categories and adhoc assumptions had to
1667 ! be made. This version was tested for Mexico in the dry season. The main
1668 ! land-use categories in WRF that are likely dust sources are grass, shrub,
1669 ! and savannna (that WRF has in the desert regions of NW Mexico). Having
1670 ! dust emitted from these types for other locations and other times of the
1671 ! year is not likely to be valid.
1672 ! 4) An upper bound on ustar was placed because the surface parameterizations
1673 ! in WRF can produce unrealistically high values that lead to very high
1674 ! dust emission rates.
1675 ! 5) Other departures' from Shaw et al. noted below, but are probably not as
1676 ! important as 2) and 3).
1677 ! 6) Now the dust added into dust species
1679 USE module_configure, only: grid_config_rec_type
1680 USE module_state_description, only: num_chem, param_first_scalar
1681 USE module_data_mosaic_asect
1683 IMPLICIT NONE
1685 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1687 INTEGER, INTENT(IN ) :: id, &
1688 ids,ide, jds,jde, kds,kde, &
1689 ims,ime, jms,jme, kms,kme, &
1690 its,ite, jts,jte, kts,kte
1692 REAL, INTENT(IN ) :: dtstep
1694 ! 10-m wind speed components (m/s)
1695 REAL, DIMENSION( ims:ime , jms:jme ), &
1696 INTENT(IN ) :: u10, v10, xland, slai, ust
1697 INTEGER, DIMENSION( ims:ime , jms:jme ), &
1698 INTENT(IN ) :: ivgtyp, isltyp
1700 ! trace species mixing ratios (aerosol mass = ug/kg; number = #/kg)
1701 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
1702 INTENT(INOUT ) :: chem
1704 ! alt = 1.0/(dry air density) in (m3/kg)
1705 ! dz8w = layer thickness in (m)
1706 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
1707 INTENT(IN ) :: alt, dz8w
1709 REAL, DIMENSION( ims:ime, config_flags%num_soil_layers, jms:jme ) , &
1710 INTENT(INOUT) :: smois
1712 ! local variables
1713 integer i, j, k, l, l_oin, l_ca, l_co3, n, ii
1714 integer iphase, itype, izob
1715 integer p1st
1717 real dum, dumdlo, dumdhi, dumlandfrac, dumspd10
1718 real factaa, factbb, fracoin, fracca, fracco3, fractot
1719 real ustart, ustar1, ustart0
1720 real alphamask, f8, f50, f51, f52, wetfactor, sumdelta, ftot
1721 real smois_grav, wp, pclay
1722 real :: beta(4,7)
1723 real :: gamma(4), delta(4)
1724 real :: sz(8)
1725 real :: dustflux, densdust, mass1part
1726 real :: c_const
1728 p1st = param_first_scalar
1729 !
1730 ! from Nickovic et al., JGR, 2001 and Shaw et al. 2007
1731 ! beta: fraction of clay, small silt, large silt, and sand correcsponding to Zobler class (7)
1732 ! beta (1,*) for 0.5-1 um
1733 ! beta (2,*) for 1-10 um
1734 ! beta (3,*) for 10-25 um
1735 ! beta (4,*) for 25-50 um
1736 !
1737 beta(1,1)=0.12
1738 beta(2,1)=0.04
1739 beta(3,1)=0.04
1740 beta(4,1)=0.80
1741 beta(1,2)=0.34
1742 beta(2,2)=0.28
1743 beta(3,2)=0.28
1744 beta(4,2)=0.10
1745 beta(1,3)=0.45
1746 beta(2,3)=0.15
1747 beta(3,3)=0.15
1748 beta(4,3)=0.25
1749 beta(1,4)=0.12
1750 beta(2,4)=0.09
1751 beta(3,4)=0.09
1752 beta(4,4)=0.70
1753 beta(1,5)=0.40
1754 beta(2,5)=0.05
1755 beta(3,5)=0.05
1756 beta(4,5)=0.50
1757 beta(1,6)=0.34
1758 beta(2,6)=0.18
1759 beta(3,6)=0.18
1760 beta(4,6)=0.30
1761 beta(1,7)=0.22
1762 beta(2,7)=0.09
1763 beta(3,7)=0.09
1764 beta(4,7)=0.60
1765 gamma(1)=0.08
1766 gamma(2)=1.00
1767 gamma(3)=1.00
1768 gamma(4)=0.12
1769 !
1770 ! * Mass fractions for each size bin. These values were recommended by
1771 ! Natalie Mahowold, with bins 7 and 8 the same as bins 3 and 4 from CAM.
1772 ! * Changed slightly since Natelie's estimates do not add up to 1.0
1773 ! * This would need to be made more generic for other bin sizes.
1774 ! sz(1)=0
1775 ! sz(2)=1.78751e-06
1776 ! sz(3)=0.000273786
1777 ! sz(4)=0.00847978
1778 ! sz(5)=0.056055
1779 ! sz(6)=0.0951896
1780 ! sz(7)=0.17
1781 ! sz(8)=0.67
1782 ! sz(1)=1.78751e-06
1783 ! sz(2)=0.00875357
1784 ! sz(3)=0.1512446
1785 ! sz(4)=0.84
1786 itype=1
1787 if (nsize_aer(itype) .eq. 8) then
1788 !BSINGH(12/11/2013): Based on the suggestions by Manish Shrivastva, sz variable is modified below.
1789 !Original values are commented out
1790 !sz(1)=0.0
1791 !sz(2)=0.0
1792 !sz(3)=0.0005
1793 !sz(4)=0.0095
1794 !sz(5)=0.01
1795 !sz(6)=0.06
1796 !sz(7)=0.20
1797 !sz(8)=0.72
1799 sz(1)=0
1800 sz(2)=1.78751e-06
1801 sz(3)=0.000273786
1802 sz(4)=0.00847978
1803 sz(5)=0.056055
1804 sz(6)=0.0951896
1805 sz(7)=0.17
1806 sz(8)=0.67
1807 else if (nsize_aer(itype) .eq. 4) then
1808 !BSINGH(12/11/2013): Based on the suggestions by Manish Shrivastva, sz variable is modified below.
1809 !Original values are commented out
1810 !sz(1)=0.0
1811 !sz(2)=0.01
1812 !sz(3)=0.07
1813 !sz(4)=0.92
1814 sz(1)=1.78751e-06
1815 sz(2)=0.00875357
1816 sz(3)=0.1512446
1817 sz(4)=0.84
1818 sz(5)=0.0
1819 sz(6)=0.0
1820 sz(7)=0.0
1821 sz(8)=0.0
1822 endif
1824 ! for now just do itype=1
1825 itype = 1
1826 iphase = ai_phase
1828 ! loop over i,j and apply dust emissions
1829 k = kts
1830 do 1830 j = jts, jte
1831 do 1820 i = its, ite
1833 if( xland(i,j) > 1.5 ) cycle
1835 ! compute wind speed anyway, even though ustar is used below
1837 dumlandfrac = 1.
1838 dumspd10=(u10(i,j)*u10(i,j) + v10(i,j)*v10(i,j))**(0.5)
1839 if(dumspd10 >= 5.0) then
1840 dumspd10 = dumlandfrac* &
1841 ( dumspd10*dumspd10*(dumspd10-5.0))
1842 else
1843 dumspd10=0.
1844 endif
1846 ! part1 - compute vegetation mask
1847 !
1848 ! * f8, f50, f51, f52 refer to vegetation classes from the Olsen categories
1849 ! for desert, sand desert, grass semi-desert, and shrub semi-desert
1850 ! * in WRF, vegetation types 7, 8 and 10 are grassland, shrubland, and savanna
1851 ! that are dominate types in Mexico and probably have some erodable surface
1852 ! during the dry season
1853 ! * currently modified these values so that only a small fraction of cell
1854 ! area is erodable
1855 ! * these values are highly tuneable!
1857 alphamask=0.001
1858 if (ivgtyp(i,j) .eq. 7) then
1859 f8=0.005
1860 f50=0.00
1861 f51=0.10
1862 f51=0.066
1863 f52=0.00
1864 alphamask=(f8+f50)*1.0+(f51+f52)*0.5
1865 endif
1866 if (ivgtyp(i,j) .eq. 8) then
1867 f8=0.010
1868 f50=0.00
1869 f51=0.00
1870 f52=0.15
1871 f52=0.10
1872 alphamask=(f8+f50)*1.0+(f51+f52)*0.5
1873 endif
1874 if (ivgtyp(i,j) .eq. 10) then
1875 f8=0.00
1876 f50=0.00
1877 f51=0.01
1878 f52=0.00
1879 alphamask=(f8+f50)*1.0+(f51+f52)*0.5
1880 endif
1882 ! part2 - zobler
1883 !
1884 ! * in Shaw's paper, dust is computed for 4 size ranges:
1885 ! 0.5-1 um
1886 ! 1-10 um
1887 ! 10-25 um
1888 ! 25-50 um
1889 ! * Shaw's paper also accounts for sub-grid variability in soil
1890 ! texture, but here we just assume the same soil texture for each
1891 ! grid cell
1892 ! * since MOSAIC is currently has a maximum size range up to 10 um,
1893 ! neglect upper 2 size ranges and lowest size range (assume small)
1894 ! * map WRF soil classes arbitrarily to Zolber soil textural classes
1895 ! * skip dust computations for WRF soil classes greater than 13, i.e.
1896 ! do not compute dust over water, bedrock, and other surfaces
1897 ! * should be skipping for water surface at this point anyway
1898 !
1899 izob=0
1900 if(isltyp(i,j).eq.1) izob=1
1901 if(isltyp(i,j).eq.2) izob=1
1902 if(isltyp(i,j).eq.3) izob=4
1903 if(isltyp(i,j).eq.4) izob=2
1904 if(isltyp(i,j).eq.5) izob=2
1905 if(isltyp(i,j).eq.6) izob=2
1906 if(isltyp(i,j).eq.7) izob=7
1907 if(isltyp(i,j).eq.8) izob=2
1908 if(isltyp(i,j).eq.9) izob=6
1909 if(isltyp(i,j).eq.10) izob=5
1910 if(isltyp(i,j).eq.11) izob=2
1911 if(isltyp(i,j).eq.12) izob=3
1912 if(isltyp(i,j).ge.13) izob=0
1913 if(izob.eq.0) goto 1840
1914 !
1915 ! part3 - dustprod
1916 !
1917 do ii=1,4
1918 delta(ii)=0.0
1919 enddo
1920 sumdelta=0.0
1921 do ii=1,4
1922 delta(ii)=beta(ii,izob)*gamma(ii)
1923 if(ii.lt.4) then
1924 sumdelta=sumdelta+delta(ii)
1925 endif
1926 enddo
1927 do ii=1,4
1928 delta(ii)=delta(ii)/sumdelta
1929 enddo
1931 ! part4 - wetness
1932 !
1933 ! * assume dry for now, have passed in soil moisture to this routine
1934 ! but needs to be included here
1935 ! * wetfactor less than 1 would reduce dustflux
1936 ! * convert model soil moisture (m3/m3) to gravimetric soil moisture
1937 ! (mass of water / mass of soil in %) assuming a constant density
1938 ! for soil
1939 pclay=beta(1,izob)*100.
1940 wp=0.0014*pclay*pclay+0.17*pclay
1941 smois_grav=(smois(i,1,j)/2.6)*100.
1942 if(smois_grav.gt.wp) then
1943 wetfactor=sqrt(1.0+1.21*(smois_grav-wp)**0.68)
1944 else
1945 wetfactor=1.0
1946 endif
1947 ! wetfactor=1.0
1949 ! part5 - dustflux
1950 ! lower bound on ustar = 20 cm/s as in Shaw et al, but set upper
1951 ! bound to 100 cm/s
1952 c_const=1.e-14 ! default
1954 ustar1=ust(i,j)*100.0
1955 if(ustar1.gt.100.0) ustar1=100.0
1956 ustart0=20.0
1957 ustart=ustart0*wetfactor
1958 if(ustar1.le.ustart) then
1959 dustflux=0.0
1960 else
1961 dustflux=c_const*(ustar1**4)*(1.0-(ustart/ustar1))
1962 endif
1963 dustflux=dustflux*10.0
1964 ! units kg m-2 s-1
1965 ftot=0.0
1966 do ii=1,2
1967 ftot=ftot+dustflux*alphamask*delta(ii)
1968 enddo
1969 ! convert to ug m-2 s-1
1970 ftot=ftot*1.0e+09
1972 ! apportion other inorganics only
1973 factaa = (dtstep/dz8w(i,k,j))*alt(i,k,j)
1974 factbb = factaa * ftot
1975 fracoin = 0.97
1976 fracca = 0.03*0.4
1977 fracco3 = 0.03*0.6
1978 fractot = fracoin + fracca + fracco3
1979 ! if (ivgtyp(i,j) .eq. 8) print*,'jdf',i,j,ustar1,ustart0,factaa,ftot
1980 do 1810 n = 1, nsize_aer(itype)
1981 l_oin = lptr_oin_aer(n,itype,iphase)
1982 l_ca = lptr_ca_aer(n,itype,iphase)
1983 l_co3 = lptr_co3_aer(n,itype,iphase)
1984 if (l_oin >= p1st) then
1985 chem(i,k,j,l_oin) = chem(i,k,j,l_oin) + &
1986 factbb * sz(n) * fracoin
1987 if (l_ca >= p1st) &
1988 chem(i,k,j,l_ca) = chem(i,k,j,l_ca) + &
1989 factbb * sz(n) * fracca
1990 if (l_co3 >= p1st) &
1991 chem(i,k,j,l_co3) = chem(i,k,j,l_co3) + &
1992 factbb * sz(n) * fracco3
1993 ! mass1part is mass of a single particle in ug, density of dust ~2.5 g cm-3
1994 densdust=2.5
1995 mass1part=0.523598*(dcen_sect(n,itype)**3)*densdust*1.0e06
1996 l = numptr_aer(n,itype,iphase)
1997 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + &
1998 factbb * sz(n) * fractot / mass1part
1999 end if
2000 1810 continue
2002 1840 continue
2004 1820 continue
2005 1830 continue
2007 return
2009 END subroutine mosaic_dust_emiss
2011 !====================================================================================
2012 ! add another dust emission scheme following GOCART mechanism --czhao 09/17/2009
2013 ! Ukhov A. 02/28/2022. Code clean up and refactoring. Removed dublicate mosaic_source_du().
2014 ! Instead original source_du() is called from module_gocart_dust.F
2015 !====================================================================================
2016 subroutine mosaic_dust_gocartemis (dt,start_month,num_soil_layers,alt, &
2017 u_phy,v_phy,chem,rho_phy,dz8w,smois,u10,v10,erod, &
2018 isltyp,xland,g, &
2019 emis_dust, &
2020 ids,ide, jds,jde, kds,kde, &
2021 ims,ime, jms,jme, kms,kme, &
2022 its,ite, jts,jte, kts,kte)
2023 USE module_data_gocart_dust
2024 USE module_configure
2025 USE module_state_description
2026 USE module_model_constants, ONLY: mwdry
2027 USE module_data_mosaic_asect
2028 USE gocart_dust, ONLY: source_du
2029 IMPLICIT NONE
2031 INTEGER, INTENT(IN ) :: start_month,num_soil_layers, &
2032 ids,ide, jds,jde, kds,kde, &
2033 ims,ime, jms,jme, kms,kme, &
2034 its,ite, jts,jte, kts,kte
2035 INTEGER,DIMENSION( ims:ime , jms:jme ) , &
2036 INTENT(IN ) :: isltyp
2038 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
2039 INTENT(INOUT ) :: chem
2041 REAL, DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL, &
2042 INTENT(INOUT ) :: emis_dust
2045 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
2046 INTENT(INOUT) :: smois
2047 REAL, DIMENSION( ims:ime , jms:jme, 3 ) , &
2048 INTENT(IN ) :: erod
2049 REAL, DIMENSION( ims:ime , jms:jme ) , &
2050 INTENT(IN ) :: &
2051 u10, &
2052 v10, &
2053 xland
2054 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
2055 INTENT(IN ) :: &
2056 alt, &
2057 dz8w, &
2058 u_phy,v_phy,rho_phy
2060 REAL, INTENT(IN ) :: dt,g
2061 !
2062 ! local variables
2063 !
2064 integer :: nmx,i,j,k,imx,jmx,lmx
2065 integer,dimension (1,1) :: ilwi
2066 real*8, DIMENSION (3) :: erodin
2067 real*8, DIMENSION (5) :: tc,bems
2068 real*8, dimension (1,1) :: w10m,gwet,airden
2069 real*8 mass1part
2070 real factbb, fracoin, fracca, fracco3, fractot
2071 real*8 converi
2072 integer p1st,l_oin,l,n,l_ca, l_co3
2073 real :: densdust
2074 real totalemis
2075 integer iphase, itype
2076 real :: dz_lowest
2077 real sz(8)
2079 p1st = param_first_scalar
2081 !czhao give an example of MOSAIC 8 size bins here (unit: um) from Jerome 2005 paper
2082 ! 0.039-0.078; 0.078-0.156; 0.156-0.312; 0.312-0.625; 0.625-1.25; 1.25-2.5; 2.5-5.0; 5.0-10.0
2083 !in the model the size bound for 8 bins is calculated in the code using dlo_sec and dhi_sec
2084 !size3
2085 ! sz(1)=0.0
2086 ! sz(2)=0.0
2087 ! sz(3)=0.0003
2088 ! sz(4)=0.0048
2089 ! sz(5)=0.0130
2090 ! sz(6)=0.0250
2091 ! sz(7)=0.1400
2092 ! sz(8)=0.8169
2093 !the size is following the best match2 of modal distribution
2094 !we have ~3.6% within 1 um radius and ~7.4% within 1.25 um radius
2095 !versus ~6.0% within 1 um radius following GOCART original and ~9.3% within 1.25 um radius
2096 ! totally ~82% are within the MOSAIC size range 0.03-10 um diameter
2097 ! now the dust distribution follows the parameter defined in module_data_sorgam.F
2098 sz(1)=1.0e-8
2099 sz(2)=1.0e-6
2100 sz(3)=3.0e-4
2101 sz(4)=3.5e-3
2102 sz(5)=0.018
2103 sz(6)=0.070
2104 sz(7)=0.2595
2105 sz(8)=0.42
2106 ! frac_10um=sum(sz(1:8))
2107 ! sz(:)=sz(:)/sum(sz(1:8)) ! relative ratio in diameter range of 0-10 um
2109 ! for now just do itype=1
2110 itype = 1
2111 iphase = ai_phase
2113 ! added option for 4bin WRF
2114 IF (nsize_aer(itype) == 4) THEN
2115 sz(1) = sz(1) + sz(2)
2116 sz(2) = sz(3) + sz(4)
2117 sz(3) = sz(5) + sz(6)
2118 sz(4) = sz(7) + sz(8)
2119 ENDIF
2121 converi=1.e9
2122 !
2123 ! number of dust bins
2124 !
2125 imx=1
2126 jmx=1
2127 lmx=1
2128 nmx=ndust
2129 k=kts
2130 do j=jts,jte
2131 do i=its,ite
2133 ! no dust over water!!!
2134 !
2135 if(xland(i,j).lt.1.5)then
2136 ilwi(1,1)=1
2138 ! tc is dummy parameter. Just to conform with source_du() call from module_gocart_dust.F
2139 tc(:) = 0.0
2141 w10m=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j))
2143 ! we don't trust the u10,v10 values, if model layers are very thin near surface
2144 if(dz8w(i,kts,j).lt.12.)w10m=sqrt(u_phy(i,kts,j)*u_phy(i,kts,j)+v_phy(i,kts,j)*v_phy(i,kts,j))
2145 erodin(1)=erod(i,j,1)
2146 erodin(2)=erod(i,j,2)
2147 erodin(3)=erod(i,j,3)
2148 !
2149 ! volumetric soil moisture over porosity
2150 gwet(1,1)=smois(i,k,j)/porosity(isltyp(i,j))
2151 airden(1,1)=rho_phy(i,k,j)
2152 dz_lowest = dz8w(i,k,j)
2154 !call from module_gocart_dust.F
2155 call source_du( imx,jmx,lmx,nmx, dt, tc, &
2156 erodin, ilwi, w10m, gwet, airden, &
2157 dz_lowest,bems,start_month,g)
2159 ! bems: (kg/m2/timestep)
2161 ! A. Ukhov
2162 ! GOCART radii range for 1..5 bins: 0.1-10 um, where 4th bin: 3-6 um, 5th bin: 6-10 um.
2163 ! MOSAIC radii range for 1..8 bins: 0.0195-5 um, where 8th bin: 2.5-5 um.
2164 ! Therefore account only 1-4 GOCART bins, when calculating total dust flux from the surface.
2165 ! ln(5/3)/ln(6/3)=0.74 is contribution from the 4 GOCART bin.
2166 totalemis=(sum(bems(1:3))+bems(4)*0.74)*converi !ug/m2/timestep
2168 ! A. Ukhov
2169 ! for output diagnostics
2170 ! bems (kg/m2) per dt
2171 ! p_edust1...5 is accumulated dust emission (kg/m2)
2172 emis_dust(i,1,j,p_edust1)=emis_dust(i,1,j,p_edust1)+bems(1)
2173 emis_dust(i,1,j,p_edust2)=emis_dust(i,1,j,p_edust2)+bems(2)
2174 emis_dust(i,1,j,p_edust3)=emis_dust(i,1,j,p_edust3)+bems(3)
2175 emis_dust(i,1,j,p_edust4)=emis_dust(i,1,j,p_edust4)+bems(4)
2176 emis_dust(i,1,j,p_edust5)=emis_dust(i,1,j,p_edust5)+bems(5)
2178 factbb = (totalemis/dz8w(i,k,j))*alt(i,k,j) !(ug/kg)
2179 fracoin = 0.97
2180 fracca = 0.03*0.4
2181 fracco3 = 0.03*0.6
2182 fractot = fracoin + fracca + fracco3
2184 do 1810 n = 1, nsize_aer(itype)
2185 l_oin = lptr_oin_aer(n,itype,iphase)
2186 l_ca = lptr_ca_aer(n,itype,iphase)
2187 l_co3 = lptr_co3_aer(n,itype,iphase)
2189 if (l_oin >= p1st) then
2190 chem(i,k,j,l_oin) = chem(i,k,j,l_oin) + &
2191 factbb * sz(n) * fracoin
2192 if (l_ca >= p1st) &
2193 chem(i,k,j,l_ca) = chem(i,k,j,l_ca) + &
2194 factbb * sz(n) * fracca
2195 if (l_co3 >= p1st) &
2196 chem(i,k,j,l_co3) = chem(i,k,j,l_co3) + &
2197 factbb * sz(n) * fracco3
2199 if (n <= 5) densdust=2.5
2200 if (n > 5 ) densdust=2.65
2201 ! added option for 4bin WRF
2202 if (nsize_aer(itype) == 4) then
2203 if (n <= 2) densdust=2.5
2204 if (n > 2 ) densdust=2.65
2205 endif
2207 ! mass1part is mass of a single particle in ug
2208 mass1part=0.523598*(dcen_sect(n,itype)**3)*densdust*1.0e06
2209 l = numptr_aer(n,itype,iphase)
2210 if (l >= p1st) chem(i,k,j,l) = chem(i,k,j,l) + factbb * sz(n) * fractot/ mass1part
2211 end if
2212 1810 continue
2214 endif ! xland
2215 enddo ! i
2216 enddo ! j
2219 end subroutine mosaic_dust_gocartemis
2220 !===========================================================================
2222 !----------------------------------------------------------------------
2223 subroutine addemiss_masscheck( id, config_flags, iflagaa, fromwhere, &
2224 dtstep, efact1, dz8w, chem, chem_sum, &
2225 ids,ide, jds,jde, kds,kde, &
2226 ims,ime, jms,jme, kms,kme, &
2227 its,ite, jts,jte, kts,kte, &
2228 nemit, &
2229 e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, &
2230 e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21 )
2231 !
2232 ! produces test diagnostics for "addemiss" routines
2233 !
2234 ! 1. computes {sum over i,j,k ( chem * dz8w )} before and after
2235 ! emissions tendencies are added to chem,
2236 ! then prints (sum_after - sum_before)/(dtstep*efact1)
2237 ! 2. computes {sum over i,j,k ( e_xxx )}, then prints them
2238 ! the two should be equal
2239 !
2241 USE module_configure, only: grid_config_rec_type
2242 USE module_state_description, only: num_chem
2244 IMPLICIT NONE
2246 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2248 INTEGER, INTENT(IN ) :: id, iflagaa, &
2249 ids,ide, jds,jde, kds,kde, &
2250 ims,ime, jms,jme, kms,kme, &
2251 its,ite, jts,jte, kts,kte, &
2252 nemit
2254 REAL, INTENT(IN ) :: dtstep, efact1
2256 ! trace species mixing ratios
2257 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
2258 INTENT(IN ) :: chem
2260 ! trace species integrals
2261 DOUBLE PRECISION, DIMENSION( num_chem ), &
2262 INTENT(INOUT ) :: chem_sum
2264 ! layer thickness (m)
2265 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2266 INTENT(IN ) :: dz8w
2268 ! emissions
2269 ! REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
2270 REAL, DIMENSION( ims:ime , kms:config_flags%kemit , jms:jme ), &
2271 INTENT(IN ) :: &
2272 e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, &
2273 e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21
2275 character(len=*), intent(in) :: fromwhere
2277 ! local variables
2278 integer, parameter :: nemit_maxd = 21
2279 integer :: i, j, k, l
2280 double precision :: chem_sum_prev
2281 real :: fact
2282 real :: emit_sum(nemit_maxd)
2285 ! compute column integral, summed over i-j grids
2286 ! compute {sum over i,j,k ( chem * dz8w ) }
2287 ! on second pass (iflagaa==2), subtract the pass-one sum
2288 do 1900 l = 1, num_chem
2290 chem_sum_prev = chem_sum(l)
2291 chem_sum(l) = 0.0
2293 do j = jts, jte
2294 do k = kts, kte
2295 do i = its, ite
2296 chem_sum(l) = chem_sum(l) + dble( chem(i,k,j,l)*dz8w(i,k,j) )
2297 end do
2298 end do
2299 end do
2301 if (iflagaa == 2) chem_sum(l) = (chem_sum(l) - chem_sum_prev)
2303 1900 continue
2305 if (iflagaa /= 2) return
2308 ! compute {sum over i,j,k ( e_xxx ) }
2309 emit_sum(:) = 0.0
2311 do 2900 l = 1, min(nemit,nemit_maxd)
2312 do j = jts, jte
2313 do k = kts, min(config_flags%kemit,kte)
2314 do i = its, ite
2315 if (l== 1) emit_sum(l) = emit_sum(l) + e01(i,k,j)
2316 if (l== 2) emit_sum(l) = emit_sum(l) + e02(i,k,j)
2317 if (l== 3) emit_sum(l) = emit_sum(l) + e03(i,k,j)
2318 if (l== 4) emit_sum(l) = emit_sum(l) + e04(i,k,j)
2319 if (l== 5) emit_sum(l) = emit_sum(l) + e05(i,k,j)
2320 if (l== 6) emit_sum(l) = emit_sum(l) + e06(i,k,j)
2321 if (l== 7) emit_sum(l) = emit_sum(l) + e07(i,k,j)
2322 if (l== 8) emit_sum(l) = emit_sum(l) + e08(i,k,j)
2323 if (l== 9) emit_sum(l) = emit_sum(l) + e09(i,k,j)
2324 if (l==10) emit_sum(l) = emit_sum(l) + e10(i,k,j)
2326 if (l==11) emit_sum(l) = emit_sum(l) + e11(i,k,j)
2327 if (l==12) emit_sum(l) = emit_sum(l) + e12(i,k,j)
2328 if (l==13) emit_sum(l) = emit_sum(l) + e13(i,k,j)
2329 if (l==14) emit_sum(l) = emit_sum(l) + e14(i,k,j)
2330 if (l==15) emit_sum(l) = emit_sum(l) + e15(i,k,j)
2331 if (l==16) emit_sum(l) = emit_sum(l) + e16(i,k,j)
2332 if (l==17) emit_sum(l) = emit_sum(l) + e17(i,k,j)
2333 if (l==18) emit_sum(l) = emit_sum(l) + e18(i,k,j)
2334 if (l==19) emit_sum(l) = emit_sum(l) + e19(i,k,j)
2335 if (l==20) emit_sum(l) = emit_sum(l) + e20(i,k,j)
2337 if (l==21) emit_sum(l) = emit_sum(l) + e21(i,k,j)
2338 end do
2339 end do
2340 end do
2341 2900 continue
2343 ! output the chem_sum and emit_sum
2344 print 9110, fromwhere, its, ite, jts, jte
2345 print 9100, 'chem_sum'
2346 fact = 1.0/(dtstep*efact1)
2347 print 9120, (l, fact*chem_sum(l), l=1,num_chem)
2348 print 9100, 'emit_sum'
2349 print 9120, (l, emit_sum(l), l=1,min(nemit,nemit_maxd))
2351 9100 format( a )
2352 9110 format( / 'addemiss_masscheck output, fromwhere = ', a / &
2353 'its, ite, jts, jte =', 4i5 )
2354 9120 format( 5( i5, 1pe11.3 ) )
2357 return
2358 END subroutine addemiss_masscheck