Line data Source code
1 : module mo_setsoa
2 :
3 : use shr_kind_mod, only : r8 => shr_kind_r8
4 : use cam_logfile, only : iulog
5 : use ppgrid, only : pcols, pver, begchunk, endchunk
6 : use cam_abortutils, only : endrun
7 : use mo_constants, only : avogadro, Rgas
8 :
9 : implicit none
10 : private
11 : public :: soa_inti, setsoa, has_soa
12 : public :: soa_register
13 :
14 : save
15 :
16 : integer, parameter :: NRX = 10 ! number of SOA forming reactions
17 : integer, parameter :: NPR = 2 ! number of products (always 2!)
18 : integer, target :: spc_ndx(22)
19 : integer, pointer :: soam_ndx, soai_ndx, soab_ndx, soat_ndx, soax_ndx
20 : integer, pointer :: sogm_ndx, sogi_ndx, sogb_ndx, sogt_ndx, sogx_ndx
21 : integer, pointer :: oc1_ndx, oc2_ndx, c10h16_ndx, isop_ndx, o3_ndx, oh_ndx
22 : integer, pointer :: no3_ndx, no_ndx, ho2_ndx, tolo2_ndx, beno2_ndx, xylo2_ndx
23 :
24 : integer :: rxn_soa(nrx),react_ndx(NRX,NPR)
25 : real(r8), dimension(NRX,NPR) :: alpha ! mass-based stoichiometric coefficients
26 : real(r8), dimension(NRX,NPR) :: k_om ! equilibrium gas-particule partition
27 : real(r8), dimension(NRX) :: T1, delH ! Clausium Clayperson parameters (K, J/mol)
28 : real(r8), dimension(nrx,npr) :: fracsog_init ! mass fraction of each SOA class from each reaction
29 : real(r8), dimension(nrx,npr) :: fracsoa_init ! mass fraction of each SOG class from each reaction
30 :
31 : integer :: fracsog_ndx = -1
32 : integer :: fracsoa_ndx = -1
33 :
34 : integer, pointer :: soa_ndx
35 : integer, pointer :: bigalk_ndx, toluene_ndx
36 :
37 : real(r8), dimension(6) :: bulk_yield ! total yield of condensable compound (ug/m3/ppm)
38 : real(r8), dimension(6) :: fraction ! fraction of VOC used in reaction
39 :
40 : real(r8), parameter :: OMscale = 2.1_r8 ! scaling factor for OM:OC [Turpin and Lim, 2001]
41 : logical :: has_soa = .false.
42 : logical :: has_soa_equil = .false.
43 :
44 : contains
45 :
46 : !===============================================================================
47 : !===============================================================================
48 1536 : subroutine soa_inti(pbuf2d)
49 : use physics_buffer, only : physics_buffer_desc
50 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
51 1536 : if ( has_soa_equil) then
52 0 : call soa_inti_equil(pbuf2d)
53 : else
54 1536 : call soa_inti_old()
55 : endif
56 1536 : endsubroutine soa_inti
57 : !===============================================================================
58 : !===============================================================================
59 0 : subroutine setsoa( ncol, lchnk, dt, reaction_rates, tfld, xhnm, vmr, pbuf)
60 1536 : use physics_buffer, only : physics_buffer_desc
61 : use chem_mods, only : gas_pcnst, rxntot
62 : !-----------------------------------------------------------------------
63 : ! ... dummy arguments
64 : !-----------------------------------------------------------------------
65 : integer, intent(in) :: ncol ! number columns in chunkx
66 : integer, intent(in) :: lchnk ! chunk index
67 : real(r8), intent(in) :: dt ! time step
68 : real(r8), intent(in) :: reaction_rates(:,:,:) ! reaction rates
69 : real(r8), intent(in) :: tfld(:,:) ! temperature (K)
70 : real(r8), intent(in) :: xhnm(:,:) ! total atms density (molec/cm**3)
71 : real(r8), intent(inout) :: vmr(:,:,:) ! xported species ( vmr )
72 : type(physics_buffer_desc), pointer :: pbuf(:)
73 :
74 0 : if ( has_soa_equil) then
75 0 : call setsoa_equil(dt,reaction_rates,tfld,vmr,xhnm,ncol,lchnk,pbuf)
76 : else
77 0 : call setsoa_old(dt,reaction_rates,tfld,vmr,xhnm,ncol,lchnk)
78 : endif
79 :
80 0 : end subroutine setsoa
81 : !===============================================================================
82 : !===============================================================================
83 1536 : subroutine soa_inti_old
84 :
85 0 : use mo_chem_utls, only : get_spc_ndx, get_rxt_ndx
86 : use cam_history, only : addfld
87 : use ppgrid, only : pver
88 : use spmd_utils, only : masterproc
89 :
90 : implicit none
91 :
92 : !-----------------------------------------------------------------------
93 : ! ... check if this is an aerosol simulation
94 : !-----------------------------------------------------------------------
95 1536 : if( .not. has_soa ) then
96 : return
97 : end if
98 :
99 : !-----------------------------------------------------------------------
100 : ! ... set reaction indicies
101 : !-----------------------------------------------------------------------
102 0 : rxn_soa(1) = get_rxt_ndx( 'soa1' )
103 0 : if ( rxn_soa(1) <= 0 ) then
104 0 : rxn_soa(1) = get_rxt_ndx( 'C10H16_O3' )
105 : end if
106 0 : rxn_soa(2) = get_rxt_ndx( 'soa2' )
107 0 : if ( rxn_soa(2) <= 0 ) then
108 0 : rxn_soa(2) = get_rxt_ndx( 'C10H16_OH' )
109 : end if
110 0 : rxn_soa(3) = get_rxt_ndx( 'soa3' )
111 0 : if ( rxn_soa(3) <= 0 ) then
112 0 : rxn_soa(3) = get_rxt_ndx( 'C10H16_NO3' )
113 : end if
114 0 : rxn_soa(4) = get_rxt_ndx( 'soa4' )
115 0 : if ( rxn_soa(4) <= 0 ) then
116 0 : rxn_soa(4) = get_rxt_ndx( 'TOLUENE_OH' )
117 : end if
118 0 : rxn_soa(5) = get_rxt_ndx( 'soa4' ) ! TOLUENE is a lumped species and there are two sets of pathways
119 0 : if ( rxn_soa(5) <= 0 ) then
120 0 : rxn_soa(5) = get_rxt_ndx( 'TOLUENE_OH' )
121 : end if
122 0 : rxn_soa(6) = get_rxt_ndx( 'soa5' )
123 0 : if ( rxn_soa(6) <= 0 ) then
124 0 : rxn_soa(6) = get_rxt_ndx( 'BIGALK_OH' )
125 : end if
126 0 : if( all( rxn_soa(:) < 1 ) ) then
127 0 : has_soa = .false.
128 0 : return
129 : else
130 0 : if (masterproc) then
131 0 : write(iulog,*) '-----------------------------------------'
132 0 : write(iulog,*) 'soa_inti_old: mozart will do soa aerosols'
133 0 : write(iulog,*) '-----------------------------------------'
134 : endif
135 : end if
136 :
137 : !
138 : ! define reactants
139 : !
140 0 : react_ndx(1,1) = c10h16_ndx
141 0 : react_ndx(1,2) = o3_ndx
142 0 : react_ndx(2,1) = c10h16_ndx
143 0 : react_ndx(2,2) = oh_ndx
144 0 : react_ndx(3,1) = c10h16_ndx
145 0 : react_ndx(3,2) = no3_ndx
146 0 : react_ndx(4,1) = toluene_ndx
147 0 : react_ndx(4,2) = oh_ndx
148 0 : react_ndx(5,1) = toluene_ndx
149 0 : react_ndx(5,2) = oh_ndx
150 0 : react_ndx(6,1) = bigalk_ndx
151 0 : react_ndx(6,2) = oh_ndx
152 :
153 0 : if ( masterproc ) then
154 0 : write(iulog,*)'soa_inti ',c10h16_ndx, o3_ndx, oh_ndx, no3_ndx, bigalk_ndx, toluene_ndx
155 0 : write(iulog,*)'soa_inti ',soa_ndx, oc1_ndx, oc2_ndx
156 0 : write(iulog,*)'soa_inti ',react_ndx
157 : endif
158 : !
159 : ! define partitioning coefficients for each reaction
160 : ! bulk yields are from Seinfeld and Pandis (1998)
161 : !
162 : ! c10h16 + o3 (from Chung and Seinfeld, JGR, 107, 2002)
163 : !
164 0 : alpha(1,1) = 0.067_r8
165 0 : alpha(1,2) = 0.354_r8
166 0 : k_om (1,1) = 0.184_r8
167 0 : k_om (1,2) = 0.0043_r8
168 0 : fraction(1) = 1._r8
169 0 : bulk_yield(1) = 762._r8
170 : !
171 : ! c10h16 + oh (from Chung and Seinfeld, JGR, 107, 2002)
172 : !
173 0 : alpha(2,1) = 0.067_r8
174 0 : alpha(2,2) = 0.354_r8
175 0 : k_om (2,1) = 0.184_r8
176 0 : k_om (2,2) = 0.0043_r8
177 0 : fraction(2) = 1._r8
178 0 : bulk_yield(2) = 762._r8
179 : !
180 : ! c10h16 + no3 (from Chung and Seinfeld, JGR, 107, 2002)
181 : !
182 0 : alpha(3,1) = 1.000_r8
183 0 : alpha(3,2) = 0.000_r8
184 0 : k_om (3,1) = 0.0163_r8
185 0 : k_om (3,2) = 0.0000_r8
186 0 : fraction(3) = 1._r8
187 0 : bulk_yield(3) = 762._r8
188 : !
189 : ! toluene + oh : toluene (from Odum et al., Environ. Sci. Technol., 1892, 1997)
190 : !
191 0 : alpha(4,1) = 0.038_r8
192 0 : alpha(4,2) = 0.167_r8
193 0 : k_om (4,1) = 0.042_r8
194 0 : k_om (4,2) = 0.0014_r8
195 0 : fraction(4) = 0.7_r8
196 0 : bulk_yield(4) = 424._r8
197 : !
198 : ! toluene + oh : m-xylene (from Cocker et al., Atmos. Environ., 6079, 2001)
199 : !
200 0 : alpha(5,1) = 0.120_r8
201 0 : alpha(5,2) = 0.019_r8
202 0 : k_om (5,1) = 0.060_r8
203 0 : k_om (5,2) = 0.010_r8
204 0 : fraction(5) = 0.2_r8
205 0 : bulk_yield(5) = 419._r8
206 : !
207 : ! bigalk + oh : only for alkanes >= heptane (assume low-yield aromatics as in Lack et al.)
208 : ! (from Odum et al., Environ. Sci. Technol., 1892, 1997)
209 : !
210 0 : alpha(6,1) = 0.071_r8
211 0 : alpha(6,2) = 0.138_r8
212 0 : k_om (6,1) = 0.053_r8
213 0 : k_om (6,2) = 0.0019_r8
214 0 : fraction(6) = 0.1_r8
215 0 : bulk_yield(6) = 200._r8
216 : !
217 0 : call addfld( 'SOA_PROD', (/ 'lev' /), 'A', 'kg/kg/s', 'production of SOA' )
218 :
219 0 : return
220 1536 : end subroutine soa_inti_old
221 : !
222 0 : subroutine setsoa_old(dt,reaction_rates,tfld,vmr,xhnm,ncol,lchnk)
223 : !
224 : ! secondary organic aerosol for mozart v2.5
225 : !
226 : ! based on Lack et al., JGR, 109, D03203, 2004
227 : !
228 : ! rewritten by Jean-Francois Lamarque for updated chemical
229 : ! mechanism (March 2004)
230 : !
231 : ! adapted to CAM (May 2004)
232 : !
233 1536 : use ppgrid, only : pcols, pver
234 : use chem_mods, only : adv_mass, gas_pcnst, rxntot
235 : use mo_chem_utls, only : get_spc_ndx, get_rxt_ndx
236 : use cam_history, only : outfld
237 : use cam_abortutils, only : endrun
238 : !
239 : implicit none
240 : !
241 : !-----------------------------------------------------------------------
242 : ! ... dummy arguments
243 : !-----------------------------------------------------------------------
244 : integer, intent(in) :: ncol ! number columns in chunkx
245 : integer, intent(in) :: lchnk ! chunk index
246 : real(r8), intent(in) :: dt ! time step
247 : real(r8), intent(in) :: reaction_rates(:,:,:) ! reaction rates
248 : real(r8), intent(inout) :: vmr(:,:,:) ! xported species ( vmr )
249 : real(r8), intent(in) :: tfld(:,:) ! temperature (K)
250 : real(r8), intent(in) :: xhnm(:,:) ! total atms density (mol/cm**3)
251 :
252 : !-----------------------------------------------------------------------
253 : ! ... local variables
254 : !-----------------------------------------------------------------------
255 : integer :: i,k,n
256 : real(r8) :: m_0
257 : real(r8) :: mw_soa,yield,prod,soa_mass
258 0 : real(r8) :: soa_prod(ncol,pver)
259 : !
260 : ! find molecular weight of SOA
261 : !
262 0 : mw_soa = adv_mass(soa_ndx)
263 : !
264 : do k=1,pver
265 : do i=1,ncol
266 : !
267 : ! calculate initial mass of organic aerosols from OC1 and OC2
268 : ! and convert to ug/m3
269 : !
270 : m_0 = (vmr(i,k,oc1_ndx)+vmr(i,k,oc2_ndx)) * xhnm(i,k) * adv_mass(oc1_ndx)/avogadro * 1.e12_r8
271 : !
272 : ! switch based on a minimum value of m_0. The bulk approach is
273 : ! used to initiate the process
274 : !
275 : if ( m_0 <= 0.2_r8 ) then
276 : !
277 : ! bulk theory
278 : !
279 : soa_mass = 0._r8
280 : do n=1,6
281 : !
282 : if ( rxn_soa(n) <= 0 ) cycle
283 : !
284 : yield = bulk_yield(n)
285 : !
286 : ! define chemical production from gas-phase chemistry
287 : !
288 : prod = reaction_rates(i,k,rxn_soa(n)) * fraction(n) &
289 : * vmr(i,k,react_ndx(n,1)) * vmr(i,k,react_ndx(n,2)) * dt
290 : !
291 : ! convert from mixing ratio to ppm
292 : !
293 : prod = 1e6_r8 * prod
294 : !
295 : ! collect into total SOA mass
296 : !
297 : soa_mass = soa_mass + yield * prod
298 : !
299 : end do
300 : !
301 : else
302 : !
303 : ! partitioning theory
304 : !
305 : soa_mass = 0._r8
306 : do n=1,6
307 : !
308 : if ( rxn_soa(n) <= 0 ) cycle
309 : !
310 : ! define yield from available m_0
311 : !
312 : yield = soa_yield(m_0,alpha(n,1:2),k_om(n,1:2))
313 : !
314 : ! define chemical production from gas-phase chemistry
315 : !
316 : prod = reaction_rates(i,k,rxn_soa(n)) * fraction(n) &
317 : * vmr(i,k,react_ndx(n,1)) * vmr(i,k,react_ndx(n,2)) * dt
318 : !
319 : ! convert from mixing ratio to mass (ug/m3)
320 : !
321 : prod = prod * xhnm(i,k) * mw_soa/avogadro * 1.e12_r8
322 : !
323 : ! collect into total SOA mass
324 : !
325 : soa_mass = soa_mass + yield * prod
326 : !
327 : end do
328 : !
329 : endif
330 : !
331 : ! convert from ug/m3 to mixing ratio and update vmr
332 : !
333 : vmr(i,k,soa_ndx) = vmr(i,k,soa_ndx) + soa_mass * 1.e-12_r8 * avogadro/(mw_soa*xhnm(i,k))
334 : if ( vmr(i,k,soa_ndx) > 1.e0_r8 ) then
335 : write(iulog,*)i,k,soa_mass,m_0
336 : call endrun('soa_yield: vmr(i,k,soa_ndx) > 1.e0_r8')
337 : endif
338 : !
339 : soa_prod(i,k) = soa_mass*1.e-12_r8*avogadro/(28.966_r8*xhnm(i,k)*dt)
340 : end do
341 : end do
342 : !
343 : call outfld('SOA_PROD',soa_prod(:ncol,:),ncol, lchnk)
344 : return
345 : end subroutine setsoa_old
346 : !
347 : real(r8) function soa_yield(m_0,xalpha,xk)
348 : !
349 : implicit none
350 : !
351 : real(r8), intent(in) :: m_0
352 : real(r8), intent(in), dimension(2) :: xalpha, xk
353 : !
354 : soa_yield = m_0 * ( ((xalpha(1)*xk(1))/(1._r8+xk(1)*m_0)) &
355 : + ((xalpha(2)*xk(2))/(1._r8+xk(2)*m_0)) )
356 : !
357 : return
358 0 : end function soa_yield
359 : !
360 :
361 : !===============================================================================
362 : !===============================================================================
363 1536 : subroutine soa_register
364 : use physics_buffer, only : pbuf_add_field, dtype_r8
365 : use mo_chem_utls, only : get_spc_ndx
366 :
367 : !-----------------------------------------------------------------------
368 : ! ... set species indices
369 : !-----------------------------------------------------------------------
370 :
371 1536 : oc1_ndx => spc_ndx(1)
372 1536 : oc2_ndx => spc_ndx(2)
373 1536 : soam_ndx => spc_ndx(3)
374 1536 : soai_ndx => spc_ndx(4)
375 1536 : soat_ndx => spc_ndx(5)
376 1536 : soab_ndx => spc_ndx(6)
377 1536 : soax_ndx => spc_ndx(7)
378 1536 : c10h16_ndx => spc_ndx(8)
379 1536 : isop_ndx => spc_ndx(9)
380 1536 : tolo2_ndx => spc_ndx(10)
381 1536 : beno2_ndx => spc_ndx(11)
382 1536 : xylo2_ndx => spc_ndx(12)
383 1536 : ho2_ndx => spc_ndx(13)
384 1536 : no_ndx => spc_ndx(14)
385 1536 : o3_ndx => spc_ndx(15)
386 1536 : oh_ndx => spc_ndx(16)
387 1536 : no3_ndx => spc_ndx(17)
388 1536 : sogm_ndx => spc_ndx(18)
389 1536 : sogi_ndx => spc_ndx(19)
390 1536 : sogt_ndx => spc_ndx(20)
391 1536 : sogb_ndx => spc_ndx(21)
392 1536 : sogx_ndx => spc_ndx(22)
393 :
394 1536 : oc1_ndx = get_spc_ndx( 'OC1' )
395 1536 : oc2_ndx = get_spc_ndx( 'OC2' )
396 1536 : c10h16_ndx = get_spc_ndx( 'C10H16')
397 1536 : isop_ndx = get_spc_ndx( 'ISOP' )
398 1536 : tolo2_ndx = get_spc_ndx( 'TOLO2' )
399 1536 : beno2_ndx = get_spc_ndx( 'BENO2' )
400 1536 : xylo2_ndx = get_spc_ndx( 'XYLO2' )
401 1536 : ho2_ndx = get_spc_ndx( 'HO2' )
402 1536 : no_ndx = get_spc_ndx( 'NO' )
403 1536 : o3_ndx = get_spc_ndx( 'OX' )
404 1536 : if( o3_ndx < 1 ) then
405 1536 : o3_ndx = get_spc_ndx( 'O3' )
406 : end if
407 1536 : oh_ndx = get_spc_ndx( 'OH' )
408 1536 : no3_ndx = get_spc_ndx( 'NO3' )
409 :
410 1536 : soam_ndx = get_spc_ndx( 'SOAM' )
411 1536 : soai_ndx = get_spc_ndx( 'SOAI' )
412 1536 : soat_ndx = get_spc_ndx( 'SOAT' )
413 1536 : soab_ndx = get_spc_ndx( 'SOAB' )
414 1536 : soax_ndx = get_spc_ndx( 'SOAX' )
415 :
416 1536 : sogm_ndx = get_spc_ndx( 'SOGM' )
417 1536 : sogi_ndx = get_spc_ndx( 'SOGI' )
418 1536 : sogt_ndx = get_spc_ndx( 'SOGT' )
419 1536 : sogb_ndx = get_spc_ndx( 'SOGB' )
420 1536 : sogx_ndx = get_spc_ndx( 'SOGX' )
421 :
422 1536 : has_soa_equil = all( spc_ndx(1:22) > 0 )
423 1536 : has_soa = has_soa_equil
424 :
425 1536 : if ( has_soa_equil ) then
426 : ! fracsog and fracsoa are added to pbuffer for persistence across restarts
427 0 : call pbuf_add_field( 'FRACSOG' ,'global',dtype_r8,(/pcols,pver,nrx,npr/), fracsog_ndx)
428 0 : call pbuf_add_field( 'FRACSOA' ,'global',dtype_r8,(/pcols,pver,nrx,npr/), fracsoa_ndx)
429 : else
430 : ! reassign these ndx pointers...
431 1536 : soa_ndx => spc_ndx(1)
432 1536 : oc1_ndx => spc_ndx(2)
433 1536 : oc2_ndx => spc_ndx(3)
434 1536 : c10h16_ndx => spc_ndx(4)
435 1536 : o3_ndx => spc_ndx(5)
436 1536 : oh_ndx => spc_ndx(6)
437 1536 : no3_ndx => spc_ndx(7)
438 1536 : bigalk_ndx => spc_ndx(8)
439 1536 : toluene_ndx => spc_ndx(9)
440 :
441 1536 : soa_ndx = get_spc_ndx( 'SOA' )
442 1536 : oc1_ndx = get_spc_ndx( 'OC1' )
443 1536 : oc2_ndx = get_spc_ndx( 'OC2' )
444 1536 : c10h16_ndx = get_spc_ndx( 'C10H16')
445 1536 : o3_ndx = get_spc_ndx( 'OX' )
446 1536 : if( o3_ndx < 1 ) then
447 1536 : o3_ndx = get_spc_ndx( 'O3' )
448 : end if
449 1536 : oh_ndx = get_spc_ndx( 'OH' )
450 1536 : no3_ndx = get_spc_ndx( 'NO3' )
451 :
452 1536 : bigalk_ndx = get_spc_ndx( 'BIGALK' )
453 1536 : toluene_ndx = get_spc_ndx( 'TOLUENE' )
454 :
455 1536 : has_soa = all( spc_ndx(1:7) > 0 )
456 : endif
457 :
458 1536 : end subroutine soa_register
459 :
460 : !===============================================================================
461 : !===============================================================================
462 0 : subroutine soa_inti_equil(pbuf2d)
463 :
464 1536 : use mo_chem_utls, only : get_spc_ndx, get_rxt_ndx
465 : use cam_history, only : addfld
466 : use ppgrid, only : pver
467 : use spmd_utils, only : masterproc
468 : use cam_control_mod, only: initial_run
469 : use physics_buffer, only : physics_buffer_desc, pbuf_set_field, pbuf_get_chunk, pbuf_get_field
470 :
471 : implicit none
472 :
473 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
474 0 : type(physics_buffer_desc), pointer :: pbuf_ptr(:)
475 :
476 : integer :: astat
477 : integer :: i,k,j, c
478 :
479 0 : real(r8), pointer :: fracsog(:,:,:,:) ! mass fraction of each SOA class from each reaction
480 0 : real(r8), pointer :: fracsoa(:,:,:,:) ! mass fraction of each SOG class from each reaction
481 :
482 : !-----------------------------------------------------------------------
483 : ! ... check if this is an aerosol simulation
484 : !-----------------------------------------------------------------------
485 0 : if( .not. has_soa ) then
486 : return
487 : end if
488 :
489 : !-----------------------------------------------------------------------
490 : ! ... set reaction indicies
491 : !-----------------------------------------------------------------------
492 0 : rxn_soa(1) = get_rxt_ndx( 'C10H16_O3' )
493 0 : rxn_soa(2) = get_rxt_ndx( 'C10H16_OH' )
494 0 : rxn_soa(3) = get_rxt_ndx( 'C10H16_NO3' )
495 0 : rxn_soa(4) = get_rxt_ndx( 'ISOP_OH' )
496 0 : rxn_soa(5) = get_rxt_ndx( 'TOLO2_HO2' )
497 0 : rxn_soa(6) = get_rxt_ndx( 'TOLO2_NO' )
498 0 : if (rxn_soa(6)<0) then
499 0 : rxn_soa(6) = get_rxt_ndx( 'ox_p12' )
500 : endif
501 0 : rxn_soa(7) = get_rxt_ndx( 'BENO2_HO2' )
502 0 : rxn_soa(8) = get_rxt_ndx( 'BENO2_NO' )
503 0 : rxn_soa(9) = get_rxt_ndx( 'XYLO2_HO2' )
504 0 : rxn_soa(10) = get_rxt_ndx( 'XYLO2_NO' )
505 0 : if( all( rxn_soa(:) < 1 ) ) then
506 0 : has_soa = .false.
507 0 : return
508 : else
509 0 : if (masterproc) then
510 0 : write(iulog,*) '-----------------------------------------'
511 0 : write(iulog,*) 'soa_inti_equil: mozart will do soa aerosols'
512 0 : write(iulog,*) '-----------------------------------------'
513 : endif
514 : end if
515 :
516 : !
517 : ! define reactants
518 : !
519 0 : react_ndx(1,1) = c10h16_ndx
520 0 : react_ndx(1,2) = o3_ndx
521 0 : react_ndx(2,1) = c10h16_ndx
522 0 : react_ndx(2,2) = oh_ndx
523 0 : react_ndx(3,1) = c10h16_ndx
524 0 : react_ndx(3,2) = no3_ndx
525 0 : react_ndx(4,1) = isop_ndx
526 0 : react_ndx(4,2) = oh_ndx
527 0 : react_ndx(5,1) = tolo2_ndx
528 0 : react_ndx(5,2) = ho2_ndx
529 0 : react_ndx(6,1) = tolo2_ndx
530 0 : react_ndx(6,2) = no_ndx
531 0 : react_ndx(7,1) = beno2_ndx
532 0 : react_ndx(7,2) = ho2_ndx
533 0 : react_ndx(8,1) = beno2_ndx
534 0 : react_ndx(8,2) = no_ndx
535 0 : react_ndx(9,1) = xylo2_ndx
536 0 : react_ndx(9,2) = ho2_ndx
537 0 : react_ndx(10,1) = xylo2_ndx
538 0 : react_ndx(10,2) = no_ndx
539 :
540 0 : if ( masterproc ) then
541 0 : print *,'soa_inti ',c10h16_ndx, isop_ndx, tolo2_ndx, beno2_ndx, xylo2_ndx,o3_ndx, oh_ndx, no3_ndx
542 0 : print *,'soa_inti ',soam_ndx, soai_ndx, soab_ndx, soat_ndx, soax_ndx, oc1_ndx, oc2_ndx
543 0 : print *,'soa_inti ',sogm_ndx, sogi_ndx, sogb_ndx, sogt_ndx, sogx_ndx
544 0 : print *,'soa_inti ',react_ndx
545 : endif
546 : !
547 : ! define partitioning coefficients for each reaction
548 : ! bulk yields are from Seinfeld and Pandis (1998)
549 : !
550 : ! c10h16 + o3 (from Chung and Seinfeld, JGR, 107, 2002)
551 : !
552 0 : alpha(1,1) = 0.067_r8
553 0 : alpha(1,2) = 0.354_r8
554 0 : k_om(1,1) = 0.184_r8
555 0 : k_om(1,2) = 0.0043_r8
556 0 : T1(1) = 310._r8
557 0 : delH(1) = 42.e3_r8
558 :
559 : ! c10h16 + oh (from Chung and Seinfeld, JGR, 107, 2002)
560 : !
561 0 : alpha(2,1) = 0.067_r8
562 0 : alpha(2,2) = 0.354_r8
563 0 : k_om(2,1) = 0.184_r8
564 0 : k_om(2,2) = 0.0043_r8
565 0 : T1(2) = 310._r8
566 0 : delH(2) = 42.e3_r8
567 : !
568 : ! c10h16 + no3 (from Chung and Seinfeld, JGR, 107, 2002)
569 : !
570 0 : alpha(3,1) = 1.000_r8
571 0 : alpha(3,2) = 0.000_r8
572 0 : k_om(3,1) = 0.0163_r8
573 0 : k_om(3,2) = 0.0000_r8
574 0 : T1(3) = 310._r8
575 0 : delH(3) = 42.e3_r8
576 : !
577 : !! isop + oh (from Henze and Seinfeld, GRL, 2006): low NOx
578 : !!
579 : ! alpha(4,1) = 0.232_r8
580 : ! alpha(4,2) = 0.0288_r8
581 : ! k_om(4,1) = 0.00862_r8
582 : ! k_om(4,2) = 1.62_r8
583 : ! T1(4) = 295._r8
584 : ! delH(4) = 42.e3_r8
585 : !!
586 : ! isop + oh (from Henze and Seinfeld, GRL, 2006): high NOx
587 : !
588 0 : alpha(4,1) = 0.264_r8
589 0 : alpha(4,2) = 0.0173_r8
590 0 : k_om(4,1) = 0.00115_r8
591 0 : k_om(4,2) = 1.52_r8
592 0 : T1(4) = 295._r8
593 0 : delH(4) = 42.e3_r8
594 : !
595 : ! toluene + oh (pers comm Seinfeld and Henze): low NOx (TOLO2 + HO2)
596 : !
597 0 : alpha(5,1) = 0.2349_r8
598 0 : alpha(5,2) = 0.0_r8
599 0 : k_om(5,1) = 1000.0_r8
600 0 : k_om(5,2) = 0.0_r8
601 0 : T1(5) = 295._r8
602 0 : delH(5) = 42.e3_r8
603 : !
604 : ! toluene + oh (pers comm Seinfeld and Henze): high NOx (TOLO2 + NO)
605 : !
606 0 : alpha(6,1) = 0.0378_r8
607 0 : alpha(6,2) = 0.0737_r8
608 0 : k_om(6,1) = 0.4300_r8
609 0 : k_om(6,2) = 0.0470_r8
610 0 : T1(6) = 295._r8
611 0 : delH(6) = 42.e3_r8
612 : !
613 : ! benzene + oh (pers comm Seinfeld and Henze): low NOx (BENO2 + HO2)
614 : !
615 0 : alpha(7,1) = 0.2272_r8
616 0 : alpha(7,2) = 0.0_r8
617 0 : k_om (7,1) = 1000.0_r8
618 0 : k_om (7,2) = 0.0_r8
619 0 : T1(7) = 295._r8
620 0 : delH(7) = 42.e3_r8
621 : !
622 : ! benzene + oh (pers comm Seinfeld and Henze): high NOx (BENO2 + NO)
623 : !
624 0 : alpha(8,1) = 0.0442_r8
625 0 : alpha(8,2) = 0.5454_r8
626 0 : k_om(8,1) = 3.3150_r8
627 0 : k_om(8,2) = 0.0090_r8
628 0 : T1(8) = 295._r8
629 0 : delH(8) = 42.e3_r8
630 : !
631 : ! xylene + oh (pers comm Seinfeld and Henze): low NOx (XYLO2 + HO2)
632 : !
633 0 : alpha(9,1) = 0.2052_r8
634 0 : alpha(9,2) = 0.0_r8
635 0 : k_om(9,1) = 1000.0_r8
636 0 : k_om(9,2) = 0.0_r8
637 0 : T1(9) = 295._r8
638 0 : delH(9) = 42.e3_r8
639 : !
640 : ! xylene + oh (pers comm Seinfeld and Henze): high NOx (XYLO2 + NO)
641 : !
642 0 : alpha(10,1) = 0.0212_r8
643 0 : alpha(10,2) = 0.0615_r8
644 0 : k_om(10,1) = 0.7610_r8
645 0 : k_om(10,2) = 0.0290_r8
646 0 : T1(10) = 295._r8
647 0 : delH(10) = 42.e3_r8
648 : !
649 0 : call addfld( 'SOAM_PROD', (/ 'lev' /), 'A', 'molec/molec/s', 'production of SOAM' )
650 0 : call addfld( 'SOAI_PROD', (/ 'lev' /), 'A', 'molec/molec/s', 'production of SOAI' )
651 0 : call addfld( 'SOAT_PROD', (/ 'lev' /), 'A', 'molec/molec/s', 'production of SOAT' )
652 0 : call addfld( 'SOAB_PROD', (/ 'lev' /), 'A', 'molec/molec/s', 'production of SOAB' )
653 0 : call addfld( 'SOAX_PROD', (/ 'lev' /), 'A', 'molec/molec/s', 'production of SOAX' )
654 :
655 0 : call addfld( 'SOAM_dens', (/ 'lev' /), 'A', 'ug/m3', 'density of SOAM' )
656 0 : call addfld( 'SOAI_dens', (/ 'lev' /), 'A', 'ug/m3', 'density of SOAI' )
657 0 : call addfld( 'SOAT_dens', (/ 'lev' /), 'A', 'ug/m3', 'density of SOAT' )
658 0 : call addfld( 'SOAB_dens', (/ 'lev' /), 'A', 'ug/m3', 'density of SOAB' )
659 0 : call addfld( 'SOAX_dens', (/ 'lev' /), 'A', 'ug/m3', 'density of SOAX' )
660 :
661 : !
662 : !initialize fracsoa for first timestep and store values for future
663 0 : fracsoa_init(1:3,:)=0.2_r8
664 0 : fracsoa_init(3,2)=0._r8
665 0 : fracsoa_init(4,:)=0.5_r8
666 0 : fracsoa_init(5:10,:)=0.33_r8
667 0 : fracsoa_init(5,2)=0._r8
668 0 : fracsoa_init(7,2)=0._r8
669 0 : fracsoa_init(9,2)=0._r8
670 :
671 : !initialize fracsog for first timestep and store values for future
672 0 : fracsog_init(1:3,:)=0.2_r8
673 0 : fracsog_init(3,2)=0._r8
674 0 : fracsog_init(4,:)=0.5_r8
675 0 : fracsog_init(5:10,:)=0.33_r8
676 0 : fracsog_init(5,2)=0._r8
677 0 : fracsog_init(7,2)=0._r8
678 0 : fracsog_init(9,2)=0._r8
679 :
680 0 : if (initial_run) then
681 0 : do c=begchunk, endchunk
682 0 : pbuf_ptr=>pbuf_get_chunk(pbuf2d, c)
683 0 : call pbuf_get_field(pbuf_ptr, fracsoa_ndx, fracsoa )
684 0 : call pbuf_get_field(pbuf_ptr, fracsog_ndx, fracsog )
685 0 : do i = 1,pcols
686 0 : do k = 1,pver
687 0 : fracsoa( i,k, :,: ) = fracsoa_init(:,:)
688 0 : fracsog( i,k, :,: ) = fracsog_init(:,:)
689 : enddo
690 : enddo
691 : enddo
692 : endif
693 :
694 : return
695 0 : end subroutine soa_inti_equil
696 : !===============================================================================
697 : ! clh (08/01/06): added T dependence of gas-aerosol phase partitioning
698 : ! added isoprene as SOA precursor [Henze and Seinfeld, 2006]
699 : ! added anthropogenic SOA according to [Henze et al., 2007]
700 : ! split SOA into classes (SOAM=monoterpenes,SOAI=isoprene,SOAB=benzene,SOAT=toluene
701 : ! SOAX=xylene)
702 : ! fixed error in yield calculation
703 : ! added scaling for OC to OM in pre-existing aerosol mass
704 : ! modified yield calculation to allow for re-evaporation
705 : ! Note: do not need to subtract formed aerosol product from reactants, because gas-phase
706 : ! products do not account for entire mass
707 : !===============================================================================
708 0 : subroutine setsoa_equil(dt,reaction_rates,tfld,vmr,xhnm,ncol,lchnk,pbuf)
709 : !
710 : ! updated SOA mechanism
711 : ! based on Chung and Seinfeld, JGR, 2002
712 : !
713 : ! implemented in CAM by Colette Heald (summer 2007)
714 : !
715 0 : use ppgrid, only : pcols, pver
716 : use chem_mods, only : adv_mass, gas_pcnst, rxntot
717 : use cam_history, only : outfld
718 : use physics_buffer, only : physics_buffer_desc, pbuf_get_field
719 : !
720 : implicit none
721 : !
722 : !-----------------------------------------------------------------------
723 : ! ... dummy arguments
724 : !-----------------------------------------------------------------------
725 : integer, intent(in) :: ncol ! number columns in chunkx
726 : integer, intent(in) :: lchnk ! chunk index
727 : real(r8), intent(in) :: dt ! time step
728 : real(r8), intent(in) :: reaction_rates(:,:,:) ! reaction rates
729 : real(r8), intent(inout) :: vmr(:,:,:) ! xported species ( vmr )
730 : real(r8), intent(in) :: tfld(:,:) ! temperature (K)
731 : real(r8), intent(in) :: xhnm(:,:) ! total atms density (molec/cm**3)
732 : type(physics_buffer_desc), pointer :: pbuf(:)
733 :
734 : !-----------------------------------------------------------------------
735 : ! ... local variables
736 : !-----------------------------------------------------------------------
737 : integer :: i,k,n,p,iter
738 : real(r8) :: T_fac
739 : real(r8) :: delHC, sumorg, poa, mnew, numer, maxM, minM, tol, mw_soa,m_air
740 : real(r8) :: k_om_T(NRX,NPR) ! equilibrium gas-particule partition (T dependent)
741 : real(r8) :: prod(NRX,NPR) ! oxidized produced (alpha_i*delHC)
742 : real(r8) :: sog0(NRX,NPR) ! pre-existing SOG
743 : real(r8) :: soa0(NRX,NPR) ! pre-exisiting SOA
744 : real(r8) :: orggas(NRX,NPR) ! intermediate (SOG0+prod of ox)
745 : real(r8) :: sog(NRX,NPR) ! final SOG
746 : real(r8) :: soa(NRX,NPR) ! final SOA
747 : real(r8), dimension(pcols,pver) :: soam_mass,soai_mass,soat_mass,soab_mass,soax_mass
748 : real(r8), dimension(pcols,pver) :: sogm_mass,sogi_mass,sogt_mass,sogb_mass,sogx_mass
749 0 : real(r8) :: soam_prod(ncol,pver), soai_prod(ncol,pver),soat_prod(ncol,pver),soab_prod(ncol,pver),soax_prod(ncol,pver)
750 :
751 0 : real(r8), pointer, dimension(:,:,:,:) :: fracsog ! mass fraction of each SOA class from each reaction
752 0 : real(r8), pointer, dimension(:,:,:,:) :: fracsoa ! mass fraction of each SOG class from each reaction
753 :
754 0 : soam_mass(:,:)=0._r8
755 0 : soai_mass(:,:)=0._r8
756 0 : soat_mass(:,:)=0._r8
757 0 : soab_mass(:,:)=0._r8
758 0 : soax_mass(:,:)=0._r8
759 0 : sogm_mass(:,:)=0._r8
760 0 : sogi_mass(:,:)=0._r8
761 0 : sogt_mass(:,:)=0._r8
762 0 : sogb_mass(:,:)=0._r8
763 0 : sogx_mass(:,:)=0._r8
764 :
765 0 : call pbuf_get_field(pbuf, fracsoa_ndx, fracsoa )
766 0 : call pbuf_get_field(pbuf, fracsog_ndx, fracsog )
767 :
768 0 : do i=1,ncol
769 0 : do k=1,pver
770 : !
771 : ! INIALIZATION AND LUMPING
772 : !
773 : ! calculate mass concentration of air (ug/m3)
774 0 : m_air = xhnm(i,k)*28.966_r8/avogadro*1.e12_r8
775 : !
776 : ! calculate initial mass of POA from OC1 and OC2 (in ug/m3)
777 0 : poa = (vmr(i,k,oc1_ndx)+vmr(i,k,oc2_ndx)) * OMscale * xhnm(i,k) * adv_mass(oc1_ndx)/avogadro * 1.e12_r8
778 : !
779 : ! specify pre-existing SOG/SOA for each class (in ug/m3)
780 :
781 :
782 : sog0(1,1)=vmr(i,k,sogm_ndx) * xhnm(i,k) * adv_mass(sogm_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 1,1)
783 : sog0(1,2)=vmr(i,k,sogm_ndx) * xhnm(i,k) * adv_mass(sogm_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 1,2)
784 : sog0(2,1)=vmr(i,k,sogm_ndx) * xhnm(i,k) * adv_mass(sogm_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 2,1)
785 : sog0(2,2)=vmr(i,k,sogm_ndx) * xhnm(i,k) * adv_mass(sogm_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 2,2)
786 : sog0(3,1)=vmr(i,k,sogm_ndx) * xhnm(i,k) * adv_mass(sogm_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 3,1)
787 : sog0(3,2)=0._r8
788 : sog0(4,1)=vmr(i,k,sogi_ndx) * xhnm(i,k) * adv_mass(sogi_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 4,1)
789 : sog0(4,2)=vmr(i,k,sogi_ndx) * xhnm(i,k) * adv_mass(sogi_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 4,2)
790 : sog0(5,1)=vmr(i,k,sogt_ndx) * xhnm(i,k) * adv_mass(sogt_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 5,1)
791 : sog0(5,2)=0._r8
792 : sog0(6,1)=vmr(i,k,sogt_ndx) * xhnm(i,k) * adv_mass(sogt_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 6,1)
793 : sog0(6,2)=vmr(i,k,sogt_ndx) * xhnm(i,k) * adv_mass(sogt_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 6,2)
794 : sog0(7,1)=vmr(i,k,sogb_ndx) * xhnm(i,k) * adv_mass(sogb_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 7,1)
795 : sog0(7,2)=0._r8
796 : sog0(8,1)=vmr(i,k,sogb_ndx) * xhnm(i,k) * adv_mass(sogb_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 8,1)
797 : sog0(8,2)=vmr(i,k,sogb_ndx) * xhnm(i,k) * adv_mass(sogb_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 8,2)
798 : sog0(9,1)=vmr(i,k,sogx_ndx) * xhnm(i,k) * adv_mass(sogx_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 9,1)
799 : sog0(9,2)=0._r8
800 : sog0(10,1)=vmr(i,k,sogx_ndx) * xhnm(i,k) * adv_mass(sogx_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 10,1)
801 : sog0(10,2)=vmr(i,k,sogx_ndx) * xhnm(i,k) * adv_mass(sogx_ndx)/avogadro * 1.e12_r8 * fracsog(i,k, 10,2)
802 : !
803 : soa0(1,1)=vmr(i,k,soam_ndx) * xhnm(i,k) * adv_mass(soam_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 1,1)
804 : soa0(1,2)=vmr(i,k,soam_ndx) * xhnm(i,k) * adv_mass(soam_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 1,2)
805 : soa0(2,1)=vmr(i,k,soam_ndx) * xhnm(i,k) * adv_mass(soam_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 2,1)
806 : soa0(2,2)=vmr(i,k,soam_ndx) * xhnm(i,k) * adv_mass(soam_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 2,2)
807 : soa0(3,1)=vmr(i,k,soam_ndx) * xhnm(i,k) * adv_mass(soam_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 3,1)
808 : soa0(3,2)=0._r8
809 : soa0(4,1)=vmr(i,k,soai_ndx) * xhnm(i,k) * adv_mass(soai_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 4,1)
810 : soa0(4,2)=vmr(i,k,soai_ndx) * xhnm(i,k) * adv_mass(soai_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 4,2)
811 : soa0(5,1)=vmr(i,k,soat_ndx) * xhnm(i,k) * adv_mass(soat_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 5,1)
812 : soa0(5,2)=0._r8
813 : soa0(6,1)=vmr(i,k,soat_ndx) * xhnm(i,k) * adv_mass(soat_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 6,1)
814 : soa0(6,2)=vmr(i,k,soat_ndx) * xhnm(i,k) * adv_mass(soat_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 6,2)
815 : soa0(7,1)=vmr(i,k,soab_ndx) * xhnm(i,k) * adv_mass(soab_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 7,1)
816 : soa0(7,2)=0._r8
817 : soa0(8,1)=vmr(i,k,soab_ndx) * xhnm(i,k) * adv_mass(soab_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 8,1)
818 : soa0(8,2)=vmr(i,k,soab_ndx) * xhnm(i,k) * adv_mass(soab_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 8,2)
819 : soa0(9,1)=vmr(i,k,soax_ndx) * xhnm(i,k) * adv_mass(soax_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 9,1)
820 : soa0(9,2)=0._r8
821 : soa0(10,1)=vmr(i,k,soax_ndx) * xhnm(i,k) * adv_mass(soax_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 10,1)
822 : soa0(10,2)=vmr(i,k,soax_ndx) * xhnm(i,k) * adv_mass(soax_ndx)/avogadro * 1.e12_r8 * fracsoa(i,k, 10,2)
823 : !
824 : !
825 : !--------------------------------------------------------------
826 : ! CHEMISTRY
827 : !
828 : sumorg=0._r8
829 : !
830 : do n=1,NRX
831 : !
832 : ! temperature dependence of paritioning via Clausiaus Clayperon [C&S, 2002]: clh (02/24/06)
833 : T_fac = tfld(i,k)/T1(n) * exp(delH(n)/Rgas*(1/tfld(i,k) - 1/T1(n)))
834 : k_om_T(n,:)=k_om(n,:)*T_fac
835 : !
836 : ! find molecular weight of SOA
837 : if ( (n >=1) .and. (n < 4) ) then
838 : mw_soa = adv_mass(soam_ndx)
839 : else if (n == 4) then
840 : mw_soa = adv_mass(soai_ndx)
841 : else if ( (n > 4) .and. (n < 7) ) then
842 : mw_soa = adv_mass(soat_ndx)
843 : else if ( (n >= 7) .and. (n < 9) ) then
844 : mw_soa = adv_mass(soab_ndx)
845 : else if ( (n >= 9) .and. (n < 11) ) then
846 : mw_soa = adv_mass(soax_ndx)
847 : end if
848 : !
849 : ! define chemical production from gas-phase chemistry (and convert to ug/m3)
850 : delHC = reaction_rates(i,k,rxn_soa(n)) &
851 : * vmr(i,k,react_ndx(n,1)) * vmr(i,k,react_ndx(n,2))* dt &
852 : * xhnm(i,k) * mw_soa/avogadro * 1.e12_r8
853 : !
854 : ! specify the new total of gas-phase products (before re-partitioning)
855 : ! and total up all SOA in gas and aerosol phase
856 : do p=1,NPR
857 : orggas(n,p) = sog0(n,p)+alpha(n,p)*delHC
858 : sumorg=sumorg+orggas(n,p)+soa0(n,p)
859 : prod(n,p)=alpha(n,p)*delHC
860 : enddo
861 : !
862 : enddo
863 : !
864 : !--------------------------------------------------------------
865 : ! check to see if no organics no partitioning!
866 : ! (set here to previous timestep concentration+oxprod to preserve mass)
867 : if (sumorg < 1.e-10_r8) then
868 : do n=1,NRX
869 : do p=1,NPR
870 : sog(n,p)=sog0(n,p)
871 : soa(n,p)=soa0(n,p)+prod(n,p)
872 : enddo
873 : enddo
874 : else
875 : ! PARTITION (Need to iteratively solve for MNEW, set tolerances to 0.1 ng/m3 or 1% of MNEW)
876 : !
877 : numer=0._r8
878 : maxM=0._r8
879 : ! If POA is essentially zero, equations simplify
880 : if (POA < 1.e-10_r8) then
881 : do n=1,NRX
882 : do p=1,NPR
883 : numer = numer + k_om_T(n,p)* (orggas(n,p)+soa0(n,p))
884 : enddo
885 : enddo
886 : ! if numerator is less than 1 then MNEW must be zero
887 : if (numer <= 1._r8) then
888 : mnew=0._r8
889 : iter=0
890 : else
891 : minM = 1.e-40_r8
892 : maxM = poa + sumorg
893 : tol = 1.e-4_r8
894 : mnew = zeroin(minM,maxM,tol,poa,soa0,orggas,k_om_T,sumorg,iter)
895 : end if
896 : !
897 : ! if additional organic mass is less than 1% of POA, or difference between the
898 : ! two is less than 0.1 ng/m3 (the tolerance) then POA is partitioning mass
899 : else if ( (sumorg < 0.01_r8*poa) .or. (abs(sumorg-poa) < 1.e-4_r8) ) then
900 : mnew = poa
901 : iter = 0
902 : ! otherwise solve for MNEW iteratively on interval [poa, poa+sumorg]
903 : else
904 : maxM = poa+sumorg
905 : minM = poa
906 : tol = 1.e-4_r8
907 : mnew = zeroin(minM,maxM,tol,poa,soa0,orggas,k_om_T,sumorg,iter)
908 : end if
909 : !
910 : !
911 : ! Now equilibrium partitioning with new MNEW
912 : ! If no MNEW then all the SOA evaporates to gas-phase
913 : if (mnew > 0._r8) then
914 : do n=1,NRX
915 : do p=1,NPR
916 : soa(n,p)=k_om_T(n,p)*mnew*(orggas(n,p)+soa0(n,p))/(1._r8+k_om_T(n,p)*mnew)
917 : if (k_om_T(n,p) /= 0._r8) then
918 : sog(n,p)=soa(n,p)/(k_om_T(n,p)*mnew)
919 : else
920 : sog(n,p)=0._r8
921 : end if
922 : enddo
923 : enddo
924 :
925 : else
926 : do n=1,NRX
927 : do p=1,NPR
928 : sog(n,p)=orggas(n,p)+soa0(n,p)
929 : soa(n,p)=1.e-20_r8
930 : enddo
931 : enddo
932 : end if
933 : !
934 : end if
935 : !
936 : !--------------------------------------------------------------
937 : ! LUMP INTO ARRAYS
938 : do n=1,NRX
939 : do p=1,NPR
940 : if ( (n >=1) .and. (n < 4) ) then
941 : soam_mass(i,k) = soam_mass(i,k) + soa(n,p)
942 : sogm_mass(i,k) = sogm_mass(i,k) + sog(n,p)
943 : else if (n == 4) then
944 : soai_mass(i,k) = soai_mass(i,k) + soa(n,p)
945 : sogi_mass(i,k) = sogi_mass(i,k) + sog(n,p)
946 : else if ( (n > 4) .and. (n < 7) ) then
947 : soat_mass(i,k) = soat_mass(i,k) + soa(n,p)
948 : sogt_mass(i,k) = sogt_mass(i,k) + sog(n,p)
949 : else if ( (n >= 7) .and. (n < 9) ) then
950 : soab_mass(i,k) = soab_mass(i,k) + soa(n,p)
951 : sogb_mass(i,k) = sogb_mass(i,k) + sog(n,p)
952 : else if ( (n >= 9) .and. (n < 11) ) then
953 : soax_mass(i,k) = soax_mass(i,k) + soa(n,p)
954 : sogx_mass(i,k) = sogx_mass(i,k) + sog(n,p)
955 : end if
956 : enddo
957 : enddo
958 : !
959 : ! Save mass fraction of each SOA rxn to SOA class
960 : ! (but if sumorg essentially zero revert to init fracs OR
961 : ! if mnew = 0 (all evap) then fracsoa revert to old, calculate fracsog only)
962 : if (sumorg < 1.e-10_r8) then
963 : do n=1,NRX
964 : do p=1,NPR
965 : fracsoa(i,k,n,p)=fracsoa_init(n,p)
966 : fracsog(i,k,n,p)=fracsog_init(n,p)
967 : enddo
968 : enddo
969 : else
970 : if (mnew < 1.e-20_r8) then
971 : do n=1,NRX
972 : do p=1,NPR
973 : fracsoa(i,k,n,p)=fracsoa_init(n,p)
974 : if ( (n >=1) .and. (n < 4) ) then
975 : if (sogm_mass(i,k) == 0._r8) then
976 : fracsog(i,k,n,p)=fracsog_init(n,p)
977 : else
978 : fracsog(i,k,n,p)=sog(n,p)/sogm_mass(i,k)
979 : end if
980 : else if (n == 4) then
981 : if (sogi_mass(i,k) == 0._r8) then
982 : fracsog(i,k,n,p)=fracsog_init(n,p)
983 : else
984 : fracsog(i,k,n,p)=sog(n,p)/sogi_mass(i,k)
985 : end if
986 : else if ( (n > 4) .and. (n < 7) ) then
987 : if (sogt_mass(i,k) == 0._r8) then
988 : fracsog(i,k,n,p)=fracsog_init(n,p)
989 : else
990 : fracsog(i,k,n,p)=sog(n,p)/sogt_mass(i,k)
991 : end if
992 : else if ( (n >= 7) .and. (n < 9) ) then
993 : if (sogb_mass(i,k) == 0._r8) then
994 : fracsog(i,k,n,p)=fracsog_init(n,p)
995 : else
996 : fracsog(i,k,n,p)=sog(n,p)/sogb_mass(i,k)
997 : end if
998 : else if ( (n >= 9) .and. (n < 11) ) then
999 : if (sogx_mass(i,k) == 0._r8) then
1000 : fracsog(i,k,n,p)=fracsog_init(n,p)
1001 : else
1002 : fracsog(i,k,n,p)=sog(n,p)/sogx_mass(i,k)
1003 : end if
1004 : end if
1005 : if ( (p==2) .and. (n==3 .or. n==5 .or. n==7 .or. n==9) ) then
1006 : fracsog(i,k,n,p)=0._r8
1007 : end if
1008 : enddo
1009 : enddo
1010 : else
1011 : do n=1,NRX
1012 : do p=1,NPR
1013 : if ( (n >=1) .and. (n < 4) ) then
1014 : if (soam_mass(i,k) == 0._r8) then
1015 : fracsoa(i,k,n,p)=fracsoa_init(n,p)
1016 : else
1017 : fracsoa(i,k,n,p)=soa(n,p)/soam_mass(i,k)
1018 : end if
1019 : if (sogm_mass(i,k) == 0._r8) then
1020 : fracsog(i,k,n,p)=fracsog_init(n,p)
1021 : else
1022 : fracsog(i,k,n,p)=sog(n,p)/sogm_mass(i,k)
1023 : end if
1024 : else if (n == 4) then
1025 : if (soai_mass(i,k) == 0._r8) then
1026 : fracsoa(i,k,n,p)=fracsoa_init(n,p)
1027 : else
1028 : fracsoa(i,k,n,p)=soa(n,p)/soai_mass(i,k)
1029 : end if
1030 : if (sogi_mass(i,k) == 0._r8) then
1031 : fracsog(i,k,n,p)=fracsog_init(n,p)
1032 : else
1033 : fracsog(i,k,n,p)=sog(n,p)/sogi_mass(i,k)
1034 : end if
1035 : else if ( (n > 4) .and. (n < 7) ) then
1036 : if (soat_mass(i,k) == 0._r8) then
1037 : fracsoa(i,k,n,p)=fracsoa_init(n,p)
1038 : else
1039 : fracsoa(i,k,n,p)=soa(n,p)/soat_mass(i,k)
1040 : end if
1041 : if (sogt_mass(i,k) == 0._r8) then
1042 : fracsog(i,k,n,p)=fracsog_init(n,p)
1043 : else
1044 : fracsog(i,k,n,p)=sog(n,p)/sogt_mass(i,k)
1045 : end if
1046 : else if ( (n >= 7) .and. (n < 9) ) then
1047 : if (soab_mass(i,k) == 0._r8) then
1048 : fracsoa(i,k,n,p)=fracsoa_init(n,p)
1049 : else
1050 : fracsoa(i,k,n,p)=soa(n,p)/soab_mass(i,k)
1051 : end if
1052 : if (sogb_mass(i,k) == 0._r8) then
1053 : fracsog(i,k,n,p)=fracsog_init(n,p)
1054 : else
1055 : fracsog(i,k,n,p)=sog(n,p)/sogb_mass(i,k)
1056 : end if
1057 : else if ( (n >= 9) .and. (n < 11) ) then
1058 : if (soax_mass(i,k) == 0._r8) then
1059 : fracsoa(i,k,n,p)=fracsoa_init(n,p)
1060 : else
1061 : fracsoa(i,k,n,p)=soa(n,p)/soax_mass(i,k)
1062 : end if
1063 : if (sogx_mass(i,k) == 0._r8) then
1064 : fracsog(i,k,n,p)=fracsog_init(n,p)
1065 : else
1066 : fracsog(i,k,n,p)=sog(n,p)/sogx_mass(i,k)
1067 : end if
1068 : end if
1069 : if ( (p==2) .and. (n==3 .or. n==5 .or. n==7 .or. n==9) ) then
1070 : fracsoa(i,k,n,p)=0._r8
1071 : fracsog(i,k,n,p)=0._r8
1072 : end if
1073 : enddo
1074 : enddo
1075 : end if
1076 : end if
1077 : !
1078 : !--------------------------------------------------------------
1079 : !
1080 : ! calculate NET production in kg/kg/s (subtract initial mass)
1081 : !
1082 : soam_prod(i,k) = ( soam_mass(i,k)*1.e-12_r8*avogadro/(adv_mass(soam_ndx)*xhnm(i,k)) - &
1083 : vmr(i,k,soam_ndx) )/dt
1084 : soai_prod(i,k) = ( soai_mass(i,k)*1.e-12_r8*avogadro/(adv_mass(soai_ndx)*xhnm(i,k)) - &
1085 : vmr(i,k,soai_ndx) )/dt
1086 : soat_prod(i,k) = ( soat_mass(i,k)*1.e-12_r8*avogadro/(adv_mass(soat_ndx)*xhnm(i,k)) - &
1087 : vmr(i,k,soat_ndx) )/dt
1088 : soab_prod(i,k) = ( soab_mass(i,k)*1.e-12_r8*avogadro/(adv_mass(soab_ndx)*xhnm(i,k)) - &
1089 : vmr(i,k,soab_ndx) )/dt
1090 : soax_prod(i,k) = ( soax_mass(i,k)*1.e-12_r8*avogadro/(adv_mass(soax_ndx)*xhnm(i,k)) - &
1091 : vmr(i,k,soax_ndx) )/dt
1092 : !
1093 : ! convert from ug/m3 to mixing ratio and update vmr
1094 : !
1095 : vmr(i,k,soam_ndx) = soam_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(soam_ndx)*xhnm(i,k))
1096 : vmr(i,k,soai_ndx) = soai_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(soai_ndx)*xhnm(i,k))
1097 : vmr(i,k,soat_ndx) = soat_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(soat_ndx)*xhnm(i,k))
1098 : vmr(i,k,soab_ndx) = soab_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(soab_ndx)*xhnm(i,k))
1099 : vmr(i,k,soax_ndx) = soax_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(soax_ndx)*xhnm(i,k))
1100 : vmr(i,k,sogm_ndx) = sogm_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(sogm_ndx)*xhnm(i,k))
1101 : vmr(i,k,sogi_ndx) = sogi_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(sogi_ndx)*xhnm(i,k))
1102 : vmr(i,k,sogt_ndx) = sogt_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(sogt_ndx)*xhnm(i,k))
1103 : vmr(i,k,sogb_ndx) = sogb_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(sogb_ndx)*xhnm(i,k))
1104 : vmr(i,k,sogx_ndx) = sogx_mass(i,k) * 1.e-12_r8 * avogadro/(adv_mass(sogx_ndx)*xhnm(i,k))
1105 : enddo
1106 : enddo
1107 : !
1108 0 : call outfld('SOAM_PROD',soam_prod(:ncol,:),ncol, lchnk)
1109 0 : call outfld('SOAI_PROD',soai_prod(:ncol,:),ncol, lchnk)
1110 0 : call outfld('SOAT_PROD',soat_prod(:ncol,:),ncol, lchnk)
1111 0 : call outfld('SOAB_PROD',soab_prod(:ncol,:),ncol, lchnk)
1112 0 : call outfld('SOAX_PROD',soax_prod(:ncol,:),ncol, lchnk)
1113 :
1114 0 : call outfld('SOAM_dens',soam_mass(:ncol,:), ncol, lchnk)
1115 0 : call outfld('SOAI_dens',soai_mass(:ncol,:), ncol, lchnk)
1116 0 : call outfld('SOAT_dens',soat_mass(:ncol,:), ncol, lchnk)
1117 0 : call outfld('SOAB_dens',soab_mass(:ncol,:), ncol, lchnk)
1118 0 : call outfld('SOAX_dens',soax_mass(:ncol,:), ncol, lchnk)
1119 :
1120 : !
1121 : !
1122 0 : return
1123 0 : end subroutine setsoa_equil
1124 :
1125 : !===============================================================================
1126 : !===============================================================================
1127 : real(r8) function zeroin(x1,x2,tol,poa,aer,gas,k,totorg,iter)
1128 : ! function to iteratively solve function using bilinear method
1129 : !
1130 : implicit none
1131 : !
1132 : integer :: iter
1133 : real(r8),intent(in) :: x1, x2 ! min/max of interval
1134 : real(r8),intent(in) :: tol ! tolerance (interval of uncertainty)
1135 : real(r8),intent(in) :: poa,totorg
1136 : real(r8),intent(in) :: aer(NRX,NPR), gas(NRX,NPR) ! aerosol and gas phase concentrations
1137 : real(r8),intent(in) :: k(NRX,NPR) ! partitioning coeff
1138 : ! local vars
1139 : real(r8) :: xa,xb,xm,fa,fb,fm
1140 : !
1141 : xa=x1
1142 : xb=x2
1143 : xm=0._r8
1144 : fa=soa_function(xa,poa,aer,gas,k)
1145 : fb=soa_function(xb,poa,aer,gas,k)
1146 : !
1147 : ! check that functions have opposite signs
1148 : if (fa >= 0._r8) then
1149 : if (fb >=0._r8) then
1150 : write(iulog,*) 'ABORT IN ZEROIN: SAME SIGN ON FUNCTION',poa,totorg,x1,x2,fa,fb,aer,gas,k
1151 : write(iulog,*) 'ABORT IN ZEROIN: ERROR1: fa, fb ',fa, fb
1152 : write(iulog,*) 'ABORT IN ZEROIN: ERROR1: maxval(aer),minval(aer),maxval(gas),minval(gas) ',&
1153 : maxval(aer),minval(aer),maxval(gas),minval(gas)
1154 : call endrun('ABORT IN ZEROIN: ERROR1')
1155 : end if
1156 : else
1157 : if (fb <=0._r8) then
1158 : write(iulog,*) 'ABORT IN ZEROIN: SAME SIGN ON FUNCTION',poa,totorg,x1,x2,fa,fb,aer,gas,k
1159 : write(iulog,*) 'ABORT IN ZEROIN: ERROR2: fa, fb ',fa, fb
1160 : write(iulog,*) 'ABORT IN ZEROIN: ERROR2: maxval(aer),minval(aer),maxval(gas),minval(gas) ',&
1161 : maxval(aer),minval(aer),maxval(gas),minval(gas)
1162 : call endrun('ABORT IN ZEROIN: ERROR2')
1163 : end if
1164 : end if
1165 : !
1166 : iter=0
1167 : do while ((abs(xa-xb) > 2._r8*tol) .and. (abs(xa-xb) > 0.01_r8*xm) )
1168 : xm=(xa+xb)/2
1169 : fm=soa_function(xm,poa,aer,gas,k)
1170 : if (fa >=0._r8) then
1171 : if (fm >=0._r8) then
1172 : xa=xm
1173 : fa=fm
1174 : else
1175 : xb=xm
1176 : fb=fm
1177 : end if
1178 : else
1179 : if (fm < 0._r8) then
1180 : xa=xm
1181 : fa=fm
1182 : else
1183 : xb=xm
1184 : fb=fm
1185 : end if
1186 : end if
1187 : iter=iter+1
1188 : enddo
1189 : !
1190 : zeroin = (xa+xb)/2
1191 : !
1192 : return
1193 0 : end function zeroin
1194 : !===============================================================================
1195 : !===============================================================================
1196 : real(r8) function soa_function(m0,poa,aer,gas,k)
1197 : ! function which calculates SOAeqn (trying to minimize to zero)
1198 : !
1199 : implicit none
1200 : !
1201 : real(r8),intent(in) :: m0,poa
1202 : real(r8),intent(in) :: aer(NRX,NPR), gas(NRX,NPR) ! aerosol and gas phase concentrations
1203 : real(r8),intent(in) :: k(NRX,NPR) ! partitioning coeff
1204 : ! local vars
1205 : integer :: n,p
1206 : real(r8) :: value
1207 : !
1208 : value=0._r8
1209 : !
1210 : do n=1,NRX
1211 : do p=1,NPR
1212 : value = value + k(n,p)*(gas(n,p)+aer(n,p))/(1._r8 + k(n,p)*m0)
1213 : enddo
1214 : enddo
1215 : soa_function = value + (poa/m0) - 1._r8
1216 : !
1217 : ! write(iulog,*) 'clh soa_function out: ',m0,soa_function
1218 : return
1219 : end function soa_function
1220 : !===============================================================================
1221 :
1222 : end module mo_setsoa
|
