Line data Source code
1 :
2 : module MO_SETSOX
3 :
4 : use shr_kind_mod, only : r8 => shr_kind_r8
5 : use cam_logfile, only : iulog
6 :
7 : private
8 : public :: sox_inti, setsox
9 : public :: has_sox
10 :
11 : save
12 : logical :: inv_o3
13 : integer :: id_msa
14 :
15 : integer :: id_so2, id_nh3, id_hno3, id_h2o2, id_o3, id_ho2
16 : integer :: id_so4, id_h2so4
17 :
18 : logical :: has_sox = .true.
19 : logical :: inv_so2, inv_nh3, inv_hno3, inv_h2o2, inv_ox, inv_nh4no3, inv_ho2
20 :
21 : logical :: cloud_borne = .false.
22 : logical :: modal_aerosols = .false.
23 :
24 : contains
25 :
26 : !-----------------------------------------------------------------------
27 : !-----------------------------------------------------------------------
28 1490712 : subroutine sox_inti
29 : !-----------------------------------------------------------------------
30 : ! ... initialize the hetero sox routine
31 : !-----------------------------------------------------------------------
32 :
33 : use mo_chem_utls, only : get_spc_ndx, get_inv_ndx
34 : use spmd_utils, only : masterproc
35 : use phys_control, only : phys_getopts
36 : use sox_cldaero_mod, only : sox_cldaero_init
37 :
38 : implicit none
39 :
40 :
41 : call phys_getopts( &
42 1536 : prog_modal_aero_out=modal_aerosols )
43 :
44 1536 : cloud_borne = modal_aerosols
45 :
46 : !-----------------------------------------------------------------
47 : ! ... get species indicies
48 : !-----------------------------------------------------------------
49 :
50 1536 : if (cloud_borne) then
51 1536 : id_h2so4 = get_spc_ndx( 'H2SO4' )
52 : else
53 0 : id_so4 = get_spc_ndx( 'SO4' )
54 : endif
55 1536 : id_msa = get_spc_ndx( 'MSA' )
56 :
57 1536 : inv_so2 = .false.
58 1536 : id_so2 = get_inv_ndx( 'SO2' )
59 1536 : inv_so2 = id_so2 > 0
60 1536 : if ( .not. inv_so2 ) then
61 1536 : id_so2 = get_spc_ndx( 'SO2' )
62 : endif
63 :
64 1536 : inv_NH3 = .false.
65 1536 : id_NH3 = get_inv_ndx( 'NH3' )
66 1536 : inv_NH3 = id_NH3 > 0
67 1536 : if ( .not. inv_NH3 ) then
68 1536 : id_NH3 = get_spc_ndx( 'NH3' )
69 : endif
70 :
71 1536 : inv_HNO3 = .false.
72 1536 : id_HNO3 = get_inv_ndx( 'HNO3' )
73 1536 : inv_HNO3 = id_hno3 > 0
74 1536 : if ( .not. inv_HNO3 ) then
75 1536 : id_HNO3 = get_spc_ndx( 'HNO3' )
76 : endif
77 :
78 1536 : inv_H2O2 = .false.
79 1536 : id_H2O2 = get_inv_ndx( 'H2O2' )
80 1536 : inv_H2O2 = id_H2O2 > 0
81 1536 : if ( .not. inv_H2O2 ) then
82 1536 : id_H2O2 = get_spc_ndx( 'H2O2' )
83 : endif
84 :
85 1536 : inv_HO2 = .false.
86 1536 : id_HO2 = get_inv_ndx( 'HO2' )
87 1536 : inv_HO2 = id_HO2 > 0
88 1536 : if ( .not. inv_HO2 ) then
89 0 : id_HO2 = get_spc_ndx( 'HO2' )
90 : endif
91 :
92 1536 : inv_o3 = get_inv_ndx( 'O3' ) > 0
93 1536 : if (inv_o3) then
94 1536 : id_o3 = get_inv_ndx( 'O3' )
95 : else
96 0 : id_o3 = get_spc_ndx( 'O3' )
97 : endif
98 1536 : inv_ho2 = get_inv_ndx( 'HO2' ) > 0
99 1536 : if (inv_ho2) then
100 1536 : id_ho2 = get_inv_ndx( 'HO2' )
101 : else
102 0 : id_ho2 = get_spc_ndx( 'HO2' )
103 : endif
104 :
105 1536 : has_sox = (id_so2>0) .and. (id_h2o2>0) .and. (id_o3>0) .and. (id_ho2>0)
106 1536 : if (cloud_borne) then
107 1536 : has_sox = has_sox .and. (id_h2so4>0)
108 : else
109 0 : has_sox = has_sox .and. (id_so4>0) .and. (id_nh3>0)
110 : endif
111 :
112 1536 : if (masterproc) then
113 2 : write(iulog,*) 'sox_inti: has_sox = ',has_sox
114 : endif
115 :
116 1536 : if( has_sox ) then
117 1536 : if (masterproc) then
118 2 : write(iulog,*) '-----------------------------------------'
119 2 : write(iulog,*) 'mozart will do sox aerosols'
120 2 : write(iulog,*) '-----------------------------------------'
121 : endif
122 : else
123 : return
124 : end if
125 :
126 1536 : call sox_cldaero_init()
127 :
128 1536 : end subroutine sox_inti
129 :
130 : !-----------------------------------------------------------------------
131 : !-----------------------------------------------------------------------
132 1489176 : subroutine SETSOX( &
133 : ncol, &
134 : lchnk, &
135 : loffset,&
136 : dtime, &
137 2978352 : press, &
138 2978352 : pdel, &
139 1489176 : tfld, &
140 1489176 : mbar, &
141 2978352 : lwc, &
142 2978352 : cldfrc, &
143 1489176 : cldnum, &
144 1489176 : xhnm, &
145 1489176 : invariants, &
146 1489176 : qcw, &
147 1489176 : qin, &
148 1489176 : xphlwc, &
149 1489176 : aqso4, &
150 1489176 : aqh2so4,&
151 1489176 : aqso4_h2o2, &
152 1489176 : aqso4_o3, &
153 : yph_in, &
154 1489176 : aqso4_h2o2_3d, &
155 1489176 : aqso4_o3_3d &
156 : )
157 :
158 : !-----------------------------------------------------------------------
159 : ! ... Compute heterogeneous reactions of SOX
160 : !
161 : ! (0) using initial PH to calculate PH
162 : ! (a) HENRYs law constants
163 : ! (b) PARTIONING
164 : ! (c) PH values
165 : !
166 : ! (1) using new PH to repeat
167 : ! (a) HENRYs law constants
168 : ! (b) PARTIONING
169 : ! (c) REACTION rates
170 : ! (d) PREDICTION
171 : !-----------------------------------------------------------------------
172 : !
173 1536 : use ppgrid, only : pcols, pver
174 : use chem_mods, only : gas_pcnst, nfs
175 : use chem_mods, only : adv_mass
176 : use physconst, only : mwdry, gravit
177 : use mo_constants, only : pi
178 : use sox_cldaero_mod, only : sox_cldaero_update, sox_cldaero_create_obj, sox_cldaero_destroy_obj
179 : use cldaero_mod, only : cldaero_conc_t
180 :
181 : !
182 : implicit none
183 : !
184 : !-----------------------------------------------------------------------
185 : ! ... Dummy arguments
186 : !-----------------------------------------------------------------------
187 : integer, intent(in) :: ncol ! num of columns in chunk
188 : integer, intent(in) :: lchnk ! chunk id
189 : integer, intent(in) :: loffset ! offset of chem tracers in the advected tracers array
190 : real(r8), intent(in) :: dtime ! time step (sec)
191 : real(r8), intent(in) :: press(:,:) ! midpoint pressure ( Pa )
192 : real(r8), intent(in) :: pdel(:,:) ! pressure thickness of levels (Pa)
193 : real(r8), intent(in) :: tfld(:,:) ! temperature
194 : real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu )
195 : real(r8), target, intent(in) :: lwc(:,:) ! cloud liquid water content (kg/kg)
196 : real(r8), target, intent(in) :: cldfrc(:,:) ! cloud fraction
197 : real(r8), intent(in) :: cldnum(:,:) ! droplet number concentration (#/kg)
198 : real(r8), intent(in) :: xhnm(:,:) ! total atms density ( /cm**3)
199 : real(r8), intent(in) :: invariants(:,:,:)
200 : real(r8), target, intent(inout) :: qcw(:,:,:) ! cloud-borne aerosol (vmr)
201 : real(r8), intent(inout) :: qin(:,:,:) ! transported species ( vmr )
202 : real(r8), intent(out) :: xphlwc(:,:) ! pH value multiplied by lwc
203 :
204 : real(r8), intent(out) :: aqso4(:,:) ! aqueous phase chemistry
205 : real(r8), intent(out) :: aqh2so4(:,:) ! aqueous phase chemistry
206 : real(r8), intent(out) :: aqso4_h2o2(:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2)
207 : real(r8), intent(out) :: aqso4_o3(:) ! SO4 aqueous phase chemistry due to O3 (kg/m2)
208 : real(r8), intent(in), optional :: yph_in ! ph value
209 : real(r8), intent(out), optional :: aqso4_h2o2_3d(:, :) ! 3D SO4 aqueous phase chemistry due to H2O2 (kg/m2)
210 : real(r8), intent(out), optional :: aqso4_o3_3d(:, :) ! 3D SO4 aqueous phase chemistry due to O3 (kg/m2)
211 :
212 :
213 : !-----------------------------------------------------------------------
214 : ! ... Local variables
215 : !
216 : ! xhno3 ... in mixing ratio
217 : !-----------------------------------------------------------------------
218 : integer, parameter :: itermax = 20
219 : real(r8), parameter :: ph0 = 5.0_r8 ! INITIAL PH VALUES
220 : real(r8), parameter :: const0 = 1.e3_r8/6.023e23_r8
221 : real(r8), parameter :: xa0 = 11._r8
222 : real(r8), parameter :: xb0 = -.1_r8
223 : real(r8), parameter :: xa1 = 1.053_r8
224 : real(r8), parameter :: xb1 = -4.368_r8
225 : real(r8), parameter :: xa2 = 1.016_r8
226 : real(r8), parameter :: xb2 = -2.54_r8
227 : real(r8), parameter :: xa3 = .816e-32_r8
228 : real(r8), parameter :: xb3 = .259_r8
229 :
230 : real(r8), parameter :: kh0 = 9.e3_r8 ! HO2(g) -> Ho2(a)
231 : real(r8), parameter :: kh1 = 2.05e-5_r8 ! HO2(a) -> H+ + O2-
232 : real(r8), parameter :: kh2 = 8.6e5_r8 ! HO2(a) + ho2(a) -> h2o2(a) + o2
233 : real(r8), parameter :: kh3 = 1.e8_r8 ! HO2(a) + o2- -> h2o2(a) + o2
234 : real(r8), parameter :: Ra = 8314._r8/101325._r8 ! universal constant (atm)/(M-K)
235 : real(r8), parameter :: xkw = 1.e-14_r8 ! water acidity
236 :
237 : !
238 2978352 : real(r8) :: xdelso4hp(ncol,pver)
239 :
240 : integer :: k, i, iter, file
241 : real(r8) :: wrk, delta
242 : real(r8) :: xph0, aden, xk, xe, x2
243 : real(r8) :: tz, xl, px, qz, pz, es, qs, patm
244 : real(r8) :: Eso2, Eso4, Ehno3, Eco2, Eh2o, Enh3
245 : real(r8) :: so2g, h2o2g, co2g, o3g
246 : real(r8) :: hno3a, nh3a, so2a, h2o2a, co2a, o3a
247 : real(r8) :: rah2o2, rao3, pso4, ccc
248 : real(r8) :: cnh3, chno3, com, com1, com2, xra
249 :
250 2978352 : real(r8) :: hno3g(ncol,pver), nh3g(ncol,pver)
251 : !
252 : !-----------------------------------------------------------------------
253 : ! for Ho2(g) -> H2o2(a) formation
254 : ! schwartz JGR, 1984, 11589
255 : !-----------------------------------------------------------------------
256 : real(r8) :: kh4 ! kh2+kh3
257 : real(r8) :: xam ! air density /cm3
258 : real(r8) :: ho2s ! ho2s = ho2(a)+o2-
259 : real(r8) :: r1h2o2 ! prod(h2o2) by ho2 in mole/L(w)/s
260 : real(r8) :: r2h2o2 ! prod(h2o2) by ho2 in mix/s
261 :
262 : real(r8), dimension(ncol,pver) :: &
263 2978352 : xhno3, xh2o2, xso2, xso4, xno3, &
264 2978352 : xnh3, xnh4, xo3, &
265 2978352 : cfact, &
266 2978352 : xph, xho2, &
267 2978352 : xh2so4, xmsa, xso4_init, &
268 2978352 : hehno3, & ! henry law const for hno3
269 2978352 : heh2o2, & ! henry law const for h2o2
270 2978352 : heso2, & ! henry law const for so2
271 2978352 : henh3, & ! henry law const for nh3
272 2978352 : heo3 !!, & ! henry law const for o3
273 :
274 : real(r8) :: patm_x
275 :
276 2978352 : real(r8), dimension(ncol) :: work1
277 : logical :: converged
278 :
279 1489176 : real(r8), pointer :: xso4c(:,:)
280 1489176 : real(r8), pointer :: xnh4c(:,:)
281 1489176 : real(r8), pointer :: xno3c(:,:)
282 : type(cldaero_conc_t), pointer :: cldconc
283 :
284 : real(r8) :: fact1_hno3, fact2_hno3, fact3_hno3
285 : real(r8) :: fact1_so2, fact2_so2, fact3_so2, fact4_so2
286 : real(r8) :: fact1_nh3, fact2_nh3, fact3_nh3
287 : real(r8) :: tmp_hp, tmp_hso3, tmp_hco3, tmp_nh4, tmp_no3
288 : real(r8) :: tmp_oh, tmp_so3, tmp_so4
289 : real(r8) :: tmp_neg, tmp_pos
290 : real(r8) :: yph, yph_lo, yph_hi
291 : real(r8) :: ynetpos, ynetpos_lo, ynetpos_hi
292 :
293 : !-----------------------------------------------------------------
294 : ! ... NOTE: The press array is in pascals and must be
295 : ! mutiplied by 10 to yield dynes/cm**2.
296 : !-----------------------------------------------------------------
297 : !==================================================================
298 : ! ... First set the PH
299 : !==================================================================
300 : ! ... Initial values
301 : ! The values of so2, so4 are after (1) SLT, and CHEM
302 : !-----------------------------------------------------------------
303 1489176 : xph0 = 10._r8**(-ph0) ! initial PH value
304 :
305 139982544 : do k = 1,pver
306 138493368 : cfact(:,k) = xhnm(:,k) & ! /cm3(a)
307 : * 1.e6_r8 & ! /m3(a)
308 : * 1.38e-23_r8/287._r8 & ! Kg(a)/m3(a)
309 2452499712 : * 1.e-3_r8 ! Kg(a)/L(a)
310 : end do
311 :
312 1489176 : cldconc => sox_cldaero_create_obj( cldfrc,qcw,lwc, cfact, ncol, loffset )
313 1489176 : xso4c => cldconc%so4c
314 1489176 : xnh4c => cldconc%nh4c
315 1489176 : xno3c => cldconc%no3c
316 :
317 2314006344 : xso4(:,:) = 0._r8
318 2314006344 : xno3(:,:) = 0._r8
319 2314006344 : xnh4(:,:) = 0._r8
320 :
321 139982544 : do k = 1,pver
322 2312517168 : xph(:,k) = xph0 ! initial PH value
323 :
324 138493368 : if ( inv_so2 ) then
325 0 : xso2 (:,k) = invariants(:,k,id_so2)/xhnm(:,k) ! mixing ratio
326 : else
327 2312517168 : xso2 (:,k) = qin(:,k,id_so2) ! mixing ratio
328 : endif
329 :
330 138493368 : if (id_hno3 > 0) then
331 0 : xhno3(:,k) = qin(:,k,id_hno3)
332 : else
333 2312517168 : xhno3(:,k) = 0.0_r8
334 : endif
335 :
336 138493368 : if ( inv_h2o2 ) then
337 0 : xh2o2 (:,k) = invariants(:,k,id_h2o2)/xhnm(:,k) ! mixing ratio
338 : else
339 2312517168 : xh2o2 (:,k) = qin(:,k,id_h2o2) ! mixing ratio
340 : endif
341 :
342 138493368 : if (id_nh3 > 0) then
343 0 : xnh3 (:,k) = qin(:,k,id_nh3)
344 : else
345 2312517168 : xnh3 (:,k) = 0.0_r8
346 : endif
347 :
348 138493368 : if ( inv_o3 ) then
349 2312517168 : xo3 (:,k) = invariants(:,k,id_o3)/xhnm(:,k) ! mixing ratio
350 : else
351 0 : xo3 (:,k) = qin(:,k,id_o3) ! mixing ratio
352 : endif
353 138493368 : if ( inv_ho2 ) then
354 2312517168 : xho2 (:,k) = invariants(:,k,id_ho2)/xhnm(:,k)! mixing ratio
355 : else
356 0 : xho2 (:,k) = qin(:,k,id_ho2) ! mixing ratio
357 : endif
358 :
359 138493368 : if (cloud_borne) then
360 2312517168 : xh2so4(:,k) = qin(:,k,id_h2so4)
361 : else
362 0 : xso4 (:,k) = qin(:,k,id_so4) ! mixing ratio
363 : endif
364 139982544 : if (id_msa > 0) xmsa (:,k) = qin(:,k,id_msa)
365 :
366 : end do
367 :
368 : !-----------------------------------------------------------------
369 : ! ... Temperature dependent Henry constants
370 : !-----------------------------------------------------------------
371 139982544 : ver_loop0: do k = 1,pver !! pver loop for STEP 0
372 2314006344 : col_loop0: do i = 1,ncol
373 :
374 2174023800 : if (cloud_borne .and. cldfrc(i,k)>0._r8) then
375 464987798 : xso4(i,k) = xso4c(i,k) / cldfrc(i,k)
376 464987798 : xnh4(i,k) = xnh4c(i,k) / cldfrc(i,k)
377 464987798 : xno3(i,k) = xno3c(i,k) / cldfrc(i,k)
378 : endif
379 2174023800 : xl = cldconc%xlwc(i,k)
380 :
381 2312517168 : if( xl >= 1.e-8_r8 ) then
382 91242790 : work1(i) = 1._r8 / tfld(i,k) - 1._r8 / 298._r8
383 :
384 : !-----------------------------------------------------------------
385 : ! 21-mar-2011 changes by rce
386 : ! ph calculation now uses bisection method to solve the electro-neutrality equation
387 : ! 3-mode aerosols (where so4 is assumed to be nh4hso4)
388 : ! old code set xnh4c = so4c
389 : ! new code sets xnh4c = 0, then uses a -1 charge (instead of -2)
390 : ! for so4 when solving the electro-neutrality equation
391 : !-----------------------------------------------------------------
392 :
393 : !-----------------------------------------------------------------
394 : ! calculations done before iterating
395 : !-----------------------------------------------------------------
396 :
397 : !-----------------------------------------------------------------
398 91242790 : pz = .01_r8*press(i,k) !! pressure in mb
399 91242790 : tz = tfld(i,k)
400 91242790 : patm = pz/1013._r8
401 91242790 : xam = press(i,k)/(1.38e-23_r8*tz) !air density /M3
402 :
403 : !-----------------------------------------------------------------
404 : ! ... hno3
405 : !-----------------------------------------------------------------
406 : ! previous code
407 : ! hehno3(i,k) = xk*(1._r8 + xe/xph(i,k))
408 : ! px = hehno3(i,k) * Ra * tz * xl
409 : ! hno3g = xhno3(i,k)/(1._r8 + px)
410 : ! Ehno3 = xk*xe*hno3g *patm
411 : ! equivalent new code
412 : ! hehno3 = xk + xk*xe/hplus
413 : ! hno3g = xhno3/(1 + px)
414 : ! = xhno3/(1 + hehno3*ra*tz*xl)
415 : ! = xhno3/(1 + xk*ra*tz*xl*(1 + xe/hplus)
416 : ! ehno3 = hno3g*xk*xe*patm
417 : ! = xk*xe*patm*xhno3/(1 + xk*ra*tz*xl*(1 + xe/hplus)
418 : ! = ( fact1_hno3 )/(1 + fact2_hno3 *(1 + fact3_hno3/hplus)
419 : ! [hno3-] = ehno3/hplus
420 91242790 : xk = 2.1e5_r8 *EXP( 8700._r8*work1(i) )
421 91242790 : xe = 15.4_r8
422 91242790 : fact1_hno3 = xk*xe*patm*xhno3(i,k)
423 91242790 : fact2_hno3 = xk*ra*tz*xl
424 91242790 : fact3_hno3 = xe
425 :
426 : !-----------------------------------------------------------------
427 : ! ... so2
428 : !-----------------------------------------------------------------
429 : ! previous code
430 : ! heso2(i,k) = xk*(1._r8 + wrk*(1._r8 + x2/xph(i,k)))
431 : ! px = heso2(i,k) * Ra * tz * xl
432 : ! so2g = xso2(i,k)/(1._r8+ px)
433 : ! Eso2 = xk*xe*so2g *patm
434 : ! equivalent new code
435 : ! heso2 = xk + xk*xe/hplus * xk*xe*x2/hplus**2
436 : ! so2g = xso2/(1 + px)
437 : ! = xso2/(1 + heso2*ra*tz*xl)
438 : ! = xso2/(1 + xk*ra*tz*xl*(1 + (xe/hplus)*(1 + x2/hplus))
439 : ! eso2 = so2g*xk*xe*patm
440 : ! = xk*xe*patm*xso2/(1 + xk*ra*tz*xl*(1 + (xe/hplus)*(1 + x2/hplus))
441 : ! = ( fact1_so2 )/(1 + fact2_so2 *(1 + (fact3_so2/hplus)*(1 + fact4_so2/hplus)
442 : ! [hso3-] + 2*[so3--] = (eso2/hplus)*(1 + 2*x2/hplus)
443 91242790 : xk = 1.23_r8 *EXP( 3120._r8*work1(i) )
444 91242790 : xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) )
445 91242790 : x2 = 6.0e-8_r8*EXP( 1120._r8*work1(i) )
446 91242790 : fact1_so2 = xk*xe*patm*xso2(i,k)
447 91242790 : fact2_so2 = xk*ra*tz*xl
448 91242790 : fact3_so2 = xe
449 91242790 : fact4_so2 = x2
450 :
451 : !-----------------------------------------------------------------
452 : ! ... nh3
453 : !-----------------------------------------------------------------
454 : ! previous code
455 : ! henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/xkw)
456 : ! px = henh3(i,k) * Ra * tz * xl
457 : ! nh3g = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px)
458 : ! Enh3 = xk*xe*nh3g/xkw *patm
459 : ! equivalent new code
460 : ! henh3 = xk + xk*xe*hplus/xkw
461 : ! nh3g = xnh34/(1 + px)
462 : ! = xnh34/(1 + henh3*ra*tz*xl)
463 : ! = xnh34/(1 + xk*ra*tz*xl*(1 + xe*hplus/xkw)
464 : ! enh3 = nh3g*xk*xe*patm/xkw
465 : ! = ((xk*xe*patm/xkw)*xnh34)/(1 + xk*ra*tz*xl*(1 + xe*hplus/xkw)
466 : ! = ( fact1_nh3 )/(1 + fact2_nh3 *(1 + fact3_nh3*hplus)
467 : ! [nh4+] = enh3*hplus
468 91242790 : xk = 58._r8 *EXP( 4085._r8*work1(i) )
469 91242790 : xe = 1.7e-5_r8*EXP( -4325._r8*work1(i) )
470 :
471 91242790 : fact1_nh3 = (xk*xe*patm/xkw)*(xnh3(i,k)+xnh4(i,k))
472 91242790 : fact2_nh3 = xk*ra*tz*xl
473 91242790 : fact3_nh3 = xe/xkw
474 :
475 : !-----------------------------------------------------------------
476 : ! ... h2o effects
477 : !-----------------------------------------------------------------
478 91242790 : Eh2o = xkw
479 :
480 : !-----------------------------------------------------------------
481 : ! ... co2 effects
482 : !-----------------------------------------------------------------
483 91242790 : co2g = 330.e-6_r8 !330 ppm = 330.e-6 atm
484 91242790 : xk = 3.1e-2_r8*EXP( 2423._r8*work1(i) )
485 91242790 : xe = 4.3e-7_r8*EXP(-913._r8 *work1(i) )
486 91242790 : Eco2 = xk*xe*co2g *patm
487 :
488 : !-----------------------------------------------------------------
489 : ! ... so4 effect
490 : !-----------------------------------------------------------------
491 91242790 : Eso4 = xso4(i,k)*xhnm(i,k) & ! /cm3(a)
492 91242790 : *const0/xl
493 :
494 :
495 : !-----------------------------------------------------------------
496 : ! now use bisection method to solve electro-neutrality equation
497 : !
498 : ! during the iteration loop,
499 : ! yph_lo = lower ph value that brackets the root (i.e., correct ph)
500 : ! yph_hi = upper ph value that brackets the root (i.e., correct ph)
501 : ! yph = current ph value
502 : ! yposnet_lo and yposnet_hi = net positive ions for
503 : ! yph_lo and yph_hi
504 : !-----------------------------------------------------------------
505 1094912699 : do iter = 1,itermax
506 :
507 1094912699 : if (.not. present(yph_in)) then
508 1094912699 : if (iter == 1) then
509 : ! 1st iteration ph = lower bound value
510 : yph_lo = 2.0_r8
511 : yph_hi = yph_lo
512 : yph = yph_lo
513 1003669909 : else if (iter == 2) then
514 : ! 2nd iteration ph = upper bound value
515 : yph_hi = 7.0_r8
516 : yph = yph_hi
517 : else
518 : ! later iteration ph = mean of the two bracketing values
519 912427190 : yph = 0.5_r8*(yph_lo + yph_hi)
520 : end if
521 : else
522 0 : yph = yph_in
523 : end if
524 :
525 : ! calc current [H+] from ph
526 1094912699 : xph(i,k) = 10.0_r8**(-yph)
527 :
528 :
529 : !-----------------------------------------------------------------
530 : ! ... hno3
531 : !-----------------------------------------------------------------
532 1094912699 : Ehno3 = fact1_hno3/(1.0_r8 + fact2_hno3*(1.0_r8 + fact3_hno3/xph(i,k)))
533 :
534 : !-----------------------------------------------------------------
535 : ! ... so2
536 : !-----------------------------------------------------------------
537 : Eso2 = fact1_so2/(1.0_r8 + fact2_so2*(1.0_r8 + (fact3_so2/xph(i,k)) &
538 1094912699 : *(1.0_r8 + fact4_so2/xph(i,k))))
539 :
540 : !-----------------------------------------------------------------
541 : ! ... nh3
542 : !-----------------------------------------------------------------
543 1094912699 : Enh3 = fact1_nh3/(1.0_r8 + fact2_nh3*(1.0_r8 + fact3_nh3*xph(i,k)))
544 :
545 1094912699 : tmp_nh4 = Enh3 * xph(i,k)
546 1094912699 : tmp_hso3 = Eso2 / xph(i,k)
547 1094912699 : tmp_so3 = tmp_hso3 * 2.0_r8*fact4_so2/xph(i,k)
548 1094912699 : tmp_hco3 = Eco2 / xph(i,k)
549 1094912699 : tmp_oh = Eh2o / xph(i,k)
550 1094912699 : tmp_no3 = Ehno3 / xph(i,k)
551 1094912699 : tmp_so4 = cldconc%so4_fact*Eso4
552 1094912699 : tmp_pos = xph(i,k) + tmp_nh4
553 1094912699 : tmp_neg = tmp_oh + tmp_hco3 + tmp_no3 + tmp_hso3 + tmp_so3 + tmp_so4
554 :
555 1094912699 : ynetpos = tmp_pos - tmp_neg
556 :
557 :
558 : ! yposnet = net positive ions/charge
559 : ! if the correct ph is bracketed by yph_lo and yph_hi (with yph_lo < yph_hi),
560 : ! then you will have yposnet_lo > 0 and yposnet_hi < 0
561 1094912699 : converged = .false.
562 1094912699 : if (iter > 2) then
563 912427190 : if (ynetpos == 0.0_r8) then
564 : ! the exact solution was found (very unlikely)
565 91242790 : tmp_hp = xph(i,k)
566 : converged = .true.
567 : exit
568 912427190 : else if (ynetpos >= 0.0_r8) then
569 : ! net positive ions are >= 0 for both yph and yph_lo
570 : ! so replace yph_lo with yph
571 : yph_lo = yph
572 : ynetpos_lo = ynetpos
573 : else
574 : ! net positive ions are <= 0 for both yph and yph_hi
575 : ! so replace yph_hi with yph
576 441423040 : yph_hi = yph
577 441423040 : ynetpos_hi = ynetpos
578 : end if
579 :
580 912427190 : if (abs(yph_hi - yph_lo) .le. 0.005_r8) then
581 : ! |yph_hi - yph_lo| <= convergence criterion, so set
582 : ! final ph to their midpoint and exit
583 : ! (.005 absolute error in pH gives .01 relative error in H+)
584 91242719 : tmp_hp = xph(i,k)
585 91242719 : yph = 0.5_r8*(yph_hi + yph_lo)
586 91242719 : xph(i,k) = 10.0_r8**(-yph)
587 91242719 : converged = .true.
588 91242719 : exit
589 : else
590 : ! do another iteration
591 : converged = .false.
592 : end if
593 :
594 182485509 : else if (iter == 1) then
595 91242790 : if (ynetpos <= 0.0_r8) then
596 : ! the lower and upper bound ph values (2.0 and 7.0) do not bracket
597 : ! the correct ph, so use the lower bound
598 : tmp_hp = xph(i,k)
599 : converged = .true.
600 : exit
601 : end if
602 : ynetpos_lo = ynetpos
603 :
604 : else ! (iter == 2)
605 91242719 : if (ynetpos >= 0.0_r8) then
606 : ! the lower and upper bound ph values (2.0 and 7.0) do not bracket
607 : ! the correct ph, so use they upper bound
608 : tmp_hp = xph(i,k)
609 : converged = .true.
610 : exit
611 : end if
612 1003669909 : ynetpos_hi = ynetpos
613 : end if
614 :
615 : end do ! iter
616 :
617 91242790 : if( .not. converged ) then
618 0 : write(iulog,*) 'SETSOX: pH failed to converge @ (',i,',',k,'), % change=', &
619 0 : 100._r8*delta
620 : end if
621 : else
622 2082781010 : xph(i,k) = 1.e-7_r8
623 : end if
624 : end do col_loop0
625 : end do ver_loop0 ! end pver loop for STEP 0
626 :
627 : !==============================================================
628 : ! ... Now use the actual PH
629 : !==============================================================
630 139982544 : ver_loop1: do k = 1,pver
631 2314006344 : col_loop1: do i = 1,ncol
632 2174023800 : work1(i) = 1._r8 / tfld(i,k) - 1._r8 / 298._r8
633 2174023800 : tz = tfld(i,k)
634 :
635 2174023800 : xl = cldconc%xlwc(i,k)
636 :
637 2174023800 : patm = press(i,k)/101300._r8 ! press is in pascal
638 2174023800 : xam = press(i,k)/(1.38e-23_r8*tz) ! air density /M3
639 :
640 : !-----------------------------------------------------------------------
641 : ! ... hno3
642 : !-----------------------------------------------------------------------
643 2174023800 : xk = 2.1e5_r8 *EXP( 8700._r8*work1(i) )
644 2174023800 : xe = 15.4_r8
645 2174023800 : hehno3(i,k) = xk*(1._r8 + xe/xph(i,k))
646 :
647 : !-----------------------------------------------------------------
648 : ! ... h2o2
649 : !-----------------------------------------------------------------
650 2174023800 : xk = 7.4e4_r8 *EXP( 6621._r8*work1(i) )
651 2174023800 : xe = 2.2e-12_r8 *EXP(-3730._r8*work1(i) )
652 2174023800 : heh2o2(i,k) = xk*(1._r8 + xe/xph(i,k))
653 :
654 : !-----------------------------------------------------------------
655 : ! ... so2
656 : !-----------------------------------------------------------------
657 2174023800 : xk = 1.23_r8 *EXP( 3120._r8*work1(i) )
658 2174023800 : xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) )
659 2174023800 : x2 = 6.0e-8_r8*EXP( 1120._r8*work1(i) )
660 :
661 2174023800 : wrk = xe/xph(i,k)
662 2174023800 : heso2(i,k) = xk*(1._r8 + wrk*(1._r8 + x2/xph(i,k)))
663 :
664 : !-----------------------------------------------------------------
665 : ! ... nh3
666 : !-----------------------------------------------------------------
667 2174023800 : xk = 58._r8 *EXP( 4085._r8*work1(i) )
668 2174023800 : xe = 1.7e-5_r8*EXP(-4325._r8*work1(i) )
669 2174023800 : henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/xkw)
670 :
671 : !-----------------------------------------------------------------
672 : ! ... o3
673 : !-----------------------------------------------------------------
674 2174023800 : xk = 1.15e-2_r8 *EXP( 2560._r8*work1(i) )
675 2174023800 : heo3(i,k) = xk
676 :
677 : !------------------------------------------------------------------------
678 : ! ... for Ho2(g) -> H2o2(a) formation
679 : ! schwartz JGR, 1984, 11589
680 : !------------------------------------------------------------------------
681 2174023800 : kh4 = (kh2 + kh3*kh1/xph(i,k)) / ((1._r8 + kh1/xph(i,k))**2)
682 2174023800 : ho2s = kh0*xho2(i,k)*patm*(1._r8 + kh1/xph(i,k)) ! ho2s = ho2(a)+o2-
683 2174023800 : r1h2o2 = kh4*ho2s*ho2s ! prod(h2o2) in mole/L(w)/s
684 :
685 2174023800 : if ( cloud_borne ) then
686 : r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s
687 : / const0*1.e+6_r8 & ! correct a bug here ????
688 2174023800 : / xam
689 : else
690 : r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s
691 : * const0 & ! mole/fm3(a)/s * 1.e-3 = mole/cm3(a)/s
692 0 : / xam ! /cm3(a)/s / air-den = mix-ratio/s
693 : endif
694 :
695 2174023800 : if ( .not. modal_aerosols ) then
696 0 : xh2o2(i,k) = xh2o2(i,k) + r2h2o2*dtime ! updated h2o2 by het production
697 : endif
698 :
699 : !-----------------------------------------------
700 : ! ... Partioning
701 : !-----------------------------------------------
702 :
703 : !-----------------------------------------------------------------
704 : ! ... hno3
705 : !-----------------------------------------------------------------
706 2174023800 : px = hehno3(i,k) * Ra * tz * xl
707 2174023800 : hno3g(i,k) = (xhno3(i,k)+xno3(i,k))/(1._r8 + px)
708 :
709 : !------------------------------------------------------------------------
710 : ! ... h2o2
711 : !------------------------------------------------------------------------
712 2174023800 : px = heh2o2(i,k) * Ra * tz * xl
713 2174023800 : h2o2g = xh2o2(i,k)/(1._r8+ px)
714 :
715 : !------------------------------------------------------------------------
716 : ! ... so2
717 : !------------------------------------------------------------------------
718 2174023800 : px = heso2(i,k) * Ra * tz * xl
719 2174023800 : so2g = xso2(i,k)/(1._r8+ px)
720 :
721 : !------------------------------------------------------------------------
722 : ! ... o3
723 : !------------------------------------------------------------------------
724 2174023800 : px = heo3(i,k) * Ra * tz * xl
725 2174023800 : o3g = xo3(i,k)/(1._r8+ px)
726 :
727 : !------------------------------------------------------------------------
728 : ! ... nh3
729 : !------------------------------------------------------------------------
730 2174023800 : px = henh3(i,k) * Ra * tz * xl
731 2174023800 : if (id_nh3>0) then
732 0 : nh3g(i,k) = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px)
733 : else
734 2174023800 : nh3g(i,k) = 0._r8
735 : endif
736 :
737 : !-----------------------------------------------
738 : ! ... Aqueous phase reaction rates
739 : ! SO2 + H2O2 -> SO4
740 : ! SO2 + O3 -> SO4
741 : !-----------------------------------------------
742 :
743 : !------------------------------------------------------------------------
744 : ! ... S(IV) (HSO3) + H2O2
745 : !------------------------------------------------------------------------
746 : rah2o2 = 8.e4_r8 * EXP( -3650._r8*work1(i) ) &
747 2174023800 : / (.1_r8 + xph(i,k))
748 :
749 : !------------------------------------------------------------------------
750 : ! ... S(IV)+ O3
751 : !------------------------------------------------------------------------
752 : rao3 = 4.39e11_r8 * EXP(-4131._r8/tz) &
753 2174023800 : + 2.56e3_r8 * EXP(-996._r8 /tz) /xph(i,k)
754 :
755 : !-----------------------------------------------------------------
756 : ! ... Prediction after aqueous phase
757 : ! so4
758 : ! When Cloud is present
759 : !
760 : ! S(IV) + H2O2 = S(VI)
761 : ! S(IV) + O3 = S(VI)
762 : !
763 : ! reference:
764 : ! (1) Seinfeld
765 : ! (2) Benkovitz
766 : !-----------------------------------------------------------------
767 :
768 : !............................
769 : ! S(IV) + H2O2 = S(VI)
770 : !............................
771 :
772 2312517168 : IF (XL .ge. 1.e-8_r8) THEN !! WHEN CLOUD IS PRESENTED
773 :
774 91242790 : if (cloud_borne) then
775 : patm_x = patm
776 : else
777 0 : patm_x = 1._r8
778 : endif
779 :
780 91242790 : if (modal_aerosols) then
781 :
782 : pso4 = rah2o2 * 7.4e4_r8*EXP(6621._r8*work1(i)) * h2o2g * patm_x &
783 91242790 : * 1.23_r8 *EXP(3120._r8*work1(i)) * so2g * patm_x
784 : else
785 : pso4 = rah2o2 * heh2o2(i,k) * h2o2g * patm_x &
786 0 : * heso2(i,k) * so2g * patm_x ! [M/s]
787 :
788 : endif
789 :
790 : pso4 = pso4 & ! [M/s] = [mole/L(w)/s]
791 : * xl & ! [mole/L(a)/s]
792 : / const0 & ! [/L(a)/s]
793 91242790 : / xhnm(i,k)
794 :
795 :
796 91242790 : ccc = pso4*dtime
797 91242790 : ccc = max(ccc, 1.e-30_r8)
798 :
799 91242790 : xso4_init(i,k)=xso4(i,k)
800 :
801 91242790 : IF (xh2o2(i,k) .gt. xso2(i,k)) THEN
802 63561755 : if (ccc .gt. xso2(i,k)) then
803 8878 : xso4(i,k)=xso4(i,k)+xso2(i,k)
804 8878 : if (cloud_borne) then
805 8878 : xh2o2(i,k)=xh2o2(i,k)-xso2(i,k)
806 8878 : xso2(i,k)=1.e-20_r8
807 : else ! ???? bug ????
808 0 : xso2(i,k)=1.e-20_r8
809 0 : xh2o2(i,k)=xh2o2(i,k)-xso2(i,k)
810 : endif
811 : else
812 63552877 : xso4(i,k) = xso4(i,k) + ccc
813 63552877 : xh2o2(i,k) = xh2o2(i,k) - ccc
814 63552877 : xso2(i,k) = xso2(i,k) - ccc
815 : end if
816 :
817 : ELSE
818 27681035 : if (ccc .gt. xh2o2(i,k)) then
819 2398781 : xso4(i,k)=xso4(i,k)+xh2o2(i,k)
820 2398781 : xso2(i,k)=xso2(i,k)-xh2o2(i,k)
821 2398781 : xh2o2(i,k)=1.e-20_r8
822 : else
823 25282254 : xso4(i,k) = xso4(i,k) + ccc
824 25282254 : xh2o2(i,k) = xh2o2(i,k) - ccc
825 25282254 : xso2(i,k) = xso2(i,k) - ccc
826 : end if
827 : END IF
828 :
829 91242790 : if (modal_aerosols) then
830 91242790 : xdelso4hp(i,k) = xso4(i,k) - xso4_init(i,k)
831 : endif
832 : !...........................
833 : ! S(IV) + O3 = S(VI)
834 : !...........................
835 :
836 91242790 : pso4 = rao3 * heo3(i,k)*o3g*patm_x * heso2(i,k)*so2g*patm_x ! [M/s]
837 :
838 : pso4 = pso4 & ! [M/s] = [mole/L(w)/s]
839 : * xl & ! [mole/L(a)/s]
840 : / const0 & ! [/L(a)/s]
841 91242790 : / xhnm(i,k) ! [mixing ratio/s]
842 :
843 91242790 : ccc = pso4*dtime
844 91242790 : ccc = max(ccc, 1.e-30_r8)
845 :
846 91242790 : xso4_init(i,k)=xso4(i,k)
847 :
848 91242790 : if (ccc .gt. xso2(i,k)) then
849 17218788 : xso4(i,k) = xso4(i,k) + xso2(i,k)
850 17218788 : xso2(i,k) = 1.e-20_r8
851 : else
852 74024002 : xso4(i,k) = xso4(i,k) + ccc
853 74024002 : xso2(i,k) = xso2(i,k) - ccc
854 : end if
855 :
856 : END IF !! WHEN CLOUD IS PRESENTED
857 :
858 : end do col_loop1
859 : end do ver_loop1
860 :
861 : call sox_cldaero_update( &
862 : ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, cldconc%xlwc, &
863 : xdelso4hp, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, &
864 4467528 : aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d=aqso4_h2o2_3d, aqso4_o3_3d=aqso4_o3_3d )
865 :
866 2314006344 : xphlwc(:,:) = 0._r8
867 139982544 : do k = 1, pver
868 2314006344 : do i = 1, ncol
869 2312517168 : if (cldfrc(i,k)>=1.e-5_r8 .and. lwc(i,k)>=1.e-8_r8) then
870 97551055 : xphlwc(i,k) = -1._r8*log10(xph(i,k)) * lwc(i,k)
871 : endif
872 : end do
873 : end do
874 :
875 1489176 : call sox_cldaero_destroy_obj(cldconc)
876 :
877 2978352 : end subroutine SETSOX
878 :
879 : end module MO_SETSOX
|
