Line data Source code
1 : !----------------------------------------------------------------------------------
2 : ! CARMA 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 cam_logfile, only : iulog
12 : use chem_mods, only : adv_mass
13 : use physconst, only : gravit
14 : use phys_control, only : phys_getopts
15 : use cldaero_mod, only : cldaero_uptakerate
16 : use chem_mods, only : gas_pcnst
17 : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_info_by_bin, rad_cnst_get_bin_props_by_idx
18 :
19 : implicit none
20 : private
21 :
22 : public :: sox_cldaero_init
23 : public :: sox_cldaero_create_obj
24 : public :: sox_cldaero_update
25 : public :: sox_cldaero_destroy_obj
26 :
27 : integer :: id_msa, id_h2so4, id_so2, id_h2o2, id_nh3
28 :
29 : real(r8), parameter :: small_value = 1.e-20_r8
30 :
31 : ! description of bin aerosols
32 : integer, public, protected :: nspec_max = 0
33 : integer, public, protected :: nbins = 0
34 : integer, public, protected, allocatable :: nspec(:)
35 :
36 : ! local indexing for bins
37 : integer, allocatable :: bin_idx(:,:) ! table for local indexing of modal aero number and mmr
38 : integer :: ncnst_tot ! total number of mode number conc + mode species
39 :
40 : contains
41 :
42 : !----------------------------------------------------------------------------------
43 : !----------------------------------------------------------------------------------
44 :
45 1536 : subroutine sox_cldaero_init
46 :
47 : integer :: l, m, ii
48 : logical :: history_aerosol ! Output the MAM aerosol tendencies
49 :
50 1536 : id_msa = get_spc_ndx( 'MSA' )
51 1536 : id_h2so4 = get_spc_ndx( 'H2SO4' )
52 1536 : id_so2 = get_spc_ndx( 'SO2' )
53 1536 : id_h2o2 = get_spc_ndx( 'H2O2' )
54 1536 : id_nh3 = get_spc_ndx( 'NH3' )
55 :
56 1536 : if (id_h2so4<1 .or. id_so2<1 .or. id_h2o2<1) then
57 : call endrun('sox_cldaero_init:MAM mech does not include necessary species' &
58 0 : //' -- should not invoke sox_cldaero_mod ')
59 : endif
60 :
61 1536 : call phys_getopts( history_aerosol_out = history_aerosol )
62 : !
63 : ! add to history
64 : !
65 :
66 : ! get info about the modal aerosols
67 : ! get nbins
68 :
69 1536 : call rad_cnst_get_info( 0, nbins=nbins)
70 :
71 4608 : allocate( nspec(nbins) )
72 :
73 62976 : do m = 1, nbins
74 62976 : call rad_cnst_get_info_by_bin(0, m, nspec=nspec(m))
75 : end do
76 : ! add plus one to include number, total mmr and nspec
77 62976 : nspec_max = maxval(nspec)
78 :
79 1536 : ncnst_tot = nspec(1)
80 61440 : do m = 2, nbins
81 61440 : ncnst_tot = ncnst_tot + nspec(m)
82 : end do
83 :
84 6144 : allocate( bin_idx(nbins,nspec_max) )
85 :
86 :
87 : ! Local indexing compresses the mode and number/mass indicies into one index.
88 : ! This indexing is used by the pointer arrays used to reference state and pbuf
89 : ! fields.
90 : ! for CARMA we add number = 0, total mass = 1, and mass from each constituence into mm.
91 1536 : ii = 0
92 62976 : do m = 1, nbins
93 278016 : do l = 1, nspec(m) ! loop through species
94 215040 : ii = ii + 1
95 276480 : bin_idx(m,l) = ii
96 : end do
97 : end do
98 :
99 :
100 1536 : end subroutine sox_cldaero_init
101 :
102 : !----------------------------------------------------------------------------------
103 : !----------------------------------------------------------------------------------
104 72960 : function sox_cldaero_create_obj(cldfrc, qcw, lwc, cfact, ncol, loffset) result( conc_obj )
105 :
106 : real(r8), intent(in) :: cldfrc(:,:)
107 : real(r8), intent(in) :: qcw(:,:,:)
108 : real(r8), intent(in) :: lwc(:,:)
109 : real(r8), intent(in) :: cfact(:,:)
110 : integer, intent(in) :: ncol
111 : integer, intent(in) :: loffset
112 :
113 : real(r8) :: so4mmr(pcols,pver)
114 :
115 : type(cldaero_conc_t), pointer :: conc_obj
116 :
117 : character(len=32) :: spectype
118 :
119 : integer :: l,m
120 : integer :: i,k,mm
121 :
122 : ! local indexing for bins
123 : !integer, allocatable :: bin_idx(:,:) ! table for local indexing of modal aero number and mmr
124 :
125 :
126 72960 : conc_obj => cldaero_allocate()
127 :
128 5180160 : do k = 1,pver
129 78723840 : do i = 1,ncol
130 78650880 : if( cldfrc(i,k) >0._r8) then
131 10663221 : conc_obj%xlwc(i,k) = lwc(i,k) *cfact(i,k) ! cloud water L(water)/L(air)
132 10663221 : conc_obj%xlwc(i,k) = conc_obj%xlwc(i,k) / cldfrc(i,k) ! liquid water in the cloudy fraction of cell
133 : else
134 62880459 : conc_obj%xlwc(i,k) = 0._r8
135 : endif
136 : enddo
137 : enddo
138 :
139 86895360 : conc_obj%no3c(:,:) = 0._r8
140 86895360 : conc_obj%nh4c(:,:) = 0._r8
141 86895360 : conc_obj%so4c(:,:) = 0._r8
142 :
143 72960 : so4mmr(:,:) = 0._r8
144 5180160 : do k = 1,pver
145 78723840 : do i = 1,ncol
146 3020398080 : do m = 1, nbins
147 13311406080 : do l = 1, nspec(m)
148 10296115200 : mm = bin_idx(m, l)
149 10296115200 : call rad_cnst_get_bin_props_by_idx(0, m, l,spectype=spectype)
150 13237862400 : if (trim(spectype) == 'sulfate') then
151 2941747200 : so4mmr(i,k) = so4mmr(i,k) + qcw(i,k,mm)
152 : end if
153 : end do
154 : end do
155 : end do
156 : end do
157 86895360 : conc_obj%so4c = so4mmr
158 :
159 72960 : end function sox_cldaero_create_obj
160 :
161 :
162 : !----------------------------------------------------------------------------------
163 : ! Update the mixing ratios
164 : !----------------------------------------------------------------------------------
165 72960 : subroutine sox_cldaero_update( &
166 218880 : state, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, xlwc, &
167 145920 : delso4_hprxn, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, &
168 145920 : aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d, aqso4_o3_3d)
169 :
170 : use aerosol_properties_mod, only: aero_name_len
171 : use physics_types, only: physics_state
172 : use carma_intr, only: carma_get_group_by_name, carma_get_dry_radius
173 :
174 : ! args
175 :
176 : type(physics_state), intent(in) :: state ! Physics state variables
177 :
178 : integer, intent(in) :: ncol
179 : integer, intent(in) :: lchnk ! chunk id
180 : integer, intent(in) :: loffset
181 :
182 : real(r8), intent(in) :: dtime ! time step (sec)
183 :
184 : real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu )
185 : real(r8), intent(in) :: pdel(:,:)
186 : real(r8), intent(in) :: press(:,:)
187 : real(r8), intent(in) :: tfld(:,:)
188 :
189 : real(r8), intent(in) :: cldnum(:,:)
190 : real(r8), intent(in) :: cldfrc(:,:)
191 : real(r8), intent(in) :: cfact(:,:)
192 : real(r8), intent(in) :: xlwc(:,:)
193 :
194 : real(r8), intent(in) :: delso4_hprxn(:,:)
195 : real(r8), intent(in) :: xh2so4(:,:)
196 : real(r8), intent(in) :: xso4(:,:)
197 : real(r8), intent(in) :: xso4_init(:,:)
198 : real(r8), intent(in) :: nh3g(:,:)
199 : real(r8), intent(in) :: hno3g(:,:)
200 : real(r8), intent(in) :: xnh3(:,:)
201 : real(r8), intent(in) :: xhno3(:,:)
202 : real(r8), intent(in) :: xnh4c(:,:)
203 : real(r8), intent(in) :: xmsa(:,:)
204 : real(r8), intent(in) :: xso2(:,:)
205 : real(r8), intent(in) :: xh2o2(:,:)
206 : real(r8), intent(in) :: xno3c(:,:)
207 :
208 : real(r8), intent(inout) :: qcw(:,:,:) ! cloud-borne aerosol (vmr) vmrcw(ncol,pver,ncnst_tot)
209 : real(r8), intent(inout) :: qin(:,:,:) ! xported species ( vmr )
210 :
211 : real(r8), intent(out) :: aqso4(:,:) ! aqueous phase chemistry
212 : real(r8), intent(out) :: aqh2so4(:,:) ! aqueous phase chemistry
213 : real(r8), intent(out) :: aqso4_h2o2(:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2)
214 : real(r8), intent(out) :: aqso4_o3(:) ! SO4 aqueous phase chemistry due to O3 (kg/m2)
215 : real(r8), intent(out), optional :: aqso4_h2o2_3d(:,:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2)
216 : real(r8), intent(out), optional :: aqso4_o3_3d(:,:) ! SO4 aqueous phase chemistry due to O3 (kg/m2)
217 :
218 : ! local vars ...
219 : real(r8) :: dryr(pcols,pver) ! CARMA dry radius in cm
220 : real(r8) :: rho(pcols,pver) !
221 145920 : real(r8) :: dryr_n(nbins,ncol,pver) ! CARMA dry radius in cm
222 145920 : real(r8) :: dqdt_aqso4(ncol,pver,ncnst_tot), &
223 145920 : dqdt_aqh2so4(ncol,pver,ncnst_tot), &
224 145920 : dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver)
225 :
226 145920 : real(r8) :: faqgain_so4(nbins)
227 72960 : real(r8) :: wt_mass(nbins)
228 :
229 : real(r8) :: delso4_o3rxn, &
230 : dso4dt_aqrxn, dso4dt_hprxn, &
231 : dso4dt_gasuptk, dmsadt_gasuptk_toso4, &
232 : dqdt_aq, dqdt_wr, dqdt
233 :
234 : real(r8) :: fwetrem, uptkrate
235 :
236 : integer :: l, n, mm
237 : integer :: ntot_msa_c
238 :
239 : integer :: i,k
240 : real(r8) :: xl
241 : real(r8) :: wt_sum
242 : real(r8) :: specmw_so4_amode
243 :
244 : character(len=32) :: spectype
245 :
246 : character(len=*), parameter :: subname = 'sox_cldaero_update'
247 : character(len=aero_name_len) :: bin_name, shortname
248 : integer :: igroup, ibin, rc, nchr
249 :
250 : ! make sure dqdt is zero initially, for budgets
251 11021410560 : dqdt_aqso4(:,:,:) = 0.0_r8
252 11021410560 : dqdt_aqh2so4(:,:,:) = 0.0_r8
253 78723840 : dqdt_aqhprxn(:,:) = 0.0_r8
254 78723840 : dqdt_aqo3rxn(:,:) = 0.0_r8
255 3020471040 : dryr_n(:,:,:) = 0.0_r8
256 :
257 : ntot_msa_c = 0.0_r8
258 157374720 : aqso4 = 0.0_r8
259 157374720 : aqh2so4 = 0.0_r8
260 1123584 : aqso4_h2o2 = 0.0_r8
261 1123584 : aqso4_o3 = 0.0_r8
262 :
263 2991360 : do n = 1, nbins
264 2918400 : call rad_cnst_get_info_by_bin(0, n, nspec=nspec(n), bin_name=bin_name)
265 :
266 :
267 2918400 : nchr = len_trim(bin_name)-2
268 2918400 : shortname = bin_name(:nchr)
269 :
270 2918400 : call carma_get_group_by_name(shortname, igroup, rc)
271 2918400 : if (rc/=0) then
272 0 : call endrun(subname//': ERROR in carma_get_group_by_name')
273 : end if
274 :
275 2918400 : read(bin_name(nchr+1:),*) ibin
276 :
277 2918400 : call carma_get_dry_radius(state, igroup, ibin, dryr, rho, rc)
278 2918400 : if (rc/=0) then
279 0 : call endrun(subname//': ERROR in carma_get_dry_radius')
280 : end if
281 :
282 3148953600 : dryr(:ncol,:) = dryr(:ncol,:)*1.e2_r8 ! cm
283 :
284 5909760 : if (index(bin_name,'MXAER')>0) then
285 1574476800 : dryr_n(n,:ncol,:) = dryr(:ncol,:)
286 : end if
287 : end do
288 :
289 5180160 : lev_loop: do k = 1,pver
290 78723840 : col_loop: do i = 1,ncol
291 78650880 : cloud: if (cldfrc(i,k) >= 1.0e-5_r8) then
292 10467718 : xl = xlwc(i,k)
293 :
294 10467718 : if (xl .ge. 1.e-8_r8) then !! when cloud is present
295 :
296 1837329 : delso4_o3rxn = xso4(i,k) - xso4_init(i,k)
297 :
298 : ! the factors are proportional to the activated particle MR for each
299 : ! bin, which is the MR of cloud drops "associated with" the mode
300 : ! thus we are assuming the cloud drop size is independent of the
301 : ! associated aerosol mode properties (i.e., drops associated with
302 : ! Aitken and coarse sea-salt particles are same size)
303 : !
304 : ! qnum_c(n) = activated particle number MR for mode n (these are just
305 : ! used for partitioning among modes, so don't need to divide by cldfrc)
306 :
307 : !faqgain_so4(n) = fraction of total so4_c gain going to mode n
308 1837329 : wt_sum = 0._r8
309 75330489 : wt_mass(:) = 0._r8
310 75330489 : faqgain_so4(:) = 0.0_r8
311 75330489 : do n = 1, nbins
312 75330489 : if (dryr_n(n,i,k) > 0._r8) then
313 36746580 : wt_mass(n) = delso4_o3rxn / dryr_n(n,i,k) / dryr_n(n,i,k)
314 36746580 : wt_sum = wt_sum + wt_mass(n)
315 : end if
316 : end do
317 75330489 : do n = 1, nbins
318 75330489 : if (wt_mass(n) > 0._r8) then
319 36746580 : faqgain_so4(n) = wt_mass(n)/wt_sum
320 : end if
321 : end do
322 :
323 1837329 : uptkrate = cldaero_uptakerate( xl, cldnum(i,k), cfact(i,k), cldfrc(i,k), tfld(i,k), press(i,k) )
324 : ! average uptake rate over dtime
325 1837329 : uptkrate = (1.0_r8 - exp(-min(100._r8,dtime*uptkrate))) / dtime
326 :
327 1837329 : dso4dt_gasuptk = xh2so4(i,k) * uptkrate
328 :
329 : ! if no modes have msa aerosol, then "rename" scavenged msa gas to so4
330 1837329 : dmsadt_gasuptk_toso4 = 0.0_r8
331 :
332 : !-----------------------------------------------------------------------
333 : ! now compute TMR tendencies
334 : ! this includes the above aqueous so2 chemistry AND
335 : ! the uptake of highly soluble aerosol precursor gases (h2so4, msa, ...)
336 : ! AND the wetremoval of dissolved, unreacted so2 and h2o2
337 :
338 1837329 : dso4dt_aqrxn = (delso4_o3rxn + delso4_hprxn(i,k)) / dtime
339 1837329 : dso4dt_hprxn = delso4_hprxn(i,k) / dtime
340 : !write(iulog,*) 'dso4dt_aqrxn ',dso4dt_aqrxn
341 :
342 : ! fwetrem = fraction of in-cloud-water material that is wet removed
343 : ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*clwlrat(i,k)))) )
344 1837329 : fwetrem = 0.0_r8 ! don't have so4 & msa wet removal here
345 :
346 : ! compute TMR tendencies for so4, not done currently for msa aerosol-in-cloud-water
347 75330489 : do n = 1, nbins
348 332556549 : do l = 1, nspec(n)
349 257226060 : mm = bin_idx(n, l)
350 257226060 : call rad_cnst_get_bin_props_by_idx(0, n, l,spectype=spectype)
351 330719220 : if (trim(spectype) == 'sulfate') then
352 73493160 : if (faqgain_so4(n) .gt. 0.0_r8) then
353 36746580 : dqdt_aqso4(i,k,mm) = faqgain_so4(n)*dso4dt_aqrxn*cldfrc(i,k)
354 :
355 : dqdt_aqh2so4(i,k,mm) = faqgain_so4(n)* &
356 36746580 : (dso4dt_gasuptk + dmsadt_gasuptk_toso4)*cldfrc(i,k)
357 36746580 : dqdt_aq = dqdt_aqso4(i,k,mm) + dqdt_aqh2so4(i,k,mm)
358 36746580 : dqdt_wr = -fwetrem*dqdt_aq
359 36746580 : dqdt= dqdt_aq + dqdt_wr
360 : !write(iulog,*) 'qcw(i,k,mm) before ', m, qcw(i,k,mm)
361 36746580 : qcw(i,k,mm) = qcw(i,k,mm) + dqdt*dtime
362 : !write(iulog,*) 'qcw(i,k,mm) after', m, qcw(i,k,mm)
363 : end if
364 : end if
365 : end do
366 : end do
367 :
368 :
369 : ! For gas species, tendency includes
370 : ! reactive uptake to cloud water that essentially transforms the gas to
371 : ! a different species. Wet removal associated with this is applied
372 : ! to the "new" species (e.g., so4_c) rather than to the gas.
373 : ! wet removal of the unreacted gas that is dissolved in cloud water.
374 : ! Need to multiply both these parts by cldfrc
375 :
376 : ! h2so4 (g) & msa (g)
377 1837329 : qin(i,k,id_h2so4) = qin(i,k,id_h2so4) - dso4dt_gasuptk * dtime * cldfrc(i,k)
378 :
379 : ! so2 -- the first order loss rate for so2 is frso2_c*clwlrat(i,k)
380 : ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frso2_c*clwlrat(i,k)))) )
381 1837329 : fwetrem = 0.0_r8 ! don't include so2 wet removal here
382 :
383 1837329 : dqdt_wr = -fwetrem*xso2(i,k)/dtime*cldfrc(i,k)
384 1837329 : dqdt_aq = -dso4dt_aqrxn*cldfrc(i,k)
385 1837329 : dqdt = dqdt_aq + dqdt_wr
386 1837329 : qin(i,k,id_so2) = qin(i,k,id_so2) + dqdt * dtime
387 1837329 : qin(i,k,id_so2) = MAX( qin(i,k,id_so2), small_value )
388 :
389 : ! h2o2 -- the first order loss rate for h2o2 is frh2o2_c*clwlrat(i,k)
390 : ! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frh2o2_c*clwlrat(i,k)))) )
391 1837329 : fwetrem = 0.0_r8 ! don't include h2o2 wet removal here
392 :
393 1837329 : dqdt_wr = -fwetrem*xh2o2(i,k)/dtime*cldfrc(i,k)
394 1837329 : dqdt_aq = -dso4dt_hprxn*cldfrc(i,k)
395 1837329 : dqdt = dqdt_aq + dqdt_wr
396 1837329 : qin(i,k,id_h2o2) = qin(i,k,id_h2o2) + dqdt * dtime
397 1837329 : qin(i,k,id_h2o2) = MAX( qin(i,k,id_h2o2), small_value )
398 :
399 : ! for SO4 from H2O2/O3 budgets
400 1837329 : dqdt_aqhprxn(i,k) = dso4dt_hprxn*cldfrc(i,k)
401 1837329 : dqdt_aqo3rxn(i,k) = (dso4dt_aqrxn - dso4dt_hprxn)*cldfrc(i,k)
402 :
403 : endif !! when cloud is present
404 : endif cloud
405 : enddo col_loop
406 : enddo lev_loop
407 :
408 : !==============================================================
409 : ! ... Update the mixing ratios
410 : !==============================================================
411 :
412 : ! diagnostics
413 :
414 72960 : specmw_so4_amode = 96.0_r8
415 2991360 : do n = 1, nbins
416 : ! while looking through all species, only dqdt_aqso4 from sulfates is gt zero
417 13205760 : do l = 1, nspec(n)
418 10214400 : mm = bin_idx(n, l)
419 157301760 : aqso4(:,n)=0._r8
420 725222400 : do k=1,pver
421 11021337600 : do i=1,ncol
422 30888345600 : aqso4(i,n)=aqso4(i,n)+dqdt_aqso4(i,k,mm)*specmw_so4_amode/mbar(i,k) &
423 41899468800 : *pdel(i,k)/gravit ! kg/m2/s
424 : enddo
425 : enddo
426 :
427 157301760 : aqh2so4(:,n)=0._r8
428 728140800 : do k=1,pver
429 11021337600 : do i=1,ncol
430 30888345600 : aqh2so4(i,n)=aqh2so4(i,n)+dqdt_aqh2so4(i,k,mm)*specmw_so4_amode/mbar(i,k) &
431 41899468800 : *pdel(i,k)/gravit ! kg/m2/s
432 : enddo
433 : enddo
434 : end do
435 : end do
436 :
437 1123584 : aqso4_h2o2(:) = 0._r8
438 5180160 : do k=1,pver
439 78723840 : do i=1,ncol
440 147087360 : aqso4_h2o2(i)=aqso4_h2o2(i)+dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) &
441 225738240 : *pdel(i,k)/gravit ! kg SO4 /m2/s
442 : enddo
443 : enddo
444 :
445 72960 : if (present(aqso4_h2o2_3d)) then
446 0 : aqso4_h2o2_3d(:,:) = 0._r8
447 0 : do k=1,pver
448 0 : do i=1,ncol
449 0 : aqso4_h2o2_3d(i,k)=dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) &
450 0 : *pdel(i,k)/gravit ! kg SO4 /m2/s
451 : enddo
452 : enddo
453 : end if
454 :
455 1123584 : aqso4_o3(:)=0._r8
456 5180160 : do k=1,pver
457 78723840 : do i=1,ncol
458 147087360 : aqso4_o3(i)=aqso4_o3(i)+dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) &
459 225738240 : *pdel(i,k)/gravit ! kg SO4 /m2/s
460 : enddo
461 : enddo
462 :
463 72960 : if (present(aqso4_o3_3d)) then
464 0 : aqso4_o3_3d(:,:)=0._r8
465 0 : do k=1,pver
466 0 : do i=1,ncol
467 0 : aqso4_o3_3d(i,k)=dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) &
468 0 : *pdel(i,k)/gravit ! kg SO4 /m2/s
469 : enddo
470 : enddo
471 : end if
472 :
473 72960 : end subroutine sox_cldaero_update
474 :
475 : !----------------------------------------------------------------------------------
476 : !----------------------------------------------------------------------------------
477 72960 : subroutine sox_cldaero_destroy_obj( conc_obj )
478 : type(cldaero_conc_t), pointer :: conc_obj
479 :
480 72960 : call cldaero_deallocate( conc_obj )
481 :
482 72960 : end subroutine sox_cldaero_destroy_obj
483 :
484 : end module sox_cldaero_mod
|
