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