Line data Source code
1 : !----------------------------------------------------------------------------------
2 : ! Modal aerosol implementation
3 : !----------------------------------------------------------------------------------
4 : module sox_cldaero_mod
5 :
6 : use shr_kind_mod, only : r8 => shr_kind_r8
7 : use cam_abortutils, only : endrun
8 : use ppgrid, only : pcols, pver
9 : use mo_chem_utls, only : get_spc_ndx
10 : use cldaero_mod, only : cldaero_conc_t, cldaero_allocate, cldaero_deallocate
11 : use modal_aero_data, only : ntot_amode, modeptr_accum, lptr_so4_cw_amode, lptr_msa_cw_amode
12 : use modal_aero_data, only : numptrcw_amode, lptr_nh4_cw_amode
13 : use modal_aero_data, only : cnst_name_cw, specmw_so4_amode
14 : use chem_mods, only : adv_mass
15 : use physconst, only : gravit
16 : use phys_control, only : phys_getopts, cam_chempkg_is
17 : use cldaero_mod, only : cldaero_uptakerate
18 : use chem_mods, only : gas_pcnst
19 :
20 : implicit none
21 : private
22 :
23 : public :: sox_cldaero_init
24 : public :: sox_cldaero_create_obj
25 : public :: sox_cldaero_update
26 : public :: sox_cldaero_destroy_obj
27 :
28 : integer :: id_msa, id_h2so4, id_so2, id_h2o2, id_nh3
29 :
30 : real(r8), parameter :: small_value = 1.e-20_r8
31 :
32 : contains
33 :
34 : !----------------------------------------------------------------------------------
35 : !----------------------------------------------------------------------------------
36 :
37 4608 : subroutine sox_cldaero_init
38 :
39 : integer :: l, m
40 : logical :: history_aerosol ! Output the MAM aerosol tendencies
41 :
42 2304 : id_msa = get_spc_ndx( 'MSA' )
43 2304 : id_h2so4 = get_spc_ndx( 'H2SO4' )
44 2304 : id_so2 = get_spc_ndx( 'SO2' )
45 2304 : id_h2o2 = get_spc_ndx( 'H2O2' )
46 2304 : id_nh3 = get_spc_ndx( 'NH3' )
47 :
48 2304 : if (id_h2so4<1 .or. id_so2<1 .or. id_h2o2<1) then
49 : call endrun('sox_cldaero_init:MAM mech does not include necessary species' &
50 0 : //' -- should not invoke sox_cldaero_mod ')
51 : endif
52 :
53 2304 : call phys_getopts( history_aerosol_out = history_aerosol )
54 : !
55 : ! add to history
56 : !
57 :
58 2304 : end subroutine sox_cldaero_init
59 :
60 : !----------------------------------------------------------------------------------
61 : !----------------------------------------------------------------------------------
62 89784 : function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( conc_obj )
63 :
64 : real(r8), intent(in) :: cldfrc(:,:)
65 : real(r8), intent(in) :: qcw(:,:,:)
66 : real(r8), intent(in) :: lwc(:,:)
67 : real(r8), intent(in) :: cfact(:,:)
68 : integer, intent(in) :: ncol
69 : integer, intent(in) :: loffset
70 :
71 : type(cldaero_conc_t), pointer :: conc_obj
72 :
73 :
74 : integer :: id_so4_1a, id_so4_2a, id_so4_3a, id_so4_4a, id_so4_5a, id_so4_6a
75 : integer :: id_nh4_1a, id_nh4_2a, id_nh4_3a, id_nh4_4a, id_nh4_5a, id_nh4_6a
76 : integer :: l,n
77 : integer :: i,k
78 :
79 : logical :: mode7
80 :
81 89784 : mode7 = ntot_amode == 7
82 :
83 89784 : conc_obj => cldaero_allocate()
84 :
85 8439696 : do k = 1,pver
86 139513896 : do i = 1,ncol
87 139424112 : if( cldfrc(i,k) >0._r8) then
88 23445741 : conc_obj%xlwc(i,k) = lwc(i,k) *cfact(i,k) ! cloud water L(water)/L(air)
89 23445741 : conc_obj%xlwc(i,k) = conc_obj%xlwc(i,k) / cldfrc(i,k) ! liquid water in the cloudy fraction of cell
90 : else
91 107628459 : conc_obj%xlwc(i,k) = 0._r8
92 : endif
93 : enddo
94 : enddo
95 :
96 142038288 : conc_obj%no3c(:,:) = 0._r8
97 :
98 89784 : if (mode7) then
99 : #if ( defined MODAL_AERO_7MODE )
100 : !put ifdef here so ifort will compile
101 : id_so4_1a = lptr_so4_cw_amode(1) - loffset
102 : id_so4_2a = lptr_so4_cw_amode(2) - loffset
103 : id_so4_3a = lptr_so4_cw_amode(4) - loffset
104 : id_so4_4a = lptr_so4_cw_amode(5) - loffset
105 : id_so4_5a = lptr_so4_cw_amode(6) - loffset
106 : id_so4_6a = lptr_so4_cw_amode(7) - loffset
107 :
108 : id_nh4_1a = lptr_nh4_cw_amode(1) - loffset
109 : id_nh4_2a = lptr_nh4_cw_amode(2) - loffset
110 : id_nh4_3a = lptr_nh4_cw_amode(4) - loffset
111 : id_nh4_4a = lptr_nh4_cw_amode(5) - loffset
112 : id_nh4_5a = lptr_nh4_cw_amode(6) - loffset
113 : id_nh4_6a = lptr_nh4_cw_amode(7) - loffset
114 : #endif
115 0 : conc_obj%so4c(:ncol,:) &
116 0 : = qcw(:ncol,:,id_so4_1a) &
117 0 : + qcw(:ncol,:,id_so4_2a) &
118 0 : + qcw(:ncol,:,id_so4_3a) &
119 0 : + qcw(:ncol,:,id_so4_4a) &
120 0 : + qcw(:ncol,:,id_so4_5a) &
121 0 : + qcw(:ncol,:,id_so4_6a)
122 :
123 0 : conc_obj%nh4c(:ncol,:) &
124 0 : = qcw(:ncol,:,id_nh4_1a) &
125 0 : + qcw(:ncol,:,id_nh4_2a) &
126 0 : + qcw(:ncol,:,id_nh4_3a) &
127 0 : + qcw(:ncol,:,id_nh4_4a) &
128 0 : + qcw(:ncol,:,id_nh4_5a) &
129 0 : + qcw(:ncol,:,id_nh4_6a)
130 : else
131 89784 : id_so4_1a = lptr_so4_cw_amode(1) - loffset
132 89784 : id_so4_2a = lptr_so4_cw_amode(2) - loffset
133 89784 : id_so4_3a = lptr_so4_cw_amode(3) - loffset
134 0 : conc_obj%so4c(:ncol,:) &
135 89784 : = qcw(:,:,id_so4_1a) &
136 89784 : + qcw(:,:,id_so4_2a) &
137 139693464 : + qcw(:,:,id_so4_3a)
138 :
139 : ! for 3-mode, so4 is assumed to be nh4hso4
140 : ! the partial neutralization of so4 is handled by using a
141 : ! -1 charge (instead of -2) in the electro-neutrality equation
142 139513896 : conc_obj%nh4c(:ncol,:) = 0._r8
143 :
144 : ! with 3-mode, assume so4 is nh4hso4, and so half-neutralized
145 89784 : conc_obj%so4_fact = 1._r8
146 :
147 : endif
148 :
149 89784 : end function sox_cldaero_create_obj
150 :
151 : !----------------------------------------------------------------------------------
152 : ! Update the mixing ratios
153 : !----------------------------------------------------------------------------------
154 89784 : subroutine sox_cldaero_update( &
155 269352 : state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, xlwc, &
156 179568 : delso4_hprxn, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, &
157 179568 : aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d, aqso4_o3_3d)
158 :
159 : use physics_types, only: physics_state
160 :
161 : ! args
162 :
163 : type(physics_state), intent(in) :: state ! Physics state variables
164 :
165 : integer, intent(in) :: ncol
166 : integer, intent(in) :: lchnk ! chunk id
167 : integer, intent(in) :: loffset
168 :
169 : real(r8), intent(in) :: dtime ! time step (sec)
170 :
171 : real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu )
172 : real(r8), intent(in) :: pdel(:,:)
173 : real(r8), intent(in) :: press(:,:)
174 : real(r8), intent(in) :: tfld(:,:)
175 :
176 : real(r8), intent(in) :: cldnum(:,:)
177 : real(r8), intent(in) :: cldfrc(:,:)
178 : real(r8), intent(in) :: cfact(:,:)
179 : real(r8), intent(in) :: xlwc(:,:)
180 :
181 : real(r8), intent(in) :: delso4_hprxn(:,:)
182 : real(r8), intent(in) :: xh2so4(:,:)
183 : real(r8), intent(in) :: xso4(:,:)
184 : real(r8), intent(in) :: xso4_init(:,:)
185 : real(r8), intent(in) :: nh3g(:,:)
186 : real(r8), intent(in) :: hno3g(:,:)
187 : real(r8), intent(in) :: xnh3(:,:)
188 : real(r8), intent(in) :: xhno3(:,:)
189 : real(r8), intent(in) :: xnh4c(:,:)
190 : real(r8), intent(in) :: xmsa(:,:)
191 : real(r8), intent(in) :: xso2(:,:)
192 : real(r8), intent(in) :: xh2o2(:,:)
193 : real(r8), intent(in) :: xno3c(:,:)
194 :
195 : real(r8), intent(inout) :: qcw(:,:,:) ! cloud-borne aerosol (vmr)
196 : real(r8), intent(inout) :: qin(:,:,:) ! xported species ( vmr )
197 :
198 : real(r8), intent(out) :: aqso4(:,:) ! aqueous phase chemistry
199 : real(r8), intent(out) :: aqh2so4(:,:) ! aqueous phase chemistry
200 : real(r8), intent(out) :: aqso4_h2o2(:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2)
201 : real(r8), intent(out) :: aqso4_o3(:) ! SO4 aqueous phase chemistry due to O3 (kg/m2)
202 : real(r8), intent(out), optional :: aqso4_h2o2_3d(:,:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2)
203 : real(r8), intent(out), optional :: aqso4_o3_3d(:,:) ! SO4 aqueous phase chemistry due to O3 (kg/m2)
204 :
205 :
206 : ! local vars ...
207 :
208 179568 : real(r8) :: dqdt_aqso4(ncol,pver,gas_pcnst), &
209 179568 : dqdt_aqh2so4(ncol,pver,gas_pcnst), &
210 179568 : dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver), &
211 : sflx(1:ncol)
212 :
213 89784 : real(r8) :: faqgain_msa(ntot_amode), faqgain_so4(ntot_amode), qnum_c(ntot_amode)
214 :
215 : real(r8) :: delso4_o3rxn, &
216 : dso4dt_aqrxn, dso4dt_hprxn, &
217 : dso4dt_gasuptk, dmsadt_gasuptk, &
218 : dmsadt_gasuptk_tomsa, dmsadt_gasuptk_toso4, &
219 : dqdt_aq, dqdt_wr, dqdt
220 :
221 : real(r8) :: fwetrem, sumf, uptkrate
222 : real(r8) :: delnh3, delnh4
223 :
224 : integer :: l, n, m
225 : integer :: ntot_msa_c
226 :
227 : integer :: i,k
228 : real(r8) :: xl
229 :
230 : ! make sure dqdt is zero initially, for budgets
231 4325020560 : dqdt_aqso4(:,:,:) = 0.0_r8
232 4325020560 : dqdt_aqh2so4(:,:,:) = 0.0_r8
233 139513896 : dqdt_aqhprxn(:,:) = 0.0_r8
234 139513896 : dqdt_aqo3rxn(:,:) = 0.0_r8
235 :
236 : ! Avoid double counting in-cloud sulfur oxidation when running with
237 : ! GEOS-Chem. If running with GEOS-Chem then sulfur oxidation
238 : ! is performed internally to GEOS-Chem. Here, we just return to the
239 : ! parent routine and thus we do not apply tendencies calculated by MAM.
240 89784 : if ( cam_chempkg_is('geoschem_mam4') ) return
241 :
242 8439696 : lev_loop: do k = 1,pver
243 139513896 : col_loop: do i = 1,ncol
244 139424112 : cloud: if (cldfrc(i,k) >= 1.0e-5_r8) then
245 22953536 : xl = xlwc(i,k) ! / cldfrc(i,k)
246 :
247 22953536 : IF (XL .ge. 1.e-8_r8) THEN !! WHEN CLOUD IS PRESENTED
248 :
249 5010776 : delso4_o3rxn = xso4(i,k) - xso4_init(i,k)
250 :
251 5010776 : if (id_nh3>0) then
252 0 : delnh3 = nh3g(i,k) - xnh3(i,k)
253 0 : delnh4 = - delnh3
254 : endif
255 :
256 : !-------------------------------------------------------------------------
257 : ! compute factors for partitioning aerosol mass gains among modes
258 : ! the factors are proportional to the activated particle MR for each
259 : ! mode, which is the MR of cloud drops "associated with" the mode
260 : ! thus we are assuming the cloud drop size is independent of the
261 : ! associated aerosol mode properties (i.e., drops associated with
262 : ! Aitken and coarse sea-salt particles are same size)
263 : !
264 : ! qnum_c(n) = activated particle number MR for mode n (these are just
265 : ! used for partitioning among modes, so don't need to divide by cldfrc)
266 :
267 25053880 : do n = 1, ntot_amode
268 20043104 : qnum_c(n) = 0.0_r8
269 20043104 : l = numptrcw_amode(n) - loffset
270 25053880 : if (l > 0) qnum_c(n) = max( 0.0_r8, qcw(i,k,l) )
271 : end do
272 :
273 : ! force qnum_c(n) to be positive for n=modeptr_accum or n=1
274 5010776 : n = modeptr_accum
275 5010776 : if (n <= 0) n = 1
276 5010776 : qnum_c(n) = max( 1.0e-10_r8, qnum_c(n) )
277 :
278 : ! faqgain_so4(n) = fraction of total so4_c gain going to mode n
279 : ! these are proportional to the activated particle MR for each mode
280 5010776 : sumf = 0.0_r8
281 25053880 : do n = 1, ntot_amode
282 20043104 : faqgain_so4(n) = 0.0_r8
283 25053880 : if (lptr_so4_cw_amode(n) > 0) then
284 15032328 : faqgain_so4(n) = qnum_c(n)
285 15032328 : sumf = sumf + faqgain_so4(n)
286 : end if
287 : end do
288 :
289 5010776 : if (sumf > 0.0_r8) then
290 25053880 : do n = 1, ntot_amode
291 25053880 : faqgain_so4(n) = faqgain_so4(n) / sumf
292 : end do
293 : end if
294 : ! at this point (sumf <= 0.0) only when all the faqgain_so4 are zero
295 :
296 : ! faqgain_msa(n) = fraction of total msa_c gain going to mode n
297 : ntot_msa_c = 0
298 : sumf = 0.0_r8
299 25053880 : do n = 1, ntot_amode
300 20043104 : faqgain_msa(n) = 0.0_r8
301 20043104 : if (lptr_msa_cw_amode(n) > 0) then
302 0 : faqgain_msa(n) = qnum_c(n)
303 0 : ntot_msa_c = ntot_msa_c + 1
304 : end if
305 25053880 : sumf = sumf + faqgain_msa(n)
306 : end do
307 :
308 5010776 : if (sumf > 0.0_r8) then
309 0 : do n = 1, ntot_amode
310 0 : faqgain_msa(n) = faqgain_msa(n) / sumf
311 : end do
312 : end if
313 : ! at this point (sumf <= 0.0) only when all the faqgain_msa are zero
314 :
315 5010776 : uptkrate = cldaero_uptakerate( xl, cldnum(i,k), cfact(i,k), cldfrc(i,k), tfld(i,k), press(i,k) )
316 : ! average uptake rate over dtime
317 5010776 : uptkrate = (1.0_r8 - exp(-min(100._r8,dtime*uptkrate))) / dtime
318 :
319 : ! dso4dt_gasuptk = so4_c tendency from h2so4 gas uptake (mol/mol/s)
320 : ! dmsadt_gasuptk = msa_c tendency from msa gas uptake (mol/mol/s)
321 5010776 : dso4dt_gasuptk = xh2so4(i,k) * uptkrate
322 5010776 : if (id_msa > 0) then
323 0 : dmsadt_gasuptk = xmsa(i,k) * uptkrate
324 : else
325 : dmsadt_gasuptk = 0.0_r8
326 : end if
327 :
328 : ! if no modes have msa aerosol, then "rename" scavenged msa gas to so4
329 5010776 : dmsadt_gasuptk_toso4 = 0.0_r8
330 5010776 : dmsadt_gasuptk_tomsa = dmsadt_gasuptk
331 5010776 : if (ntot_msa_c == 0) then
332 5010776 : dmsadt_gasuptk_tomsa = 0.0_r8
333 5010776 : dmsadt_gasuptk_toso4 = dmsadt_gasuptk
334 : end if
335 :
336 : !-----------------------------------------------------------------------
337 : ! now compute TMR tendencies
338 : ! this includes the above aqueous so2 chemistry AND
339 : ! the uptake of highly soluble aerosol precursor gases (h2so4, msa, ...)
340 : ! AND the wetremoval of dissolved, unreacted so2 and h2o2
341 :
342 5010776 : dso4dt_aqrxn = (delso4_o3rxn + delso4_hprxn(i,k)) / dtime
343 5010776 : dso4dt_hprxn = delso4_hprxn(i,k) / dtime
344 :
345 : ! fwetrem = fraction of in-cloud-water material that is wet removed
346 : ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*clwlrat(i,k)))) )
347 5010776 : fwetrem = 0.0_r8 ! don't have so4 & msa wet removal here
348 :
349 : ! compute TMR tendencies for so4 and msa aerosol-in-cloud-water
350 25053880 : do n = 1, ntot_amode
351 20043104 : l = lptr_so4_cw_amode(n) - loffset
352 20043104 : if (l > 0) then
353 15032328 : dqdt_aqso4(i,k,l) = faqgain_so4(n)*dso4dt_aqrxn*cldfrc(i,k)
354 15032328 : dqdt_aqh2so4(i,k,l) = faqgain_so4(n)* &
355 15032328 : (dso4dt_gasuptk + dmsadt_gasuptk_toso4)*cldfrc(i,k)
356 15032328 : dqdt_aq = dqdt_aqso4(i,k,l) + dqdt_aqh2so4(i,k,l)
357 15032328 : dqdt_wr = -fwetrem*dqdt_aq
358 15032328 : dqdt= dqdt_aq + dqdt_wr
359 15032328 : qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
360 : end if
361 :
362 20043104 : l = lptr_msa_cw_amode(n) - loffset
363 20043104 : if (l > 0) then
364 0 : dqdt_aq = faqgain_msa(n)*dmsadt_gasuptk_tomsa*cldfrc(i,k)
365 0 : dqdt_wr = -fwetrem*dqdt_aq
366 0 : dqdt = dqdt_aq + dqdt_wr
367 0 : qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
368 : end if
369 :
370 20043104 : l = lptr_nh4_cw_amode(n) - loffset
371 25053880 : if (l > 0) then
372 0 : if (delnh4 > 0.0_r8) then
373 0 : dqdt_aq = faqgain_so4(n)*delnh4/dtime*cldfrc(i,k)
374 0 : dqdt = dqdt_aq
375 0 : qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
376 : else
377 0 : dqdt = (qcw(i,k,l)/max(xnh4c(i,k),1.0e-35_r8)) &
378 0 : *delnh4/dtime*cldfrc(i,k)
379 0 : qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
380 : endif
381 : end if
382 : end do
383 :
384 : ! For gas species, tendency includes
385 : ! reactive uptake to cloud water that essentially transforms the gas to
386 : ! a different species. Wet removal associated with this is applied
387 : ! to the "new" species (e.g., so4_c) rather than to the gas.
388 : ! wet removal of the unreacted gas that is dissolved in cloud water.
389 : ! Need to multiply both these parts by cldfrc
390 :
391 : ! h2so4 (g) & msa (g)
392 5010776 : qin(i,k,id_h2so4) = qin(i,k,id_h2so4) - dso4dt_gasuptk * dtime * cldfrc(i,k)
393 5010776 : if (id_msa > 0) qin(i,k,id_msa) = qin(i,k,id_msa) - dmsadt_gasuptk * dtime * cldfrc(i,k)
394 :
395 : ! so2 -- the first order loss rate for so2 is frso2_c*clwlrat(i,k)
396 : ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frso2_c*clwlrat(i,k)))) )
397 5010776 : fwetrem = 0.0_r8 ! don't include so2 wet removal here
398 :
399 5010776 : dqdt_wr = -fwetrem*xso2(i,k)/dtime*cldfrc(i,k)
400 5010776 : dqdt_aq = -dso4dt_aqrxn*cldfrc(i,k)
401 5010776 : dqdt = dqdt_aq + dqdt_wr
402 5010776 : qin(i,k,id_so2) = qin(i,k,id_so2) + dqdt * dtime
403 :
404 : ! h2o2 -- the first order loss rate for h2o2 is frh2o2_c*clwlrat(i,k)
405 : ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frh2o2_c*clwlrat(i,k)))) )
406 5010776 : fwetrem = 0.0_r8 ! don't include h2o2 wet removal here
407 :
408 5010776 : dqdt_wr = -fwetrem*xh2o2(i,k)/dtime*cldfrc(i,k)
409 5010776 : dqdt_aq = -dso4dt_hprxn*cldfrc(i,k)
410 5010776 : dqdt = dqdt_aq + dqdt_wr
411 5010776 : qin(i,k,id_h2o2) = qin(i,k,id_h2o2) + dqdt * dtime
412 :
413 : ! NH3
414 5010776 : if (id_nh3>0) then
415 0 : dqdt_aq = delnh3/dtime*cldfrc(i,k)
416 0 : dqdt = dqdt_aq
417 0 : qin(i,k,id_nh3) = qin(i,k,id_nh3) + dqdt * dtime
418 : endif
419 :
420 : ! for SO4 from H2O2/O3 budgets
421 5010776 : dqdt_aqhprxn(i,k) = dso4dt_hprxn*cldfrc(i,k)
422 5010776 : dqdt_aqo3rxn(i,k) = (dso4dt_aqrxn - dso4dt_hprxn)*cldfrc(i,k)
423 :
424 : ENDIF !! WHEN CLOUD IS PRESENTED
425 : endif cloud
426 : enddo col_loop
427 : enddo lev_loop
428 :
429 : !==============================================================
430 : ! ... Update the mixing ratios
431 : !==============================================================
432 8439696 : do k = 1,pver
433 :
434 41749560 : do n = 1, ntot_amode
435 :
436 33399648 : l = lptr_so4_cw_amode(n) - loffset
437 33399648 : if (l > 0) then
438 418272336 : qcw(:,k,l) = MAX(qcw(:,k,l), small_value )
439 : end if
440 33399648 : l = lptr_msa_cw_amode(n) - loffset
441 33399648 : if (l > 0) then
442 0 : qcw(:,k,l) = MAX(qcw(:,k,l), small_value )
443 : end if
444 33399648 : l = lptr_nh4_cw_amode(n) - loffset
445 41749560 : if (l > 0) then
446 0 : qcw(:,k,l) = MAX(qcw(:,k,l), small_value )
447 : end if
448 :
449 : end do
450 :
451 139424112 : qin(:,k,id_so2) = MAX( qin(:,k,id_so2), small_value )
452 :
453 8439696 : if ( id_nh3 > 0 ) then
454 0 : qin(:,k,id_nh3) = MAX( qin(:,k,id_nh3), small_value )
455 : endif
456 :
457 : end do
458 :
459 : ! diagnostics
460 :
461 448920 : do n = 1, ntot_amode
462 359136 : m = lptr_so4_cw_amode(n)
463 359136 : l = m - loffset
464 448920 : if (l > 0) then
465 4497552 : aqso4(:,n)=0._r8
466 25319088 : do k=1,pver
467 418541688 : do i=1,ncol
468 1179667800 : aqso4(i,n)=aqso4(i,n)+dqdt_aqso4(i,k,l)*adv_mass(l)/mbar(i,k) &
469 1597940136 : *pdel(i,k)/gravit ! kg/m2/s
470 : enddo
471 : enddo
472 :
473 4497552 : aqh2so4(:,n)=0._r8
474 25319088 : do k=1,pver
475 418541688 : do i=1,ncol
476 1179667800 : aqh2so4(i,n)=aqh2so4(i,n)+dqdt_aqh2so4(i,k,l)*adv_mass(l)/mbar(i,k) &
477 1597940136 : *pdel(i,k)/gravit ! kg/m2/s
478 : enddo
479 : enddo
480 : endif
481 : end do
482 :
483 1499184 : aqso4_h2o2(:) = 0._r8
484 8439696 : do k=1,pver
485 139513896 : do i=1,ncol
486 262148400 : aqso4_h2o2(i)=aqso4_h2o2(i)+dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) &
487 401572512 : *pdel(i,k)/gravit ! kg SO4 /m2/s
488 : enddo
489 : enddo
490 :
491 89784 : if (present(aqso4_h2o2_3d)) then
492 0 : aqso4_h2o2_3d(:,:) = 0._r8
493 0 : do k=1,pver
494 0 : do i=1,ncol
495 0 : aqso4_h2o2_3d(i,k)=dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) &
496 0 : *pdel(i,k)/gravit ! kg SO4 /m2/s
497 : enddo
498 : enddo
499 : end if
500 :
501 1499184 : aqso4_o3(:)=0._r8
502 8439696 : do k=1,pver
503 139513896 : do i=1,ncol
504 262148400 : aqso4_o3(i)=aqso4_o3(i)+dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) &
505 401572512 : *pdel(i,k)/gravit ! kg SO4 /m2/s
506 : enddo
507 : enddo
508 :
509 89784 : if (present(aqso4_o3_3d)) then
510 0 : aqso4_o3_3d(:,:)=0._r8
511 0 : do k=1,pver
512 0 : do i=1,ncol
513 0 : aqso4_o3_3d(i,k)=dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) &
514 0 : *pdel(i,k)/gravit ! kg SO4 /m2/s
515 : enddo
516 : enddo
517 : end if
518 :
519 89784 : end subroutine sox_cldaero_update
520 :
521 : !----------------------------------------------------------------------------------
522 : !----------------------------------------------------------------------------------
523 89784 : subroutine sox_cldaero_destroy_obj( conc_obj )
524 : type(cldaero_conc_t), pointer :: conc_obj
525 :
526 89784 : call cldaero_deallocate( conc_obj )
527 :
528 89784 : end subroutine sox_cldaero_destroy_obj
529 :
530 : end module sox_cldaero_mod
|
