| blob | d0d02ffe7611b0f5901313d850c7bc6875ce5f34 |
1 !Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
2 #include "MODAL_AERO_CPP_DEFINES.h"
3 #define WRF_PORT
4 #define MODAL_AERO
6 module MO_SETSOX
8 use shr_kind_mod, only : r8 => shr_kind_r8
9 #ifndef WRF_PORT
10 use cam_logfile, only : iulog
11 #else
12 use module_cam_support, only: iulog
13 #endif
14 #if (defined MODAL_AERO)
15 use modal_aero_data
16 #endif
18 private
19 public :: sox_inti, setsox
20 public :: has_sox
22 save
24 #if (defined MODAL_AERO)
25 #if ( defined MODAL_AERO_7MODE )
26 integer, target :: spc_ids(18)
27 integer, pointer :: id_so2, id_so4_1, id_so4_2, id_so4_3, id_so4_4, id_so4_5, id_so4_6, &
28 id_nh3, id_nh4_1, id_nh4_2, id_nh4_3, id_nh4_4, id_nh4_5, id_nh4_6, &
29 id_h2o2, id_ox, id_ho2, id_h2so4
30 #elif ( defined MODAL_AERO_3MODE )
31 integer, target :: spc_ids(8)
32 integer, pointer :: id_so2, id_so4_1, id_so4_2, id_so4_3, &
33 id_h2o2, id_ox, id_ho2, id_h2so4
34 integer :: id_nh3
35 #endif
36 logical :: inv_o3
37 integer :: id_hno3, id_msa
38 #else
39 integer, target :: spc_ids(8)
40 integer, pointer :: id_so2, id_so4, id_nh3, id_hno3, id_h2o2, id_ox, id_nh4no3, id_ho2
41 #endif
42 logical :: has_sox = .true.
43 logical :: inv_so2, inv_nh3, inv_hno3, inv_h2o2, inv_ox, inv_nh4no3, inv_ho2
44 contains
46 subroutine sox_inti
47 !-----------------------------------------------------------------------
48 ! ... initialize the hetero sox routine
49 !-----------------------------------------------------------------------
51 use mo_chem_utls, only : get_spc_ndx, get_inv_ndx
52 #ifndef WRF_PORT
53 use spmd_utils, only : masterproc
54 #else
55 use module_cam_support, only: masterproc
56 #endif
58 implicit none
60 #if (defined MODAL_AERO)
61 #if ( defined MODAL_AERO_7MODE )
62 id_so2 => spc_ids(1)
63 id_so4_1 => spc_ids(2)
64 id_so4_2 => spc_ids(3)
65 id_so4_3 => spc_ids(4)
66 id_so4_4 => spc_ids(5)
67 id_so4_5 => spc_ids(6)
68 id_so4_6 => spc_ids(7)
70 id_nh3 => spc_ids(8)
71 id_nh4_1 => spc_ids(9)
72 id_nh4_2 => spc_ids(10)
73 id_nh4_3 => spc_ids(11)
74 id_nh4_4 => spc_ids(12)
75 id_nh4_5 => spc_ids(13)
76 id_nh4_6 => spc_ids(14)
78 id_h2o2 => spc_ids(15)
79 id_ox => spc_ids(16)
80 id_ho2 => spc_ids(17)
81 id_h2so4 => spc_ids(18)
82 #elif ( defined MODAL_AERO_3MODE )
83 id_so2 => spc_ids(1)
84 id_so4_1 => spc_ids(2)
85 id_so4_2 => spc_ids(3)
86 id_so4_3 => spc_ids(4)
88 id_h2o2 => spc_ids(5)
89 id_ox => spc_ids(6)
90 id_ho2 => spc_ids(7)
91 id_h2so4 => spc_ids(8)
92 #endif
93 #else
94 id_so2 => spc_ids(1)
95 id_so4 => spc_ids(2)
96 id_nh3 => spc_ids(3)
97 id_hno3 => spc_ids(4)
98 id_h2o2 => spc_ids(5)
99 id_ox => spc_ids(6)
100 id_ho2 => spc_ids(7)
101 #endif
102 !-----------------------------------------------------------------
103 ! ... get species indicies
104 !-----------------------------------------------------------------
105 #if (defined MODAL_AERO)
106 #if ( defined MODAL_AERO_7MODE )
107 id_so4_1 = lptr_so4_cw_amode(1)
108 id_so4_2 = lptr_so4_cw_amode(2)
109 id_so4_3 = lptr_so4_cw_amode(4) ! no so4_c3
110 id_so4_4 = lptr_so4_cw_amode(5)
111 id_so4_5 = lptr_so4_cw_amode(6)
112 id_so4_6 = lptr_so4_cw_amode(7)
113 #elif ( defined MODAL_AERO_3MODE )
114 id_so4_1 = lptr_so4_cw_amode(1)
115 id_so4_2 = lptr_so4_cw_amode(2)
116 id_so4_3 = lptr_so4_cw_amode(3)
117 #endif
119 id_h2so4 = get_spc_ndx( 'H2SO4' )
120 id_msa = get_spc_ndx( 'MSA' )
121 #else
122 id_so4 = get_spc_ndx( 'SO4' )
123 #endif
125 inv_so2 = .false.
126 id_so2 = get_inv_ndx( 'SO2' )
127 inv_so2 = id_so2 > 0
128 if ( .not. inv_so2 ) then
129 id_so2 = get_spc_ndx( 'SO2' )
130 endif
132 inv_NH3 = .false.
133 id_NH3 = get_inv_ndx( 'NH3' )
134 inv_NH3 = id_NH3 > 0
135 if ( .not. inv_NH3 ) then
136 id_NH3 = get_spc_ndx( 'NH3' )
137 endif
139 #if (defined MODAL_AERO)
140 #if ( defined MODAL_AERO_7MODE )
141 id_nh4_1 = lptr_nh4_cw_amode(1)
142 id_nh4_2 = lptr_nh4_cw_amode(2)
143 id_nh4_3 = lptr_nh4_cw_amode(4) ! no nh4_c3
144 id_nh4_4 = lptr_nh4_cw_amode(5)
145 id_nh4_5 = lptr_nh4_cw_amode(6)
146 id_nh4_6 = lptr_nh4_cw_amode(7)
147 #endif
148 #endif
150 inv_HNO3 = .false.
151 id_HNO3 = get_inv_ndx( 'HNO3' )
152 inv_HNO3 = id_hno3 > 0
153 if ( .not. inv_HNO3 ) then
154 id_HNO3 = get_spc_ndx( 'HNO3' )
155 endif
157 inv_H2O2 = .false.
158 id_H2O2 = get_inv_ndx( 'H2O2' )
159 inv_H2O2 = id_H2O2 > 0
160 if ( .not. inv_H2O2 ) then
161 id_H2O2 = get_spc_ndx( 'H2O2' )
162 endif
164 inv_HO2 = .false.
165 id_HO2 = get_inv_ndx( 'HO2' )
166 inv_HO2 = id_HO2 > 0
167 if ( .not. inv_HO2 ) then
168 id_HO2 = get_spc_ndx( 'HO2' )
169 endif
171 #if (defined MODAL_AERO)
172 inv_o3 = get_inv_ndx( 'O3' ) > 0
173 if (inv_o3) then
174 id_ox = get_inv_ndx( 'O3' )
175 else
176 id_ox = get_spc_ndx( 'O3' )
177 endif
178 inv_ho2 = get_inv_ndx( 'HO2' ) > 0
179 if (inv_ho2) then
180 id_ho2 = get_inv_ndx( 'HO2' )
181 else
182 id_ho2 = get_spc_ndx( 'HO2' )
183 endif
184 #else
185 inv_OX = .false.
186 id_OX = get_inv_ndx( 'OX' )
187 if( id_ox < 1 ) then
188 id_ox = get_inv_ndx( 'O3' )
189 end if
190 inv_OX = id_OX > 0
191 if ( .not. inv_OX ) then
192 id_OX = get_spc_ndx( 'OX' )
193 if ( id_OX < 0 ) then
194 id_OX = get_spc_ndx( 'O3' )
195 endif
196 endif
197 #endif
199 #if (defined MODAL_AERO)
200 #if ( defined MODAL_AERO_7MODE )
201 has_sox = all( spc_ids(1:18) > 0 )
202 #elif ( defined MODAL_AERO_3MODE )
203 has_sox = all( spc_ids(1:8) > 0 )
204 #endif
205 #else
206 has_sox = all( spc_ids(1:7) > 0 )
207 #endif
209 if (masterproc) then
210 write(iulog,*) 'sox_inti: has_sox = ',has_sox
211 #ifdef WRF_PORT
212 call wrf_message(iulog)
213 #endif
214 endif
216 if( has_sox ) then
217 if (masterproc) then
218 write(iulog,*) '-----------------------------------------'
219 #ifdef WRF_PORT
220 call wrf_message(iulog)
221 #endif
222 write(iulog,*) 'mozart will do sox aerosols'
223 #ifdef WRF_PORT
224 call wrf_message(iulog)
225 #endif
226 write(iulog,*) '-----------------------------------------'
227 #ifdef WRF_PORT
228 call wrf_message(iulog)
229 #endif
230 endif
231 end if
233 end subroutine sox_inti
235 subroutine SETSOX( ncol, &
236 press, &
237 dtime, &
238 tfld, &
239 qfld, &
240 lwc, &
241 #if (defined MODAL_AERO)
242 lchnk, &
243 pdel, &
244 mbar, &
245 clwlrat,&
246 cldfrc, &
247 cldnum, &
248 qcw, &
249 loffset,&
250 #endif
251 xhnm, &
252 qin, invariants )
254 !-----------------------------------------------------------------------
255 ! ... Compute heterogeneous reactions of SOX
256 !
257 ! (0) using initial PH to calculate PH
258 ! (a) HENRYs law constants
259 ! (b) PARTIONING
260 ! (c) PH values
261 !
262 ! (1) using new PH to repeat
263 ! (a) HENRYs law constants
264 ! (b) PARTIONING
265 ! (c) REACTION rates
266 ! (d) PREDICTION
267 !-----------------------------------------------------------------------
268 !
269 #ifndef WRF_PORT
270 use ppgrid, only : pcols, pver
271 use chem_mods, only : gas_pcnst, nfs
272 #else
273 use module_cam_support, only: pcols, pver, gas_pcnst => gas_pcnst_modal_aero, &
274 nfs
275 #endif
276 #if (defined MODAL_AERO)
277 #ifndef WRF_PORT
278 use chem_mods, only : adv_mass
279 #else
280 use module_data_cam_mam_asect, only: adv_mass => mw_q_mo_array
281 #endif
282 use physconst, only : mwdry, gravit
283 #ifndef WRF_PORT
284 use mo_constants, only : pi
285 use cam_history, only : outfld
286 #else
287 use physconst, only : pi
288 use module_cam_support, only: outfld
289 #endif
290 #endif
291 !
292 implicit none
293 !
294 !-----------------------------------------------------------------------
295 ! ... Dummy arguments
296 !-----------------------------------------------------------------------
297 integer, intent(in) :: ncol
298 real(r8), intent(in) :: dtime ! time step (sec)
299 real(r8), dimension(ncol ,pver,gas_pcnst), intent(inout) :: qin ! xported species ( vmr )
300 real(r8), dimension(ncol ,pver), intent(in) :: xhnm ! total atms density ( /cm**3)
301 real(r8), dimension(pcols,pver), intent(in) :: tfld, & ! temperature
302 qfld, & ! specific humidity( kg/kg )
303 press ! midpoint pressure ( Pa )
304 real(r8), dimension(ncol ,pver), intent(in) :: lwc !!, & ! cloud liquid water content (kg/kg)
305 real(r8), intent(in) :: invariants(ncol,pver,nfs)
306 #if (defined MODAL_AERO)
307 integer, intent(in) :: lchnk ! chunk id
308 real(r8), dimension(pcols,pver), intent(in) :: pdel ! pressure thickness of levels (Pa)
309 real(r8), dimension(ncol ,pver), intent(in) :: mbar ! mean wet atmospheric mass ( amu )
310 real(r8), dimension(ncol ,pver), intent(in) :: cldfrc, & ! cloud fraction
311 clwlrat,& ! first-order loss rate of l.s. condensate to precip (1/s)
312 cldnum ! droplet number concentration (#/kg)
313 real(r8), intent(inout) :: qcw(ncol,pver,gas_pcnst) ! cloud-borne aerosol (vmr)
314 integer, intent(in) :: loffset ! offset applied to modal aero "ptrs"
315 #endif
316 !-----------------------------------------------------------------------
317 ! ... Local variables
318 !
319 ! xhno3 ... in mixing ratio
320 !-----------------------------------------------------------------------
321 integer, parameter :: itermax = 20
322 real(r8), parameter :: ph0 = 5.0_r8 ! INITIAL PH VALUES
323 real(r8), parameter :: const0 = 1.e3_r8/6.023e23_r8
324 real(r8), parameter :: xa0 = 11._r8, &
325 xb0 = -.1_r8, &
326 xa1 = 1.053_r8, &
327 xb1 = -4.368_r8,&
328 xa2 = 1.016_r8, &
329 xb2 = -2.54_r8, &
330 xa3 = .816e-32_r8, &
331 xb3 = .259_r8
333 real(r8), parameter :: kh0 = 9.e3_r8, & ! HO2(g) -> Ho2(a)
334 kh1 = 2.05e-5_r8, & ! HO2(a) -> H+ + O2-
335 kh2 = 8.6e5_r8, & ! HO2(a) + ho2(a) -> h2o2(a) + o2
336 kh3 = 1.e8_r8, & ! HO2(a) + o2- -> h2o2(a) + o2
337 Ra = 8314._r8/101325._r8, & ! universal constant (atm)/(M-K)
338 xkw = 1.e-14_r8 ! water acidity
339 real(r8), parameter :: small_value = 1.e-20_r8
340 !
341 #if (defined MODAL_AERO)
342 integer :: l, n, m
343 integer :: ntot_msa_c
345 integer :: id_so4_1a, id_so4_2a, id_so4_3a, id_so4_4a, id_so4_5a, id_so4_6a, &
346 id_nh4_1a, id_nh4_2a, id_nh4_3a, id_nh4_4a, id_nh4_5a, id_nh4_6a
348 real(r8) :: delso4_hprxn, delso4_o3rxn, &
349 dso4dt_aqrxn, dso4dt_hprxn, &
350 dso4dt_gasuptk, dmsadt_gasuptk, &
351 dmsadt_gasuptk_tomsa, dmsadt_gasuptk_toso4, &
352 dqdt_aq, dqdt_wr, dqdt, &
353 rad_cd, radxnum_cd, num_cd, &
354 gasdiffus, gasspeed, knudsen, uptkrate, &
355 fuchs_sutugin, volx34pi_cd, &
356 fwetrem, sumf, &
357 frso2_g, frso2_c, frh2o2_g, frh2o2_c
358 real(r8) :: delnh3, delnh4
360 real(r8) :: faqgain_msa(ntot_amode), faqgain_so4(ntot_amode), &
361 qnum_c(ntot_amode)
362 real(r8) :: dqdt_aqso4(ncol,pver,gas_pcnst), &
363 dqdt_aqh2so4(ncol,pver,gas_pcnst), &
364 dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver), &
365 xphlwc(ncol,pver), &
366 sflx(1:ncol)
367 #endif
368 integer :: k, i, iter, file
369 real(r8) :: wrk, delta
370 real(r8) :: xph0, aden, xk, xe, x2
371 real(r8) :: tz, xl, px, qz, pz, es, qs, patm
372 real(r8) :: Eso2, Eso4, Ehno3, Eco2, Eh2o, Enh3
373 real(r8) :: hno3g, nh3g, so2g, h2o2g, co2g, o3g
374 real(r8) :: hno3a, nh3a, so2a, h2o2a, co2a, o3a
375 real(r8) :: rah2o2, rao3, pso4, ccc
376 real(r8) :: cnh3, chno3, com, com1, com2, xra
377 !
378 !-----------------------------------------------------------------------
379 ! for Ho2(g) -> H2o2(a) formation
380 ! schwartz JGR, 1984, 11589
381 !-----------------------------------------------------------------------
382 real(r8) :: kh4 ! kh2+kh3
383 real(r8) :: xam ! air density /cm3
384 real(r8) :: ho2s ! ho2s = ho2(a)+o2-
385 real(r8) :: r1h2o2 ! prod(h2o2) by ho2 in mole/L(w)/s
386 real(r8) :: r2h2o2 ! prod(h2o2) by ho2 in mix/s
388 real(r8), dimension(ncol,pver) :: &
389 xhno3, xh2o2, xso2, xso4,&
390 xnh3, xnh4, xo3, &
391 xlwc, cfact, &
392 xph, xho2, &
393 #if (defined MODAL_AERO)
394 xh2so4, xmsa, xso4_init, xso4c, xnh4c, &
395 #endif
396 hehno3, & ! henry law const for hno3
397 heh2o2, & ! henry law const for h2o2
398 heso2, & ! henry law const for so2
399 henh3, & ! henry law const for nh3
400 heo3 !!, & ! henry law const for o3
401 real(r8), dimension(ncol) :: work1
402 logical :: converged
404 #if (defined MODAL_AERO)
405 id_so4_1a = id_so4_1 - loffset
406 id_so4_2a = id_so4_2 - loffset
407 id_so4_3a = id_so4_3 - loffset
408 #if ( defined MODAL_AERO_7MODE )
409 id_so4_4a = id_so4_4 - loffset
410 id_so4_5a = id_so4_5 - loffset
411 id_so4_6a = id_so4_6 - loffset
413 id_nh4_1a = id_nh4_1 - loffset
414 id_nh4_2a = id_nh4_2 - loffset
415 id_nh4_3a = id_nh4_3 - loffset
416 id_nh4_4a = id_nh4_4 - loffset
417 id_nh4_5a = id_nh4_5 - loffset
418 id_nh4_6a = id_nh4_6 - loffset
419 #endif
421 ! make sure dqdt is zero initially, for budgets
422 dqdt_aqso4(:,:,:) = 0.0_r8
423 dqdt_aqh2so4(:,:,:) = 0.0_r8
424 dqdt_aqhprxn(:,:) = 0.0_r8
425 dqdt_aqo3rxn(:,:) = 0.0_r8
426 #endif
427 !-----------------------------------------------------------------
428 ! ... NOTE: The press array is in pascals and must be
429 ! mutiplied by 10 to yield dynes/cm**2.
430 !-----------------------------------------------------------------
431 !==================================================================
432 ! ... First set the PH
433 !==================================================================
434 ! ... Initial values
435 ! The values of so2, so4 are after (1) SLT, and CHEM
436 !-----------------------------------------------------------------
437 xph0 = 10._r8**(-ph0) ! initial PH value
439 do k = 1,pver
440 cfact(:,k) = xhnm(:,k) & ! /cm3(a)
441 * 1.e6_r8 & ! /m3(a)
442 * 1.38e-23_r8/287._r8 & ! Kg(a)/m3(a)
443 * 1.e-3_r8 ! Kg(a)/L(a)
444 end do
446 do k = 1,pver
447 xph(:,k) = xph0 ! initial PH value
448 xlwc(:,k) = lwc(:,k) *cfact(:,k) ! cloud water L(water)/L(air)
450 if ( inv_so2 ) then
451 xso2 (:,k) = invariants(:,k,id_so2) ! mixing ratio
452 else
453 xso2 (:,k) = qin(:,k,id_so2) ! mixing ratio
454 endif
456 #if (defined MODAL_AERO)
457 if (id_hno3 > 0) then
458 xhno3(:,k) = qin(:,k,id_hno3)
459 else
460 xhno3(:,k) = 0.0_r8
461 endif
462 #else
463 if ( inv_hno3 ) then
464 xhno3 (:,k) = invariants(:,k,id_hno3) ! mixing ratio
465 else
466 xhno3 (:,k) = qin(:,k,id_hno3) ! mixing ratio
467 endif
468 #endif
470 if ( inv_h2o2 ) then
471 xh2o2 (:,k) = invariants(:,k,id_h2o2) ! mixing ratio
472 else
473 xh2o2 (:,k) = qin(:,k,id_h2o2) ! mixing ratio
474 endif
476 #if (defined MODAL_AERO)
477 if (id_nh3 > 0) then
478 xnh3 (:,k) = qin(:,k,id_nh3)
479 else
480 xnh3 (:,k) = 0.0_r8
481 endif
482 #else
483 if ( inv_nh3 ) then
484 xnh3 (:,k) = invariants(:,k,id_nh3) ! mixing ratio
485 else
486 xnh3 (:,k) = qin(:,k,id_nh3) ! mixing ratio
487 endif
488 #endif
491 #if (defined MODAL_AERO)
492 if ( inv_o3 ) then
493 xo3 (:,k) = invariants(:,k,id_ox)/xhnm(:,k) ! mixing ratio
494 else
495 xo3 (:,k) = qin(:,k,id_ox) ! mixing ratio
496 endif
497 if ( inv_ho2 ) then
498 xho2 (:,k) = invariants(:,k,id_ho2)/xhnm(:,k)! mixing ratio
499 else
500 xho2 (:,k) = qin(:,k,id_ho2) ! mixing ratio
501 endif
502 #else
503 if ( inv_ox ) then
504 xo3 (:,k) = invariants(:,k,id_ox) ! mixing ratio
505 else
506 xo3 (:,k) = qin(:,k,id_ox) ! mixing ratio
507 endif
510 if ( inv_ho2 ) then
511 xho2 (:,k) = invariants(:,k,id_ho2) ! mixing ratio
512 else
513 xho2 (:,k) = qin(:,k,id_ho2) ! mixing ratio
514 endif
515 #endif
517 #if (defined MODAL_AERO)
518 #if ( defined MODAL_AERO_7MODE )
519 xso4c(:,k) = qcw(:,k,id_so4_1a) &
520 + qcw(:,k,id_so4_2a) &
521 + qcw(:,k,id_so4_3a) &
522 + qcw(:,k,id_so4_4a) &
523 + qcw(:,k,id_so4_5a) &
524 + qcw(:,k,id_so4_6a)
525 #elif ( defined MODAL_AERO_3MODE )
526 xso4c(:,k) = qcw(:,k,id_so4_1a) &
527 + qcw(:,k,id_so4_2a) &
528 + qcw(:,k,id_so4_3a)
529 #endif
530 xh2so4(:,k)= qin(:,k,id_h2so4)
531 if (id_msa > 0) xmsa (:,k) = qin(:,k,id_msa)
532 #else
533 xso4 (:,k) = qin(:,k,id_so4) ! mixing ratio
534 #endif
536 #if (defined MODAL_AERO)
537 #if ( defined MODAL_AERO_7MODE )
538 xnh4c(:,k) = qcw(:,k,id_nh4_1a) &
539 + qcw(:,k,id_nh4_2a) &
540 + qcw(:,k,id_nh4_3a) &
541 + qcw(:,k,id_nh4_4a) &
542 + qcw(:,k,id_nh4_5a) &
543 + qcw(:,k,id_nh4_6a)
544 #elif ( defined MODAL_AERO_3MODE )
545 xnh4c(:,k) = xso4c(:,k) * 1.0_r8 !assume NH4HSO4
546 #endif
547 #endif
548 end do
550 !-----------------------------------------------------------------
551 ! ... Temperature dependent Henry constants
552 !-----------------------------------------------------------------
553 do k = 1,pver !! pver loop for STEP 0
554 do i = 1,ncol
555 #if (defined MODAL_AERO)
556 if (cldfrc(i,k) >= 1.0e-5_r8) then
557 xso4(i,k) = xso4c(i,k) / cldfrc(i,k)
558 xnh4(i,k) = xnh4c(i,k) / cldfrc(i,k)
559 xl = xlwc(i,k) / cldfrc(i,k) ! liquid water in the cloudy fraction of cell
560 #else
561 xl = xlwc(i,k)
562 #endif
563 if( xl >= 1.e-8_r8 ) then
564 work1(i) = 1._r8 / tfld(i,k) - 1._r8 / 298._r8
565 !-----------------------------------------------------------------------
566 ! ... hno3
567 !-----------------------------------------------------------------------
568 do iter = 1,itermax
569 xk = 2.1e5_r8 *EXP( 8700._r8*work1(i) )
570 xe = 15.4_r8
571 hehno3(i,k) = xk*(1._r8 + xe/xph(i,k))
572 !-----------------------------------------------------------------------
573 ! ... h2o2
574 !-----------------------------------------------------------------------
575 xk = 7.4e4_r8 *EXP( 6621._r8*work1(i) )
576 xe = 2.2e-12_r8 *EXP(-3730._r8*work1(i) )
577 heh2o2(i,k) = xk*(1._r8 + xe/xph(i,k))
578 !-----------------------------------------------------------------------
579 ! ... so2
580 !-----------------------------------------------------------------------
581 xk = 1.23_r8 *EXP( 3120._r8*work1(i) )
582 xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) )
583 x2 = 6.0e-8_r8*EXP( 1120._r8*work1(i) )
584 wrk = xe/xph(i,k)
585 heso2(i,k) = xk*(1._r8 + wrk*(1._r8 + x2/xph(i,k)))
586 !-----------------------------------------------------------------------
587 ! ... nh3
588 !-----------------------------------------------------------------------
589 xk = 58._r8 *EXP( 4085._r8*work1(i) )
590 xe = 1.7e-5_r8*EXP(-4325._r8*work1(i) )
591 henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/xkw)
592 !-----------------------------------------------------------------
593 ! ... Partioning and effect of pH
594 !-----------------------------------------------------------------
595 pz = .01_r8*press(i,k) !! pressure in mb
596 tz = tfld(i,k)
597 patm = pz/1013._r8
598 xam = press(i,k)/(1.38e-23_r8*tz) !air density /M3
599 !-----------------------------------------------------------------
600 ! ... hno3
601 !-----------------------------------------------------------------
602 px = hehno3(i,k) * Ra * tz * xl
603 hno3g = xhno3(i,k)/(1._r8 + px)
604 xk = 2.1e5_r8 *EXP( 8700._r8*work1(i) )
605 xe = 15.4_r8
606 Ehno3 = xk*xe*hno3g *patm
607 !-----------------------------------------------------------------
608 ! ... so2
609 !-----------------------------------------------------------------
610 px = heso2(i,k) * Ra * tz * xl
611 so2g = xso2(i,k)/(1._r8+ px)
612 xk = 1.23_r8 *EXP( 3120._r8*work1(i) )
613 xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) )
614 Eso2 = xk*xe*so2g *patm
615 !-----------------------------------------------------------------
616 ! ... nh3
617 !-----------------------------------------------------------------
618 px = henh3(i,k) * Ra * tz * xl
619 #if (defined MODAL_AERO)
620 #if ( defined MODAL_AERO_7MODE )
621 nh3g = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px)
622 #elif ( defined MODAL_AERO_3MODE )
623 nh3g = xnh4(i,k)/px
624 #endif
625 #else
626 nh3g = xnh3(i,k)/(1._r8+ px)
627 #endif
628 xk = 58._r8 *EXP( 4085._r8*work1(i) )
629 xe = 1.7e-5_r8*EXP( -4325._r8*work1(i) )
630 Enh3 = xk*xe*nh3g/xkw *patm
631 !-----------------------------------------------------------------
632 ! ... h2o effects
633 !-----------------------------------------------------------------
634 Eh2o = xkw
635 !-----------------------------------------------------------------
636 ! ... co2 effects
637 !-----------------------------------------------------------------
638 co2g = 330.e-6_r8 !330 ppm = 330.e-6 atm
639 xk = 3.1e-2_r8*EXP( 2423._r8*work1(i) )
640 xe = 4.3e-7_r8*EXP(-913._r8 *work1(i) )
641 Eco2 = xk*xe*co2g *patm
642 !-----------------------------------------------------------------
643 ! ... PH cal
644 !-----------------------------------------------------------------
645 com2 = (Eh2o + Ehno3 + Eso2 + Eco2) &
646 / (1._r8 + Enh3 )
647 com2 = MAX( com2,1.e-20_r8 )
648 xph(i,k) = SQRT( com2 )
649 !-----------------------------------------------------------------
650 ! ... Add so4 effect
651 !-----------------------------------------------------------------
652 Eso4 = xso4(i,k)*xhnm(i,k) & ! /cm3(a)
653 *const0/xl
654 xph(i,k) = MIN( 1.e-2_r8,MAX( 1.e-7_r8,xph(i,k) + 2._r8*Eso4 ) )
655 if( iter > 1 ) then
656 delta = ABS( (xph(i,k) - delta)/delta )
657 converged = delta < .01_r8
658 if( converged ) then
659 exit
660 else
661 delta = xph(i,k)
662 end if
663 else
664 delta = xph(i,k)
665 end if
666 end do
667 if( .not. converged ) then
668 write(iulog,*) 'SETSOX: pH failed to converge @ (',i,',',k,'), % change=', &
669 100._r8*delta
670 #ifdef WRF_PORT
671 call wrf_message(iulog)
672 #endif
673 end if
674 else
675 xph(i,k) = 1.e-7_r8
676 end if
677 #if (defined MODAL_AERO)
678 end if ! if (cldfrc(i,k) >= 1.0e-5_r8)
679 #endif
680 end do
681 end do ! end pver loop for STEP 0
683 !==============================================================
684 ! ... Now use the actual PH
685 !==============================================================
686 do k = 1,pver
687 do i = 1,ncol
688 work1(i) = 1._r8 / tfld(i,k) - 1._r8 / 298._r8
689 tz = tfld(i,k)
690 #if (defined MODAL_AERO)
691 if (cldfrc(i,k) >= 1.0e-5_r8) then
692 xl = xlwc(i,k) / cldfrc(i,k)
693 #else
694 xl = xlwc(i,k)
695 #endif
696 patm = press(i,k)/101300._r8 ! press is in pascal
697 xam = press(i,k)/(1.38e-23_r8*tz) ! air density /M3
699 !-----------------------------------------------------------------
700 ! ... h2o2
701 !-----------------------------------------------------------------
702 xk = 7.4e4_r8 *EXP( 6621._r8*work1(i) )
703 xe = 2.2e-12_r8 *EXP(-3730._r8*work1(i) )
704 heh2o2(i,k) = xk*(1._r8 + xe/xph(i,k))
706 !-----------------------------------------------------------------
707 ! ... so2
708 !-----------------------------------------------------------------
709 xk = 1.23_r8 *EXP( 3120._r8*work1(i) )
710 xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) )
711 x2 = 6.0e-8_r8*EXP( 1120._r8*work1(i) )
713 wrk = xe/xph(i,k)
714 heso2(i,k) = xk*(1._r8 + wrk*(1._r8 + x2/xph(i,k)))
716 #if (defined MODAL_AERO)
717 !-----------------------------------------------------------------
718 ! ... nh3
719 !-----------------------------------------------------------------
720 xk = 58._r8 *EXP( 4085._r8*work1(i) )
721 xe = 1.7e-5_r8*EXP(-4325._r8*work1(i) )
722 henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/xkw)
723 #endif
725 !-----------------------------------------------------------------
726 ! ... o3
727 !-----------------------------------------------------------------
728 xk = 1.15e-2_r8 *EXP( 2560._r8*work1(i) )
729 heo3(i,k) = xk
731 !------------------------------------------------------------------------
732 ! ... for Ho2(g) -> H2o2(a) formation
733 ! schwartz JGR, 1984, 11589
734 !------------------------------------------------------------------------
735 kh4 = (kh2 + kh3*kh1/xph(i,k)) / ((1._r8 + kh1/xph(i,k))**2)
736 ho2s = kh0*xho2(i,k)*patm*(1._r8 + kh1/xph(i,k)) ! ho2s = ho2(a)+o2-
737 r1h2o2 = kh4*ho2s*ho2s ! prod(h2o2) in mole/L(w)/s
738 #if (defined MODAL_AERO)
739 r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s
740 /const0*1.e+6_r8 & ! correct a bug here
741 /xam
742 r2h2o2 = 0.0_r8
743 #else
744 r2h2o2 = r1h2o2*xlwc(i,k) & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s
745 *const0 & ! mole/fm3(a)/s * 1.e-3 = mole/cm3(a)/s
746 /xam ! /cm3(a)/s / air-den = mix-ratio/s
747 #endif
748 xh2o2(i,k) = xh2o2(i,k) + r2h2o2*dtime ! updated h2o2 by het production
750 !-----------------------------------------------
751 ! ... Partioning
752 !-----------------------------------------------
753 !------------------------------------------------------------------------
754 ! ... h2o2
755 !------------------------------------------------------------------------
756 px = heh2o2(i,k) * Ra * tz * xl
757 h2o2g = xh2o2(i,k)/(1._r8+ px)
759 !------------------------------------------------------------------------
760 ! ... so2
761 !------------------------------------------------------------------------
762 px = heso2(i,k) * Ra * tz * xl
763 so2g = xso2(i,k)/(1._r8+ px)
765 !------------------------------------------------------------------------
766 ! ... o3 ============
767 !------------------------------------------------------------------------
768 px = heo3(i,k) * Ra * tz * xl
769 o3g = xo3(i,k)/(1._r8+ px)
771 #if (defined MODAL_AERO)
772 !------------------------------------------------------------------------
773 ! ... nh3
774 !------------------------------------------------------------------------
775 px = henh3(i,k) * Ra * tz * xl
776 #if ( defined MODAL_AERO_7MODE )
777 nh3g = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px)
778 #elif ( defined MODAL_AERO_3MODE )
779 if(xl .ge. 1.e-8_r8) then
780 nh3g = xnh4(i,k)/px
781 endif
782 #endif
783 #endif
784 !-----------------------------------------------
785 ! ... Aqueous phase reaction rates
786 ! SO2 + H2O2 -> SO4
787 ! SO2 + O3 -> SO4
788 !-----------------------------------------------
790 !------------------------------------------------------------------------
791 ! ... S(IV) (HSO3) + H2O2
792 !------------------------------------------------------------------------
793 rah2o2 = 8.e4_r8 * EXP( -3650._r8*work1(i) ) &
794 / (.1_r8 + xph(i,k))
796 !------------------------------------------------------------------------
797 ! ... S(IV)+ O3
798 !------------------------------------------------------------------------
799 rao3 = 4.39e11_r8 * EXP(-4131._r8/tz) &
800 + 2.56e3_r8 * EXP(-996._r8 /tz) /xph(i,k)
802 !-----------------------------------------------------------------
803 ! ... Prediction after aqueous phase
804 ! so4
805 ! When Cloud is present
806 !
807 ! S(IV) + H2O2 = S(VI)
808 ! S(IV) + O3 = S(VI)
809 !
810 ! reference:
811 ! (1) Seinfeld
812 ! (2) Benkovitz
813 !-----------------------------------------------------------------
814 !............................
815 ! S(IV) + H2O2 = S(VI)
816 !............................
818 IF (XL .ge. 1.e-8_r8) THEN !! WHEN CLOUD IS PRESENTED
820 #if (defined MODAL_AERO)
821 pso4 = rah2o2 * 7.4e4_r8*EXP(6621._r8*work1(i)) * h2o2g * patm &
822 * 1.23_r8 *EXP(3120._r8*work1(i)) * so2g * patm
823 pso4 = pso4 & ! [M/s] = [mole/L(w)/s]
824 * xl & ! [mole/L(a)/s]
825 / const0 & ! [/L(a)/s]
826 / xhnm(i,k)
827 #else
828 pso4 = rah2o2 * heh2o2(i,k)*h2o2g &
829 * heso2(i,k) *so2g ! [M/s]
830 pso4 = pso4 & ! [M/s] = [mole/L(w)/s]
831 * xlwc(i,k) & ! [mole/L(a)/s]
832 / const0 & ! [/L(a)/s]
833 / xhnm(i,k)
834 #endif
836 ccc = pso4*dtime
837 ccc = max(ccc, 1.e-30_r8)
839 #if (defined MODAL_AERO)
840 xso4_init(i,k)=xso4(i,k)
841 #endif
842 IF (xh2o2(i,k) .gt. xso2(i,k)) THEN
843 if (ccc .gt. xso2(i,k)) then
844 xso4(i,k)=xso4(i,k)+xso2(i,k)
845 #if (defined MODAL_AERO)
846 xh2o2(i,k)=xh2o2(i,k)-xso2(i,k)
847 xso2(i,k)=1.e-20_r8
848 #else
849 xso2(i,k)=1.e-20_r8
850 xh2o2(i,k)=xh2o2(i,k)-xso2(i,k)
851 #endif
852 else
853 xso4(i,k) = xso4(i,k) + ccc
854 xh2o2(i,k) = xh2o2(i,k) - ccc
855 xso2(i,k) = xso2(i,k) - ccc
856 end if
858 ELSE
859 if (ccc .gt. xh2o2(i,k)) then
860 xso4(i,k)=xso4(i,k)+xh2o2(i,k)
861 xso2(i,k)=xso2(i,k)-xh2o2(i,k)
862 xh2o2(i,k)=1.e-20_r8
863 else
864 xso4(i,k) = xso4(i,k) + ccc
865 xh2o2(i,k) = xh2o2(i,k) - ccc
866 xso2(i,k) = xso2(i,k) - ccc
867 end if
868 END IF
870 #if (defined MODAL_AERO)
871 delso4_hprxn = xso4(i,k) - xso4_init(i,k)
872 #endif
873 !...........................
874 ! S(IV) + O3 = S(VI)
875 !...........................
877 #if (defined MODAL_AERO)
878 pso4 = rao3 * heo3(i,k)*o3g*patm * heso2(i,k)*so2g*patm ! [M/s]
879 pso4 = pso4 & ! [M/s] = [mole/L(w)/s]
880 * xl & ! [mole/L(a)/s]
881 / const0 & ! [/L(a)/s]
882 / xhnm(i,k) ! [mixing ratio/s]
883 #else
884 pso4 = rao3 * heo3(i,k)*o3g * heso2(i,k)*so2g ! [M/s]
885 pso4 = pso4 & ! [M/s] = [mole/L(w)/s]
886 * xlwc(i,k) & ! [mole/L(a)/s]
887 / const0 & ! [/L(a)/s]
888 / xhnm(i,k) ! [mixing ratio/s]
889 #endif
890 ccc = pso4*dtime
891 ccc = max(ccc, 1.e-30_r8)
893 #if (defined MODAL_AERO)
894 xso4_init(i,k)=xso4(i,k)
895 #endif
897 if (ccc .gt. xso2(i,k)) then
898 xso4(i,k)=xso4(i,k)+xso2(i,k)
899 xso2(i,k)=1.e-20_r8
900 else
901 xso4(i,k) = xso4(i,k) + ccc
902 xso2(i,k) = xso2(i,k) - ccc
903 end if
905 #if (defined MODAL_AERO)
906 delso4_o3rxn = xso4(i,k) - xso4_init(i,k)
907 #endif
909 #if (defined MODAL_AERO)
910 #if ( defined MODAL_AERO_7MODE )
911 delnh3 = nh3g - xnh3(i,k)
912 delnh4 = - delnh3
913 #endif
914 #endif
915 #if (defined MODAL_AERO)
916 !-------------------------------------------------------------------------
917 ! compute factors for partitioning aerosol mass gains among modes
918 ! the factors are proportional to the activated particle MR for each
919 ! mode, which is the MR of cloud drops "associated with" the mode
920 ! thus we are assuming the cloud drop size is independent of the
921 ! associated aerosol mode properties (i.e., drops associated with
922 ! Aitken and coarse sea-salt particles are same size)
923 !
924 ! qnum_c(n) = activated particle number MR for mode n (these are just
925 ! used for partitioning among modes, so don't need to divide by cldfrc)
927 do n = 1, ntot_amode
928 qnum_c(n) = 0.0
929 l = numptrcw_amode(n) - loffset
930 if (l > 0) qnum_c(n) = max( 0.0_r8, qcw(i,k,l) )
931 end do
933 ! force qnum_c(n) to be positive for n=modeptr_accum or n=1
934 n = modeptr_accum
935 if (n <= 0) n = 1
936 qnum_c(n) = max( 1.0e-10_r8, qnum_c(n) )
938 ! faqgain_so4(n) = fraction of total so4_c gain going to mode n
939 ! these are proportional to the activated particle MR for each mode
940 sumf = 0.0
941 do n = 1, ntot_amode
942 faqgain_so4(n) = 0.0
943 if (lptr_so4_cw_amode(n) > 0) then
944 faqgain_so4(n) = qnum_c(n)
945 sumf = sumf + faqgain_so4(n)
946 end if
947 end do
949 if (sumf > 0.0) then
950 do n = 1, ntot_amode
951 faqgain_so4(n) = faqgain_so4(n) / sumf
952 end do
953 end if
954 ! at this point (sumf <= 0.0) only when all the faqgain_so4 are zero
956 ! faqgain_msa(n) = fraction of total msa_c gain going to mode n
957 ntot_msa_c = 0
958 sumf = 0.0
959 do n = 1, ntot_amode
960 faqgain_msa(n) = 0.0
961 if (lptr_msa_cw_amode(n) > 0) then
962 faqgain_msa(n) = qnum_c(n)
963 ntot_msa_c = ntot_msa_c + 1
964 end if
965 sumf = sumf + faqgain_msa(n)
966 end do
968 if (sumf > 0.0) then
969 do n = 1, ntot_amode
970 faqgain_msa(n) = faqgain_msa(n) / sumf
971 end do
972 end if
973 ! at this point (sumf <= 0.0) only when all the faqgain_msa are zero
975 !-----------------------------------------------------------------------
976 ! compute uptake of h2so4 and msa to cloud water
977 !
978 ! first-order uptake rate is
979 ! 4*pi*(drop radius)*(drop number conc)
980 ! *(gas diffusivity)*(fuchs sutugin correction)
982 ! num_cd = (drop number conc in 1/cm^3)
983 num_cd = 1.0e-3_r8*cldnum(i,k)*cfact(i,k)/cldfrc(i,k)
984 num_cd = max( num_cd, 0.0_r8 )
986 ! rad_cd = (drop radius in cm), computed from liquid water and drop number,
987 ! then bounded by 0.5 and 50.0 micrometers
988 ! radxnum_cd = (drop radius)*(drop number conc)
989 ! volx34pi_cd = (3/4*pi) * (liquid water volume in cm^3/cm^3)
991 volx34pi_cd = xl*0.75_r8/pi
993 ! following holds because volx34pi_cd = num_cd*(rad_cd**3)
994 radxnum_cd = (volx34pi_cd*num_cd*num_cd)**0.3333333_r8
996 ! apply bounds to rad_cd to avoid the occasional unphysical value
997 if (radxnum_cd .le. volx34pi_cd*4.0e4_r8) then
998 radxnum_cd = volx34pi_cd*4.0e4_r8
999 rad_cd = 50.0e-4_r8
1000 else if (radxnum_cd .ge. volx34pi_cd*4.0e8_r8) then
1001 radxnum_cd = volx34pi_cd*4.0e8_r8
1002 rad_cd = 0.5e-4_r8
1003 else
1004 rad_cd = radxnum_cd/num_cd
1005 end if
1007 ! gasdiffus = h2so4 gas diffusivity from mosaic code (cm^2/s)
1008 ! (pmid must be Pa)
1009 gasdiffus = 0.557_r8 * (tfld(i,k)**1.75_r8) / press(i,k)
1011 ! gasspeed = h2so4 gas mean molecular speed from mosaic code (cm/s)
1012 gasspeed = 1.455e4_r8 * sqrt(tfld(i,k)/98.0_r8)
1014 ! knudsen number
1015 knudsen = 3.0_r8*gasdiffus/(gasspeed*rad_cd)
1017 ! following assumes accomodation coefficient = 0.65
1018 ! (Adams & Seinfeld, 2002, JGR, and references therein)
1019 ! fuchs_sutugin = (0.75*accom*(1. + knudsen)) /
1020 ! (knudsen*(1.0 + knudsen + 0.283*accom) + 0.75*accom)
1021 fuchs_sutugin = (0.4875_r8*(1. + knudsen)) / &
1022 (knudsen*(1.184_r8 + knudsen) + 0.4875_r8)
1024 ! instantaneous uptake rate
1025 uptkrate = 12.56637_r8*radxnum_cd*gasdiffus*fuchs_sutugin
1027 ! average uptake rate over dtime
1028 uptkrate = (1.0_r8 - exp(-min(100._r8,dtime*uptkrate))) / dtime
1030 ! dso4dt_gasuptk = so4_c tendency from h2so4 gas uptake (mol/mol/s)
1031 ! dmsadt_gasuptk = msa_c tendency from msa gas uptake (mol/mol/s)
1032 dso4dt_gasuptk = xh2so4(i,k) * uptkrate
1033 if (id_msa > 0) then
1034 dmsadt_gasuptk = xmsa(i,k) * uptkrate
1035 else
1036 dmsadt_gasuptk = 0.0_r8
1037 end if
1039 ! if no modes have msa aerosol, then "rename" scavenged msa gas to so4
1040 dmsadt_gasuptk_toso4 = 0.0_r8
1041 dmsadt_gasuptk_tomsa = dmsadt_gasuptk
1042 if (ntot_msa_c == 0) then
1043 dmsadt_gasuptk_tomsa = 0.0_r8
1044 dmsadt_gasuptk_toso4 = dmsadt_gasuptk
1045 end if
1047 !-----------------------------------------------------------------------
1048 ! now compute TMR tendencies
1049 ! this includes the above aqueous so2 chemistry AND
1050 ! the uptake of highly soluble aerosol precursor gases (h2so4, msa, ...)
1051 ! AND the wetremoval of dissolved, unreacted so2 and h2o2
1053 dso4dt_aqrxn = (delso4_o3rxn + delso4_hprxn) / dtime
1054 dso4dt_hprxn = delso4_hprxn / dtime
1056 ! fwetrem = fraction of in-cloud-water material that is wet removed
1057 ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*clwlrat(i,k)))) )
1058 fwetrem = 0.0 ! don't have so4 & msa wet removal here
1060 ! compute TMR tendencies for so4 and msa aerosol-in-cloud-water
1061 do n = 1, ntot_amode
1062 l = lptr_so4_cw_amode(n) - loffset
1063 if (l > 0) then
1064 dqdt_aqso4(i,k,l) = faqgain_so4(n)*dso4dt_aqrxn*cldfrc(i,k)
1065 dqdt_aqh2so4(i,k,l) = faqgain_so4(n)* &
1066 (dso4dt_gasuptk + dmsadt_gasuptk_toso4)*cldfrc(i,k)
1067 dqdt_aq = dqdt_aqso4(i,k,l) + dqdt_aqh2so4(i,k,l)
1068 dqdt_wr = -fwetrem*dqdt_aq
1069 dqdt= dqdt_aq + dqdt_wr
1070 qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
1071 end if
1073 l = lptr_msa_cw_amode(n) - loffset
1074 if (l > 0) then
1075 dqdt_aq = faqgain_msa(n)*dmsadt_gasuptk_tomsa*cldfrc(i,k)
1076 dqdt_wr = -fwetrem*dqdt_aq
1077 dqdt = dqdt_aq + dqdt_wr
1078 qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
1079 end if
1081 l = lptr_nh4_cw_amode(n) - loffset
1082 if (l > 0) then
1083 if (delnh4 > 0.0_r8) then
1084 dqdt_aq = faqgain_so4(n)*delnh4/dtime*cldfrc(i,k)
1085 dqdt = dqdt_aq
1086 qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
1087 else
1088 dqdt = (qcw(i,k,l)/max(xnh4c(i,k),1.0e-35_r8)) &
1089 *delnh4/dtime*cldfrc(i,k)
1090 qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
1091 endif
1092 end if
1093 end do
1095 ! For gas species, tendency includes
1096 ! reactive uptake to cloud water that essentially transforms the gas to
1097 ! a different species. Wet removal associated with this is applied
1098 ! to the "new" species (e.g., so4_c) rather than to the gas.
1099 ! wet removal of the unreacted gas that is dissolved in cloud water.
1100 ! Need to multiply both these parts by cldfrc
1102 ! h2so4 (g) & msa (g)
1103 qin(i,k,id_h2so4) = qin(i,k,id_h2so4) - dso4dt_gasuptk * dtime * cldfrc(i,k)
1104 if (id_msa > 0) qin(i,k,id_msa) = qin(i,k,id_msa) - dmsadt_gasuptk * dtime * cldfrc(i,k)
1106 ! frso2_g is fraction of total so2 as gas, etc.
1107 frso2_g = 1.0_r8 / ( 1.0_r8 + heso2(i,k) * Ra * tz * xl )
1108 frh2o2_g = 1.0_r8 / ( 1.0_r8 + heh2o2(i,k) * Ra * tz * xl )
1110 ! frso2_c is fraction of total so2 in cloud water, etc.
1111 frso2_c = max( 0.0_r8, (1.0_r8-frso2_g) )
1112 frh2o2_c = max( 0.0_r8, (1.0_r8-frh2o2_g) )
1114 ! so2 -- the first order loss rate for so2 is frso2_c*clwlrat(i,k)
1115 ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frso2_c*clwlrat(i,k)))) )
1116 fwetrem = 0.0 ! don't include so2 wet removal here
1118 dqdt_wr = -fwetrem*xso2(i,k)/dtime*cldfrc(i,k)
1119 dqdt_aq = -dso4dt_aqrxn*cldfrc(i,k)
1120 dqdt = dqdt_aq + dqdt_wr
1121 qin(i,k,id_so2) = qin(i,k,id_so2) + dqdt * dtime
1123 ! h2o2 -- the first order loss rate for h2o2 is frh2o2_c*clwlrat(i,k)
1124 ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frh2o2_c*clwlrat(i,k)))) )
1125 fwetrem = 0.0 ! don't include h2o2 wet removal here
1127 dqdt_wr = -fwetrem*xh2o2(i,k)/dtime*cldfrc(i,k)
1128 dqdt_aq = -dso4dt_hprxn*cldfrc(i,k)
1129 dqdt = dqdt_aq + dqdt_wr
1130 qin(i,k,id_h2o2) = qin(i,k,id_h2o2) + dqdt * dtime
1132 ! NH3
1133 #if ( defined MODAL_AERO_7MODE )
1134 dqdt_aq = delnh3/dtime*cldfrc(i,k)
1135 dqdt = dqdt_aq
1136 qin(i,k,id_nh3) = qin(i,k,id_nh3) + dqdt * dtime
1137 #endif
1139 ! for SO4 from H2O2/O3 budgets
1140 dqdt_aqhprxn(i,k) = dso4dt_hprxn*cldfrc(i,k)
1141 dqdt_aqo3rxn(i,k) = (dso4dt_aqrxn - dso4dt_hprxn)*cldfrc(i,k)
1142 #endif
1144 END IF !! WHEN CLOUD IS PRESENTED
1146 #if (defined MODAL_AERO)
1147 end if !! if (cldfrc(i,k) >= 1.0e-5_r8) then
1148 #endif
1150 end do
1151 end do
1153 !==============================================================
1154 ! ... Update the mixing ratios
1155 !==============================================================
1156 do k = 1,pver
1157 #if (defined MODAL_AERO)
1158 qin(:,k,id_so2) = MAX( qin(:,k,id_so2), small_value )
1159 qcw(:,k,id_so4_1a) = MAX( qcw(:,k,id_so4_1a), small_value )
1160 qcw(:,k,id_so4_2a) = MAX( qcw(:,k,id_so4_2a), small_value )
1161 qcw(:,k,id_so4_3a) = MAX( qcw(:,k,id_so4_3a), small_value )
1162 #if ( defined MODAL_AERO_7MODE )
1163 qcw(:,k,id_so4_4a) = MAX( qcw(:,k,id_so4_4a), small_value )
1164 qcw(:,k,id_so4_5a) = MAX( qcw(:,k,id_so4_5a), small_value )
1165 qcw(:,k,id_so4_6a) = MAX( qcw(:,k,id_so4_6a), small_value )
1166 #endif
1168 #if ( defined MODAL_AERO_7MODE )
1169 qin(:,k,id_nh3) = MAX( qin(:,k,id_nh3), small_value )
1170 qcw(:,k,id_nh4_1a) = MAX( qcw(:,k,id_nh4_1a), small_value )
1171 qcw(:,k,id_nh4_2a) = MAX( qcw(:,k,id_nh4_2a), small_value )
1172 qcw(:,k,id_nh4_3a) = MAX( qcw(:,k,id_nh4_3a), small_value )
1173 qcw(:,k,id_nh4_4a) = MAX( qcw(:,k,id_nh4_4a), small_value )
1174 qcw(:,k,id_nh4_5a) = MAX( qcw(:,k,id_nh4_5a), small_value )
1175 qcw(:,k,id_nh4_6a) = MAX( qcw(:,k,id_nh4_6a), small_value )
1176 #endif
1177 #else
1178 if (.not. inv_so2) then
1179 qin(:,k,id_so2) = MAX( xso2(:,k), small_value )
1180 endif
1181 if (.not. inv_h2o2) then
1182 qin(:,k,id_h2o2)= MAX( xh2o2(:,k), small_value )
1183 endif
1184 qin(:,k,id_so4) = MAX( xso4(:,k), small_value )
1185 #endif
1186 end do
1188 #if (defined MODAL_AERO)
1189 do n = 1, ntot_amode
1190 m = lptr_so4_cw_amode(n)
1191 l = m - loffset
1192 if (l > 0) then
1193 sflx(:)=0._r8
1194 do k=1,pver
1195 do i=1,ncol
1196 sflx(i)=sflx(i)+dqdt_aqso4(i,k,l)*adv_mass(l)/mbar(i,k) &
1197 *pdel(i,k)/gravit ! kg/m2/s
1198 enddo
1199 enddo
1200 call outfld( trim(cnst_name_cw(m))//'AQSO4', sflx(:ncol), ncol, lchnk)
1202 sflx(:)=0._r8
1203 do k=1,pver
1204 do i=1,ncol
1205 sflx(i)=sflx(i)+dqdt_aqh2so4(i,k,l)*adv_mass(l)/mbar(i,k) &
1206 *pdel(i,k)/gravit ! kg/m2/s
1207 enddo
1208 enddo
1209 call outfld( trim(cnst_name_cw(m))//'AQH2SO4', sflx(:ncol), ncol, lchnk)
1210 endif
1211 end do
1213 sflx(:)=0._r8
1214 do k=1,pver
1215 do i=1,ncol
1216 sflx(i)=sflx(i)+dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) &
1217 *pdel(i,k)/gravit ! kg SO4 /m2/s
1218 enddo
1219 enddo
1220 call outfld( 'AQSO4_H2O2', sflx(:ncol), ncol, lchnk)
1221 sflx(:)=0._r8
1222 do k=1,pver
1223 do i=1,ncol
1224 sflx(i)=sflx(i)+dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) &
1225 *pdel(i,k)/gravit ! kg SO4 /m2/s
1226 enddo
1227 enddo
1228 call outfld( 'AQSO4_O3', sflx(:ncol), ncol, lchnk)
1230 xphlwc(:,:) = 0._r8
1231 do k = 1,pver
1232 do i = 1,ncol
1233 if (cldfrc(i,k) >= 1.0e-5_r8) then
1234 if( lwc(i,k) >= 1.0e-8_r8 ) then
1235 xphlwc(i,k) = -1._r8*log10(xph(i,k)) * lwc(i,k)
1236 endif
1237 endif
1238 enddo
1239 enddo
1240 call outfld( 'XPH_LWC', xphlwc(:ncol,:), ncol ,lchnk )
1241 #endif
1243 end subroutine SETSOX
1245 end module MO_SETSOX