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 3072 : subroutine sox_cldaero_init
38 :
39 : integer :: l, m
40 : logical :: history_aerosol ! Output the MAM aerosol tendencies
41 :
42 1536 : id_msa = get_spc_ndx( 'MSA' )
43 1536 : id_h2so4 = get_spc_ndx( 'H2SO4' )
44 1536 : id_so2 = get_spc_ndx( 'SO2' )
45 1536 : id_h2o2 = get_spc_ndx( 'H2O2' )
46 1536 : id_nh3 = get_spc_ndx( 'NH3' )
47 :
48 1536 : 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 1536 : call phys_getopts( history_aerosol_out = history_aerosol )
54 : !
55 : ! add to history
56 : !
57 :
58 1536 : end subroutine sox_cldaero_init
59 :
60 : !----------------------------------------------------------------------------------
61 : !----------------------------------------------------------------------------------
62 58824 : 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 58824 : mode7 = ntot_amode == 7
82 :
83 58824 : conc_obj => cldaero_allocate()
84 :
85 5529456 : do k = 1,pver
86 91405656 : do i = 1,ncol
87 91346832 : if( cldfrc(i,k) >0._r8) then
88 17717220 : conc_obj%xlwc(i,k) = lwc(i,k) *cfact(i,k) ! cloud water L(water)/L(air)
89 17717220 : conc_obj%xlwc(i,k) = conc_obj%xlwc(i,k) / cldfrc(i,k) ! liquid water in the cloudy fraction of cell
90 : else
91 68158980 : conc_obj%xlwc(i,k) = 0._r8
92 : endif
93 : enddo
94 : enddo
95 :
96 93059568 : conc_obj%no3c(:,:) = 0._r8
97 :
98 58824 : 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 58824 : id_so4_1a = lptr_so4_cw_amode(1) - loffset
132 58824 : id_so4_2a = lptr_so4_cw_amode(2) - loffset
133 58824 : id_so4_3a = lptr_so4_cw_amode(3) - loffset
134 0 : conc_obj%so4c(:ncol,:) &
135 58824 : = qcw(:,:,id_so4_1a) &
136 58824 : + qcw(:,:,id_so4_2a) &
137 91523304 : + 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 91405656 : conc_obj%nh4c(:ncol,:) = 0._r8
143 :
144 : ! with 3-mode, assume so4 is nh4hso4, and so half-neutralized
145 58824 : conc_obj%so4_fact = 1._r8
146 :
147 : endif
148 :
149 58824 : end function sox_cldaero_create_obj
150 :
151 : !----------------------------------------------------------------------------------
152 : ! Update the mixing ratios
153 : !----------------------------------------------------------------------------------
154 58824 : subroutine sox_cldaero_update( &
155 176472 : state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, xlwc, &
156 117648 : delso4_hprxn, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, &
157 117648 : 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 117648 : real(r8) :: dqdt_aqso4(ncol,pver,gas_pcnst), &
209 117648 : dqdt_aqh2so4(ncol,pver,gas_pcnst), &
210 117648 : dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver), &
211 : sflx(1:ncol)
212 :
213 58824 : 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 2833634160 : dqdt_aqso4(:,:,:) = 0.0_r8
232 2833634160 : dqdt_aqh2so4(:,:,:) = 0.0_r8
233 91405656 : dqdt_aqhprxn(:,:) = 0.0_r8
234 91405656 : 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 58824 : if ( cam_chempkg_is('geoschem_mam4') ) return
241 :
242 5529456 : lev_loop: do k = 1,pver
243 91405656 : col_loop: do i = 1,ncol
244 91346832 : cloud: if (cldfrc(i,k) >= 1.0e-5_r8) then
245 17319407 : xl = xlwc(i,k) ! / cldfrc(i,k)
246 :
247 17319407 : IF (XL .ge. 1.e-8_r8) THEN !! WHEN CLOUD IS PRESENTED
248 :
249 3382701 : delso4_o3rxn = xso4(i,k) - xso4_init(i,k)
250 :
251 3382701 : 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 16913505 : do n = 1, ntot_amode
268 13530804 : qnum_c(n) = 0.0_r8
269 13530804 : l = numptrcw_amode(n) - loffset
270 16913505 : 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 3382701 : n = modeptr_accum
275 3382701 : if (n <= 0) n = 1
276 3382701 : 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 3382701 : sumf = 0.0_r8
281 16913505 : do n = 1, ntot_amode
282 13530804 : faqgain_so4(n) = 0.0_r8
283 16913505 : if (lptr_so4_cw_amode(n) > 0) then
284 10148103 : faqgain_so4(n) = qnum_c(n)
285 10148103 : sumf = sumf + faqgain_so4(n)
286 : end if
287 : end do
288 :
289 3382701 : if (sumf > 0.0_r8) then
290 16913505 : do n = 1, ntot_amode
291 16913505 : 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 16913505 : do n = 1, ntot_amode
300 13530804 : faqgain_msa(n) = 0.0_r8
301 13530804 : 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 16913505 : sumf = sumf + faqgain_msa(n)
306 : end do
307 :
308 3382701 : 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 3382701 : 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 3382701 : 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 3382701 : dso4dt_gasuptk = xh2so4(i,k) * uptkrate
322 3382701 : 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 3382701 : dmsadt_gasuptk_toso4 = 0.0_r8
330 3382701 : dmsadt_gasuptk_tomsa = dmsadt_gasuptk
331 3382701 : if (ntot_msa_c == 0) then
332 3382701 : dmsadt_gasuptk_tomsa = 0.0_r8
333 3382701 : 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 3382701 : dso4dt_aqrxn = (delso4_o3rxn + delso4_hprxn(i,k)) / dtime
343 3382701 : 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 3382701 : 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 16913505 : do n = 1, ntot_amode
351 13530804 : l = lptr_so4_cw_amode(n) - loffset
352 13530804 : if (l > 0) then
353 10148103 : dqdt_aqso4(i,k,l) = faqgain_so4(n)*dso4dt_aqrxn*cldfrc(i,k)
354 10148103 : dqdt_aqh2so4(i,k,l) = faqgain_so4(n)* &
355 10148103 : (dso4dt_gasuptk + dmsadt_gasuptk_toso4)*cldfrc(i,k)
356 10148103 : dqdt_aq = dqdt_aqso4(i,k,l) + dqdt_aqh2so4(i,k,l)
357 10148103 : dqdt_wr = -fwetrem*dqdt_aq
358 10148103 : dqdt= dqdt_aq + dqdt_wr
359 10148103 : qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
360 : end if
361 :
362 13530804 : l = lptr_msa_cw_amode(n) - loffset
363 13530804 : 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 13530804 : l = lptr_nh4_cw_amode(n) - loffset
371 16913505 : 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 3382701 : qin(i,k,id_h2so4) = qin(i,k,id_h2so4) - dso4dt_gasuptk * dtime * cldfrc(i,k)
393 3382701 : 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 3382701 : fwetrem = 0.0_r8 ! don't include so2 wet removal here
398 :
399 3382701 : dqdt_wr = -fwetrem*xso2(i,k)/dtime*cldfrc(i,k)
400 3382701 : dqdt_aq = -dso4dt_aqrxn*cldfrc(i,k)
401 3382701 : dqdt = dqdt_aq + dqdt_wr
402 3382701 : 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 3382701 : fwetrem = 0.0_r8 ! don't include h2o2 wet removal here
407 :
408 3382701 : dqdt_wr = -fwetrem*xh2o2(i,k)/dtime*cldfrc(i,k)
409 3382701 : dqdt_aq = -dso4dt_hprxn*cldfrc(i,k)
410 3382701 : dqdt = dqdt_aq + dqdt_wr
411 3382701 : qin(i,k,id_h2o2) = qin(i,k,id_h2o2) + dqdt * dtime
412 :
413 : ! NH3
414 3382701 : 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 3382701 : dqdt_aqhprxn(i,k) = dso4dt_hprxn*cldfrc(i,k)
422 3382701 : 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 5529456 : do k = 1,pver
433 :
434 27353160 : do n = 1, ntot_amode
435 :
436 21882528 : l = lptr_so4_cw_amode(n) - loffset
437 21882528 : if (l > 0) then
438 274040496 : qcw(:,k,l) = MAX(qcw(:,k,l), small_value )
439 : end if
440 21882528 : l = lptr_msa_cw_amode(n) - loffset
441 21882528 : if (l > 0) then
442 0 : qcw(:,k,l) = MAX(qcw(:,k,l), small_value )
443 : end if
444 21882528 : l = lptr_nh4_cw_amode(n) - loffset
445 27353160 : if (l > 0) then
446 0 : qcw(:,k,l) = MAX(qcw(:,k,l), small_value )
447 : end if
448 :
449 : end do
450 :
451 91346832 : qin(:,k,id_so2) = MAX( qin(:,k,id_so2), small_value )
452 :
453 5529456 : 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 294120 : do n = 1, ntot_amode
462 235296 : m = lptr_so4_cw_amode(n)
463 235296 : l = m - loffset
464 294120 : if (l > 0) then
465 2946672 : aqso4(:,n)=0._r8
466 16588368 : do k=1,pver
467 274216968 : do i=1,ncol
468 772885800 : aqso4(i,n)=aqso4(i,n)+dqdt_aqso4(i,k,l)*adv_mass(l)/mbar(i,k) &
469 1046926296 : *pdel(i,k)/gravit ! kg/m2/s
470 : enddo
471 : enddo
472 :
473 2946672 : aqh2so4(:,n)=0._r8
474 16588368 : do k=1,pver
475 274216968 : do i=1,ncol
476 772885800 : aqh2so4(i,n)=aqh2so4(i,n)+dqdt_aqh2so4(i,k,l)*adv_mass(l)/mbar(i,k) &
477 1046926296 : *pdel(i,k)/gravit ! kg/m2/s
478 : enddo
479 : enddo
480 : endif
481 : end do
482 :
483 982224 : aqso4_h2o2(:) = 0._r8
484 5529456 : do k=1,pver
485 91405656 : do i=1,ncol
486 171752400 : aqso4_h2o2(i)=aqso4_h2o2(i)+dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) &
487 263099232 : *pdel(i,k)/gravit ! kg SO4 /m2/s
488 : enddo
489 : enddo
490 :
491 58824 : 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 982224 : aqso4_o3(:)=0._r8
502 5529456 : do k=1,pver
503 91405656 : do i=1,ncol
504 171752400 : aqso4_o3(i)=aqso4_o3(i)+dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) &
505 263099232 : *pdel(i,k)/gravit ! kg SO4 /m2/s
506 : enddo
507 : enddo
508 :
509 58824 : 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 58824 : end subroutine sox_cldaero_update
520 :
521 : !----------------------------------------------------------------------------------
522 : !----------------------------------------------------------------------------------
523 58824 : subroutine sox_cldaero_destroy_obj( conc_obj )
524 : type(cldaero_conc_t), pointer :: conc_obj
525 :
526 58824 : call cldaero_deallocate( conc_obj )
527 :
528 58824 : end subroutine sox_cldaero_destroy_obj
529 :
530 : end module sox_cldaero_mod
|
