| blob | a6c420d1c7826488f458c382cc62696b160985d5 |
1 #define WRF_PORT
2 #define MODAL_AERO
3 ! module_cam_mam_gasaerexch.F
4 ! adapted from cam3 modal_aero_gasaerexch.F90 by r.c.easter, august 2010
5 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
6 !
7 ! 2010-08-04 rce notes
8 ! > change pcnstxx to a subr argument because gas_pcnst is
9 ! not a constant/parameter in wrf-chem
10 !--------------------------------------------------------------
12 ! modal_aero_gasaerexch.F90
15 !----------------------------------------------------------------------
16 !----------------------------------------------------------------------
17 !BOP
18 !
19 ! !MODULE: modal_aero_gasaerexch --- does modal aerosol gas-aerosol exchange
20 !
21 ! !INTERFACE:
22 module modal_aero_gasaerexch
24 ! !USES:
25 use shr_kind_mod, only: r8 => shr_kind_r8
26 #ifndef WRF_PORT
27 use chem_mods, only: gas_pcnst
28 #else
29 use module_cam_support, only: gas_pcnst => gas_pcnst_modal_aero
30 #endif
31 use modal_aero_data, only: maxd_aspectype
33 implicit none
34 private
35 save
37 ! !PUBLIC MEMBER FUNCTIONS:
38 public modal_aero_gasaerexch_sub, modal_aero_gasaerexch_init
40 ! !PUBLIC DATA MEMBERS:
41 #ifndef WRF_PORT
42 integer, parameter :: pcnstxx = gas_pcnst
43 #endif
44 integer, parameter, public :: maxspec_pcage = maxd_aspectype
46 integer, public :: npair_pcage
47 integer, public :: modefrm_pcage
48 integer, public :: modetoo_pcage
49 integer, public :: nspecfrm_pcage
50 integer, public :: lspecfrm_pcage(maxspec_pcage)
51 integer, public :: lspectoo_pcage(maxspec_pcage)
54 ! real(r8), parameter, public :: n_so4_monolayers_pcage = 1.0_r8
55 real(r8), parameter, public :: n_so4_monolayers_pcage = 3.0_r8
56 ! number of so4(+nh4) monolayers needed to "age" a carbon particle
58 real(r8), parameter, public :: &
59 dr_so4_monolayers_pcage = n_so4_monolayers_pcage * 4.76e-10
60 ! thickness of the so4 monolayers (m)
61 ! for so4(+nh4), use bi-sulfate mw and 1.77 g/cm3,
62 ! --> 1 mol so4(+nh4) = 65 cm^3 --> 1 molecule = (4.76e-10 m)^3
63 ! aging criterion is approximate so do not try to distinguish
64 ! sulfuric acid, bisulfate, ammonium sulfate
66 real(r8), save, public :: soa_equivso4_factor
67 ! this factor converts an soa volume to a volume of so4(+nh4)
68 ! having same hygroscopicity as the soa
71 ! !DESCRIPTION: This module implements ...
72 !
73 ! !REVISION HISTORY:
74 !
75 ! RCE 07.04.13: Adapted from MIRAGE2 code
76 !
77 !EOP
78 !----------------------------------------------------------------------
79 !BOC
81 ! list private module data here
83 !EOC
84 !----------------------------------------------------------------------
87 contains
90 !----------------------------------------------------------------------
91 !----------------------------------------------------------------------
92 !BOP
93 ! !ROUTINE: modal_aero_gasaerexch_sub --- ...
94 !
95 ! !INTERFACE:
96 subroutine modal_aero_gasaerexch_sub( &
97 lchnk, ncol, nstep, &
98 loffset, deltat, &
99 latndx, lonndx, &
100 t, pmid, pdel, &
101 q, qqcw, &
102 dqdt_other, dqqcwdt_other, &
103 dgncur_a, dgncur_awet &
104 #ifdef WRF_PORT
105 , pcnstxx &
106 #endif
107 )
109 ! !USES:
110 use modal_aero_data
111 use modal_aero_rename, only: modal_aero_rename_sub
112 #ifndef WRF_PORT
113 use cam_history, only: outfld, fieldname_len
114 use chem_mods, only: adv_mass
115 use constituents, only: pcnst, cnst_name, cnst_get_ind
116 use mo_tracname, only: solsym
117 #endif
118 use physconst, only: gravit, mwdry, rair
119 #ifndef WRF_PORT
120 use ppgrid, only: pcols, pver
122 use abortutils, only : endrun
123 use spmd_utils, only : iam, masterproc
124 #else
125 use module_cam_support, only: outfld, fieldname_len, pcnst => pcnst_runtime, &
126 pcols, pver, endrun, iam, masterproc
127 use constituents, only: cnst_name, cnst_get_ind
128 use module_data_cam_mam_asect, only: adv_mass => mw_q_mo_array
129 #endif
133 implicit none
135 ! !PARAMETERS:
136 #ifdef WRF_PORT
137 integer, intent(in) :: pcnstxx ! number of species in "incoming q array"
138 #endif
139 integer, intent(in) :: lchnk ! chunk identifier
140 integer, intent(in) :: ncol ! number of atmospheric column
141 integer, intent(in) :: nstep ! model time-step number
142 integer, intent(in) :: loffset ! offset applied to modal aero "ptrs"
143 integer, intent(in) :: latndx(pcols), lonndx(pcols)
144 real(r8), intent(in) :: deltat ! time step (s)
146 real(r8), intent(inout) :: q(ncol,pver,pcnstxx) ! tracer mixing ratio (TMR) array
147 ! *** MUST BE #/kmol-air for number
148 ! *** MUST BE mol/mol-air for mass
149 ! *** NOTE ncol dimension
150 real(r8), intent(inout) :: qqcw(ncol,pver,pcnstxx)
151 ! like q but for cloud-borner tracers
152 real(r8), intent(in) :: dqdt_other(ncol,pver,pcnstxx)
153 ! TMR tendency from other continuous
154 ! growth processes (aqchem, soa??)
155 ! *** NOTE ncol dimension
156 real(r8), intent(in) :: dqqcwdt_other(ncol,pver,pcnstxx)
157 ! like dqdt_other but for cloud-borner tracers
158 real(r8), intent(in) :: t(pcols,pver) ! temperature at model levels (K)
159 real(r8), intent(in) :: pmid(pcols,pver) ! pressure at model levels (Pa)
160 real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickness of levels (Pa)
161 real(r8), intent(in) :: dgncur_a(pcols,pver,ntot_amode)
162 real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)
163 ! dry & wet geo. mean dia. (m) of number distrib.
165 ! !DESCRIPTION:
166 ! computes TMR (tracer mixing ratio) tendencies for gas condensation
167 ! onto aerosol particles
168 !
169 ! this version does condensation of H2SO4, NH3, and MSA, both treated as
170 ! completely non-volatile (gas --> aerosol, but no aerosol --> gas)
171 ! gas H2SO4 goes to aerosol SO4
172 ! gas MSA (if present) goes to aerosol SO4
173 ! aerosol MSA is not distinguished from aerosol SO4
174 ! gas NH3 (if present) goes to aerosol NH4
175 ! if gas NH3 is not present, then ????
176 !
177 !
178 ! !REVISION HISTORY:
179 ! RCE 07.04.13: Adapted from MIRAGE2 code
180 !
181 !EOP
182 !----------------------------------------------------------------------
183 !BOC
185 ! local variables
186 integer, parameter :: jsrflx_gaexch = 1
187 integer, parameter :: jsrflx_rename = 2
188 integer, parameter :: ldiag1=-1, ldiag2=-1, ldiag3=-1, ldiag4=-1
189 integer, parameter :: method_soa = 2
190 ! method_soa=0 is no uptake
191 ! method_soa=1 is irreversible uptake done like h2so4 uptake
192 ! method_soa=2 is reversible uptake using subr modal_aero_soaexch
194 integer :: i, icol_diag, iq
195 integer :: idiagss
196 integer :: ido_so4a(ntot_amode), ido_nh4a(ntot_amode), ido_soaa(ntot_amode)
197 integer :: jac, jsrf
198 integer :: k
199 integer :: l, l2, lb, lsfrm, lstoo
200 integer :: l_so4g, l_nh4g, l_msag, l_soag
201 integer :: n, niter, niter_max, ntot_soamode
203 logical :: is_dorename_atik, dorename_atik(ncol,pver)
205 character(len=fieldname_len+3) :: fieldname
207 real (r8) :: avg_uprt_nh4, avg_uprt_so4, avg_uprt_soa
208 real (r8) :: deltatxx
209 real (r8) :: dqdt_nh4(ntot_amode), dqdt_so4(ntot_amode), dqdt_soa(ntot_amode)
210 real (r8) :: fac_m2v_nh4, fac_m2v_so4, fac_m2v_soa
211 real (r8) :: fac_m2v_pcarbon(maxd_aspectype)
212 real (r8) :: fac_volsfc_pcarbon
213 real (r8) :: fgain_nh4(ntot_amode), fgain_so4(ntot_amode), fgain_soa(ntot_amode)
214 real (r8) :: g0_soa
215 real (r8) :: pdel_fac
216 real (r8) :: qmax_nh4, qnew_nh4, qnew_so4
217 real (r8) :: qold_nh4(ntot_amode), qold_so4(ntot_amode)
218 real (r8) :: qold_soa(ntot_amode), qold_poa(ntot_amode)
219 real (r8) :: sum_dqdt_msa, sum_dqdt_so4, sum_dqdt_soa
220 real (r8) :: sum_dqdt_nh4, sum_dqdt_nh4_b
221 real (r8) :: sum_uprt_msa, sum_uprt_nh4, sum_uprt_so4, sum_uprt_soa
222 real (r8) :: tmp1, tmp2
223 real (r8) :: uptkrate(ntot_amode,pcols,pver)
224 real (r8) :: uptkratebb(ntot_amode), uptkrate_soa(ntot_amode)
225 ! gas-to-aerosol mass transfer rates (1/s)
226 real (r8) :: vol_core, vol_shell
227 real (r8) :: xferfrac_pcage, xferfrac_max
228 real (r8) :: xferrate
230 logical :: do_msag ! true if msa gas is a species
231 logical :: do_nh4g ! true if nh3 gas is a species
232 logical :: do_soag ! true if soa gas is a species
234 logical :: dotend(pcnstxx) ! identifies species directly involved in
235 ! gas-aerosol exchange (gas condensation)
236 logical :: dotendqqcw(pcnstxx) ! like dotend but for cloud-borner tracers
237 logical :: dotendrn(pcnstxx), dotendqqcwrn(pcnstxx)
238 ! identifies species involved in renaming
239 ! after "continuous growth"
240 ! (gas-aerosol exchange and aqchem)
242 integer, parameter :: nsrflx = 2 ! last dimension of qsrflx
243 real(r8) :: dqdt(ncol,pver,pcnstxx) ! TMR "delta q" array - NOTE dims
244 real(r8) :: dqqcwdt(ncol,pver,pcnstxx) ! like dqdt but for cloud-borner tracers
245 real(r8) :: qsrflx(pcols,pcnstxx,nsrflx)
246 ! process-specific column tracer tendencies
247 ! (1=renaming, 2=gas condensation)
248 real(r8) :: qqcwsrflx(pcols,pcnstxx,nsrflx)
250 ! following only needed for diagnostics
251 real(r8) :: qold(ncol,pver,pcnstxx) ! NOTE dims
252 real(r8) :: qnew(ncol,pver,pcnstxx) ! NOTE dims
253 real(r8) :: qdel(ncol,pver,pcnstxx) ! NOTE dims
254 real(r8) :: dumavec(1000), dumbvec(1000), dumcvec(1000)
255 real(r8) :: qqcwold(ncol,pver,pcnstxx)
256 real(r8) :: dqdtsv1(ncol,pver,pcnstxx)
257 real(r8) :: dqqcwdtsv1(ncol,pver,pcnstxx)
260 !----------------------------------------------------------------------
263 icol_diag = -1
264 if (ldiag1 > 0) then
265 if (nstep < 3) then
266 do i = 1, ncol
267 if ((latndx(i) == 23) .and. (lonndx(i) == 37)) icol_diag = i
268 end do
269 end if
270 end if
273 ! set gas species indices
274 call cnst_get_ind( 'H2SO4', l_so4g, .false. )
275 call cnst_get_ind( 'NH3', l_nh4g, .false. )
276 call cnst_get_ind( 'MSA', l_msag, .false. )
277 call cnst_get_ind( 'SOAG', l_soag, .false. )
278 l_so4g = l_so4g - loffset
279 l_nh4g = l_nh4g - loffset
280 l_msag = l_msag - loffset
281 l_soag = l_soag - loffset
282 if ((l_so4g <= 0) .or. (l_so4g > pcnstxx)) then
283 write( *, '(/a/a,2i7)' ) &
284 '*** modal_aero_gasaerexch_sub -- cannot find H2SO4 species', &
285 ' l_so4g, loffset =', l_so4g, loffset
286 call endrun( 'modal_aero_gasaerexch_sub error' )
287 end if
288 do_nh4g = .false.
289 do_msag = .false.
290 if ((l_nh4g > 0) .and. (l_nh4g <= pcnstxx)) do_nh4g = .true.
291 if ((l_msag > 0) .and. (l_msag <= pcnstxx)) do_msag = .true.
292 do_soag = .false.
293 if ((method_soa == 1) .or. (method_soa == 2)) then
294 if ((l_soag > 0) .and. (l_soag <= pcnstxx)) do_soag = .true.
295 else if (method_soa /= 0) then
296 write(*,'(/a,1x,i10)') '*** modal_aero_gasaerexch_sub - bad method_soa =', method_soa
297 call endrun( 'modal_aero_gasaerexch_sub error' )
298 end if
300 ! set tendency flags
301 dotend(:) = .false.
302 dotendqqcw(:) = .false.
303 ido_so4a(:) = 0
304 ido_nh4a(:) = 0
305 ido_soaa(:) = 0
307 dotend(l_so4g) = .true.
308 if ( do_nh4g ) dotend(l_nh4g) = .true.
309 if ( do_msag ) dotend(l_msag) = .true.
310 if ( do_soag ) dotend(l_soag) = .true.
311 ntot_soamode = 0
312 do n = 1, ntot_amode
313 l = lptr_so4_a_amode(n)-loffset
314 if ((l > 0) .and. (l <= pcnstxx)) then
315 dotend(l) = .true.
316 ido_so4a(n) = 1
317 if ( do_nh4g ) then
318 l = lptr_nh4_a_amode(n)-loffset
319 if ((l > 0) .and. (l <= pcnstxx)) then
320 dotend(l) = .true.
321 ido_nh4a(n) = 1
322 end if
323 end if
324 end if
325 if ( do_soag ) then
326 l = lptr_soa_a_amode(n)-loffset
327 if ((l > 0) .and. (l <= pcnstxx)) then
328 dotend(l) = .true.
329 ido_soaa(n) = 1
330 ntot_soamode = n
331 end if
332 end if
333 end do
334 if ( do_soag ) ntot_soamode = max( ntot_soamode, modefrm_pcage )
336 if (modefrm_pcage > 0) then
337 ido_so4a(modefrm_pcage) = 2
338 if (ido_nh4a(modetoo_pcage) == 1) ido_nh4a(modefrm_pcage) = 2
339 if (ido_soaa(modetoo_pcage) == 1) ido_soaa(modefrm_pcage) = 2
340 do iq = 1, nspecfrm_pcage
341 lsfrm = lspecfrm_pcage(iq)-loffset
342 lstoo = lspectoo_pcage(iq)-loffset
343 if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
344 dotend(lsfrm) = .true.
345 if ((lstoo > 0) .and. (lstoo <= pcnst)) then
346 dotend(lstoo) = .true.
347 end if
348 end if
349 end do
352 fac_m2v_so4 = specmw_so4_amode / specdens_so4_amode
353 fac_m2v_nh4 = specmw_nh4_amode / specdens_nh4_amode
354 fac_m2v_soa = specmw_soa_amode / specdens_soa_amode
355 fac_m2v_pcarbon(:) = 0.0
356 n = modeptr_pcarbon
357 do l = 1, nspec_amode(n)
358 l2 = lspectype_amode(l,n)
359 ! fac_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
360 ! [m3-AP/kmol-AP] = [kg-AP/kmol-AP] / [kg-AP/m3-AP]
361 fac_m2v_pcarbon(l) = specmw_amode(l2) / specdens_amode(l2)
362 end do
363 fac_volsfc_pcarbon = exp( 2.5*(alnsg_amode(n)**2) )
364 xferfrac_max = 1.0_r8 - 10.0_r8*epsilon(1.0_r8) ! 1-eps
365 end if
368 ! zero out tendencies and other
369 dqdt(:,:,:) = 0.0
370 dqqcwdt(:,:,:) = 0.0
371 qsrflx(:,:,:) = 0.0
372 qqcwsrflx(:,:,:) = 0.0
374 ! compute gas-to-aerosol mass transfer rates
375 call gas_aer_uptkrates( ncol, loffset, &
376 q, t, pmid, &
377 dgncur_awet, uptkrate &
378 #ifdef WRF_PORT
379 , pcnstxx &
380 #endif
381 )
384 ! use this for tendency calcs to avoid generating very small negative values
385 deltatxx = deltat * (1.0d0 + 1.0d-15)
388 jsrf = jsrflx_gaexch
389 do k=1,pver
390 do i=1,ncol
392 ! fgain_so4(n) = fraction of total h2so4 uptake going to mode n
393 ! fgain_nh4(n) = fraction of total nh3 uptake going to mode n
394 sum_uprt_so4 = 0.0
395 sum_uprt_nh4 = 0.0
396 sum_uprt_soa = 0.0
397 do n = 1, ntot_amode
398 uptkratebb(n) = uptkrate(n,i,k)
399 if (ido_so4a(n) > 0) then
400 fgain_so4(n) = uptkratebb(n)
401 sum_uprt_so4 = sum_uprt_so4 + fgain_so4(n)
402 if (ido_so4a(n) == 1) then
403 qold_so4(n) = q(i,k,lptr_so4_a_amode(n)-loffset)
404 else
405 qold_so4(n) = 0.0
406 end if
407 else
408 fgain_so4(n) = 0.0
409 qold_so4(n) = 0.0
410 end if
412 if (ido_nh4a(n) > 0) then
413 ! 2.08 factor is for gas diffusivity (nh3/h2so4)
414 ! differences in fuch-sutugin and accom coef ignored
415 fgain_nh4(n) = uptkratebb(n)*2.08
416 sum_uprt_nh4 = sum_uprt_nh4 + fgain_nh4(n)
417 if (ido_nh4a(n) == 1) then
418 qold_nh4(n) = q(i,k,lptr_nh4_a_amode(n)-loffset)
419 else
420 qold_nh4(n) = 0.0
421 end if
422 else
423 fgain_nh4(n) = 0.0
424 qold_nh4(n) = 0.0
425 end if
427 if (ido_soaa(n) > 0) then
428 ! 0.81 factor is for gas diffusivity (soa/h2so4)
429 ! (differences in fuch-sutugin and accom coef ignored)
430 fgain_soa(n) = uptkratebb(n)*0.81
431 sum_uprt_soa = sum_uprt_soa + fgain_soa(n)
432 if (ido_soaa(n) == 1) then
433 qold_soa(n) = q(i,k,lptr_soa_a_amode(n)-loffset)
434 l = lptr_pom_a_amode(n)-loffset
435 if (l > 0) then
436 qold_poa(n) = q(i,k,l)
437 else
438 qold_poa(n) = 0.0
439 end if
440 else
441 qold_soa(n) = 0.0
442 qold_poa(n) = 0.0
443 end if
444 else
445 fgain_soa(n) = 0.0
446 qold_soa(n) = 0.0
447 qold_poa(n) = 0.0
448 end if
449 uptkrate_soa(n) = fgain_soa(n)
450 end do
452 if (sum_uprt_so4 > 0.0) then
453 do n = 1, ntot_amode
454 fgain_so4(n) = fgain_so4(n) / sum_uprt_so4
455 end do
456 end if
457 ! at this point (sum_uprt_so4 <= 0.0) only when all the fgain_so4 are zero
458 if (sum_uprt_nh4 > 0.0) then
459 do n = 1, ntot_amode
460 fgain_nh4(n) = fgain_nh4(n) / sum_uprt_nh4
461 end do
462 end if
463 if (sum_uprt_soa > 0.0) then
464 do n = 1, ntot_amode
465 fgain_soa(n) = fgain_soa(n) / sum_uprt_soa
466 end do
467 end if
469 ! uptake amount (fraction of gas uptaken) over deltat
470 avg_uprt_so4 = (1.0 - exp(-deltatxx*sum_uprt_so4))/deltatxx
471 avg_uprt_nh4 = (1.0 - exp(-deltatxx*sum_uprt_nh4))/deltatxx
472 avg_uprt_soa = (1.0 - exp(-deltatxx*sum_uprt_soa))/deltatxx
474 ! sum_dqdt_so4 = so4_a tendency from h2so4 gas uptake (mol/mol/s)
475 ! sum_dqdt_msa = msa_a tendency from msa gas uptake (mol/mol/s)
476 ! sum_dqdt_nh4 = nh4_a tendency from nh3 gas uptake (mol/mol/s)
477 ! sum_dqdt_soa = soa_a tendency from soa gas uptake (mol/mol/s)
478 sum_dqdt_so4 = q(i,k,l_so4g) * avg_uprt_so4
479 if ( do_msag ) then
480 sum_dqdt_msa = q(i,k,l_msag) * avg_uprt_so4
481 else
482 sum_dqdt_msa = 0.0
483 end if
484 if ( do_nh4g ) then
485 sum_dqdt_nh4 = q(i,k,l_nh4g) * avg_uprt_nh4
486 else
487 sum_dqdt_nh4 = 0.0
488 end if
489 if ( do_soag ) then
490 sum_dqdt_soa = q(i,k,l_soag) * avg_uprt_soa
491 else
492 sum_dqdt_soa = 0.0
493 end if
495 ! compute TMR tendencies for so4, nh4, msa interstial aerosol
496 ! due to simple gas uptake
497 pdel_fac = pdel(i,k)/gravit
498 sum_dqdt_nh4_b = 0.0
499 do n = 1, ntot_amode
500 dqdt_so4(n) = fgain_so4(n)*(sum_dqdt_so4 + sum_dqdt_msa)
502 if ( do_nh4g ) then
503 dqdt_nh4(n) = fgain_nh4(n)*sum_dqdt_nh4
504 qnew_nh4 = qold_nh4(n) + dqdt_nh4(n)*deltat
505 qnew_so4 = qold_so4(n) + dqdt_so4(n)*deltat
506 qmax_nh4 = 2.0*qnew_so4
507 if (qnew_nh4 > qmax_nh4) then
508 dqdt_nh4(n) = (qmax_nh4 - qold_nh4(n))/deltatxx
509 end if
510 sum_dqdt_nh4_b = sum_dqdt_nh4_b + dqdt_nh4(n)
511 end if
512 end do
514 if (( do_soag ) .and. (method_soa > 1)) then
515 ! compute TMR tendencies for soag and soa interstial aerosol
516 ! using soa parameterization
517 niter_max = 1000
518 dqdt_soa(:) = 0.0
520 ! idiagss = 0
521 ! if (ldiag2 > 0) then
522 ! if ((i == icol_diag) .and. (mod(k-1,5) == 0)) then
523 ! idiagss = 1
524 ! end if
525 ! end if
527 ! call modal_aero_soaexch( dtfull, temp, pres, &
528 ! niter, niter_max, ntot_soamode, &
529 ! g_soa_in, a_soa_in, a_poa_in, xferrate, &
530 ! g_soa_tend, a_soa_tend )
531 call modal_aero_soaexch( deltat, t(i,k), pmid(i,k), &
532 niter, niter_max, ntot_soamode, &
533 q(i,k,l_soag), qold_soa, qold_poa, uptkrate_soa, &
534 tmp1, dqdt_soa )
535 ! tmp1, dqdt_soa, g0_soa, idiagss )
536 sum_dqdt_soa = -tmp1
538 else if ( do_soag ) then
539 ! compute TMR tendencies for soa interstial aerosol
540 ! due to simple gas uptake
541 do n = 1, ntot_amode
542 dqdt_soa(n) = fgain_soa(n)*sum_dqdt_soa
543 end do
544 else
545 dqdt_soa(:) = 0.0
546 end if
548 do n = 1, ntot_amode
549 if (ido_so4a(n) == 1) then
550 l = lptr_so4_a_amode(n)-loffset
551 dqdt(i,k,l) = dqdt_so4(n)
552 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_so4(n)*pdel_fac
553 end if
555 if ( do_nh4g ) then
556 if (ido_nh4a(n) == 1) then
557 l = lptr_nh4_a_amode(n)-loffset
558 dqdt(i,k,l) = dqdt_nh4(n)
559 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_nh4(n)*pdel_fac
560 end if
561 end if
563 if ( do_soag ) then
564 if (ido_soaa(n) == 1) then
565 l = lptr_soa_a_amode(n)-loffset
566 dqdt(i,k,l) = dqdt_soa(n)
567 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_soa(n)*pdel_fac
568 end if
569 end if
570 end do
572 ! compute TMR tendencies for h2so4, nh3, and msa gas
573 ! due to simple gas uptake
574 l = l_so4g
575 dqdt(i,k,l) = -sum_dqdt_so4
576 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
578 if ( do_msag ) then
579 l = l_msag
580 dqdt(i,k,l) = -sum_dqdt_msa
581 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
582 end if
584 if ( do_nh4g ) then
585 l = l_nh4g
586 dqdt(i,k,l) = -sum_dqdt_nh4_b
587 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
588 end if
590 if ( do_soag ) then
591 l = l_soag
592 dqdt(i,k,l) = -sum_dqdt_soa
593 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
594 end if
596 ! compute TMR tendencies associated with primary carbon aging
597 if (modefrm_pcage > 0) then
598 n = modeptr_pcarbon
599 vol_shell = deltat * &
600 ( dqdt_so4(n)*fac_m2v_so4 + dqdt_nh4(n)*fac_m2v_nh4 + &
601 dqdt_soa(n)*fac_m2v_soa*soa_equivso4_factor )
602 vol_core = 0.0
603 do l = 1, nspec_amode(n)
604 vol_core = vol_core + &
605 q(i,k,lmassptr_amode(l,n)-loffset)*fac_m2v_pcarbon(l)
606 end do
607 ! ratio1 = vol_shell/vol_core =
608 ! actual hygroscopic-shell-volume/carbon-core-volume after gas uptake
609 ! ratio2 = 6.0_r8*dr_so4_monolayers_pcage/(dgncur_a*fac_volsfc_pcarbon)
610 ! = (shell-volume corresponding to n_so4_monolayers_pcage)/core-volume
611 ! The 6.0/(dgncur_a*fac_volsfc_pcarbon) = (mode-surface-area/mode-volume)
612 ! Note that vol_shell includes both so4+nh4 AND soa as "equivalent so4",
613 ! The soa_equivso4_factor accounts for the lower hygroscopicity of soa.
614 !
615 ! Define xferfrac_pcage = min( 1.0, ratio1/ratio2)
616 ! But ratio1/ratio2 == tmp1/tmp2, and coding below avoids possible overflow
617 !
618 tmp1 = vol_shell*dgncur_a(i,k,n)*fac_volsfc_pcarbon
619 tmp2 = max( 6.0_r8*dr_so4_monolayers_pcage*vol_core, 0.0_r8 )
620 if (tmp1 >= tmp2) then
621 xferfrac_pcage = xferfrac_max
622 else
623 xferfrac_pcage = min( tmp1/tmp2, xferfrac_max )
624 end if
626 if (xferfrac_pcage > 0.0_r8) then
627 do iq = 1, nspecfrm_pcage
628 lsfrm = lspecfrm_pcage(iq)-loffset
629 lstoo = lspectoo_pcage(iq)-loffset
630 xferrate = (xferfrac_pcage/deltat)*q(i,k,lsfrm)
631 dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferrate
632 qsrflx(i,lsfrm,jsrf) = qsrflx(i,lsfrm,jsrf) - xferrate*pdel_fac
633 if ((lstoo > 0) .and. (lstoo <= pcnst)) then
634 dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferrate
635 qsrflx(i,lstoo,jsrf) = qsrflx(i,lstoo,jsrf) + xferrate*pdel_fac
636 end if
637 end do
639 if (ido_so4a(modetoo_pcage) > 0) then
640 l = lptr_so4_a_amode(modetoo_pcage)-loffset
641 dqdt(i,k,l) = dqdt(i,k,l) + dqdt_so4(modefrm_pcage)
642 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_so4(modefrm_pcage)*pdel_fac
643 end if
645 if (ido_nh4a(modetoo_pcage) > 0) then
646 l = lptr_nh4_a_amode(modetoo_pcage)-loffset
647 dqdt(i,k,l) = dqdt(i,k,l) + dqdt_nh4(modefrm_pcage)
648 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_nh4(modefrm_pcage)*pdel_fac
649 end if
651 if (ido_soaa(modetoo_pcage) > 0) then
652 l = lptr_soa_a_amode(modetoo_pcage)-loffset
653 dqdt(i,k,l) = dqdt(i,k,l) + dqdt_soa(modefrm_pcage)
654 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_soa(modefrm_pcage)*pdel_fac
655 end if
656 end if
658 end if
661 ! diagnostics start -------------------------------------------------------
662 if (ldiag2 > 0) then
663 if (i == icol_diag) then
664 if (mod(k-1,5) == 0) then
665 write(*,'(a,43i5)') 'gasaerexch aaa nstep,lat,lon,k', nstep, latndx(i), lonndx(i), k
666 write(*,'(a,1p,10e12.4)') 'uptkratebb ', uptkratebb(:)
667 write(*,'(a,1p,10e12.4)') 'sum_uprt_so4 ', sum_uprt_so4
668 write(*,'(a,1p,10e12.4)') 'fgain_so4 ', fgain_so4(:)
669 write(*,'(a,1p,10e12.4)') 'sum_uprt_nh4 ', sum_uprt_nh4
670 write(*,'(a,1p,10e12.4)') 'fgain_nh4 ', fgain_nh4(:)
671 write(*,'(a,1p,10e12.4)') 'sum_uprt_soa ', sum_uprt_soa
672 write(*,'(a,1p,10e12.4)') 'fgain_soa ', fgain_soa(:)
673 write(*,'(a,1p,10e12.4)') 'so4g o,dqdt,n', q(i,k,l_so4g), sum_dqdt_so4, &
674 (q(i,k,l_so4g)-deltat*sum_dqdt_so4)
675 write(*,'(a,1p,10e12.4)') 'nh3g o,dqdt,n', q(i,k,l_nh4g), sum_dqdt_nh4, sum_dqdt_nh4_b, &
676 (q(i,k,l_nh4g)-deltat*sum_dqdt_nh4_b)
677 write(*,'(a,1p,10e12.4)') 'soag o,dqdt,n', q(i,k,l_soag), sum_dqdt_soa, &
678 (q(i,k,l_soag)-deltat*sum_dqdt_soa)
679 write(*,'(a,i12,1p,10e12.4)') &
680 'method,g0,t,p', method_soa, g0_soa, t(i,k), pmid(i,k)
681 write(*,'(a,1p,10e12.4)') 'so4 old ', qold_so4(:)
682 write(*,'(a,1p,10e12.4)') 'so4 dqdt ', dqdt_so4(:)
683 write(*,'(a,1p,10e12.4)') 'so4 new ', (qold_so4(:)+deltat*dqdt_so4(:))
684 write(*,'(a,1p,10e12.4)') 'nh4 old ', qold_nh4(:)
685 write(*,'(a,1p,10e12.4)') 'nh4 dqdt ', dqdt_nh4(:)
686 write(*,'(a,1p,10e12.4)') 'nh4 new ', (qold_nh4(:)+deltat*dqdt_nh4(:))
687 write(*,'(a,1p,10e12.4)') 'soa old ', qold_soa(:)
688 write(*,'(a,1p,10e12.4)') 'soa dqdt ', dqdt_soa(:)
689 write(*,'(a,1p,10e12.4)') 'soa new ', (qold_soa(:)+deltat*dqdt_soa(:))
690 write(*,'(a,1p,10e12.4)') 'vshell, core ', vol_shell, vol_core
691 write(*,'(a,1p,10e12.4)') 'dr_mono, ... ', dr_so4_monolayers_pcage, &
692 soa_equivso4_factor
693 write(*,'(a,1p,10e12.4)') 'dgn, ... ', dgncur_a(i,k,modefrm_pcage), &
694 fac_volsfc_pcarbon
695 write(*,'(a,1p,10e12.4)') 'tmp1, tmp2 ', tmp1, tmp2
696 write(*,'(a,1p,10e12.4)') 'xferfrac_age ', xferfrac_pcage
697 end if
698 end if
699 end if
700 ! diagnostics end ---------------------------------------------------------
703 end do ! "i = 1, ncol"
704 end do ! "k = 1, pver"
707 ! set "temporary testing arrays"
708 qold(:,:,:) = q(:,:,:)
709 qqcwold(:,:,:) = qqcw(:,:,:)
710 dqdtsv1(:,:,:) = dqdt(:,:,:)
711 dqqcwdtsv1(:,:,:) = dqqcwdt(:,:,:)
714 !
715 ! do renaming calcs
716 !
717 dotendrn(:) = .false.
718 dotendqqcwrn(:) = .false.
719 dorename_atik(1:ncol,:) = .true.
720 is_dorename_atik = .true.
721 if (ncol >= -13579) then
722 call modal_aero_rename_sub( &
723 'modal_aero_gasaerexch_sub', &
724 lchnk, ncol, nstep, &
725 loffset, deltat, &
726 latndx, lonndx, &
727 pdel, &
728 dotendrn, q, &
729 dqdt, dqdt_other, &
730 dotendqqcwrn, qqcw, &
731 dqqcwdt, dqqcwdt_other, &
732 is_dorename_atik, dorename_atik, &
733 jsrflx_rename, nsrflx, &
734 qsrflx, qqcwsrflx )
735 end if
738 !
739 ! apply the dqdt to update q (and same for qqcw)
740 !
741 do l = 1, pcnstxx
742 if ( dotend(l) .or. dotendrn(l) ) then
743 do k = 1, pver
744 do i = 1, ncol
745 q(i,k,l) = q(i,k,l) + dqdt(i,k,l)*deltat
746 end do
747 end do
748 end if
749 if ( dotendqqcw(l) .or. dotendqqcwrn(l) ) then
750 do k = 1, pver
751 do i = 1, ncol
752 qqcw(i,k,l) = qqcw(i,k,l) + dqqcwdt(i,k,l)*deltat
753 end do
754 end do
755 end if
756 end do
759 ! diagnostics start -------------------------------------------------------
760 if (ldiag3 > 0) then
761 if (icol_diag > 0) then
762 i = icol_diag
763 write(*,'(a,3i5)') 'gasaerexch ppp nstep,lat,lon', nstep, latndx(i), lonndx(i)
764 write(*,'(2i5,3(2x,a))') 0, 0, 'ppp', 'pdel for all k'
765 write(*,'(1p,7e12.4)') (pdel(i,k), k=1,pver)
767 write(*,'(a,3i5)') 'gasaerexch ddd nstep,lat,lon', nstep, latndx(i), lonndx(i)
768 do l = 1, pcnstxx
769 lb = l + loffset
771 if ( dotend(l) .or. dotendrn(l) ) then
772 write(*,'(2i5,3(2x,a))') 1, l, 'ddd1', cnst_name(lb), 'qold for all k'
773 write(*,'(1p,7e12.4)') (qold(i,k,l), k=1,pver)
774 write(*,'(2i5,3(2x,a))') 1, l, 'ddd2', cnst_name(lb), 'qnew for all k'
775 write(*,'(1p,7e12.4)') (q(i,k,l), k=1,pver)
776 write(*,'(2i5,3(2x,a))') 1, l, 'ddd3', cnst_name(lb), 'dqdt from conden for all k'
777 write(*,'(1p,7e12.4)') (dqdtsv1(i,k,l), k=1,pver)
778 write(*,'(2i5,3(2x,a))') 1, l, 'ddd4', cnst_name(lb), 'dqdt from rename for all k'
779 write(*,'(1p,7e12.4)') ((dqdt(i,k,l)-dqdtsv1(i,k,l)), k=1,pver)
780 write(*,'(2i5,3(2x,a))') 1, l, 'ddd5', cnst_name(lb), 'dqdt other for all k'
781 write(*,'(1p,7e12.4)') (dqdt_other(i,k,l), k=1,pver)
782 end if
784 if ( dotendqqcw(l) .or. dotendqqcwrn(l) ) then
785 write(*,'(2i5,3(2x,a))') 2, l, 'ddd1', cnst_name_cw(lb), 'qold for all k'
786 write(*,'(1p,7e12.4)') (qqcwold(i,k,l), k=1,pver)
787 write(*,'(2i5,3(2x,a))') 2, l, 'ddd2', cnst_name_cw(lb), 'qnew for all k'
788 write(*,'(1p,7e12.4)') (qqcw(i,k,l), k=1,pver)
789 write(*,'(2i5,3(2x,a))') 2, l, 'ddd3', cnst_name_cw(lb), 'dqdt from conden for all k'
790 write(*,'(1p,7e12.4)') (dqqcwdtsv1(i,k,l), k=1,pver)
791 write(*,'(2i5,3(2x,a))') 2, l, 'ddd4', cnst_name_cw(lb), 'dqdt from rename for all k'
792 write(*,'(1p,7e12.4)') ((dqqcwdt(i,k,l)-dqqcwdtsv1(i,k,l)), k=1,pver)
793 write(*,'(2i5,3(2x,a))') 2, l, 'ddd5', cnst_name_cw(lb), 'dqdt other for all k'
794 write(*,'(1p,7e12.4)') (dqqcwdt_other(i,k,l), k=1,pver)
795 end if
797 end do
799 write(*,'(a,3i5)') 'gasaerexch fff nstep,lat,lon', nstep, latndx(i), lonndx(i)
800 do l = 1, pcnstxx
801 lb = l + loffset
802 if ( dotend(l) .or. dotendrn(l) .or. dotendqqcw(l) .or. dotendqqcwrn(l) ) then
803 write(*,'(i5,2(2x,a,2l3))') l, &
804 cnst_name(lb), dotend(l), dotendrn(l), &
805 cnst_name_cw(lb), dotendqqcw(l), dotendqqcwrn(l)
806 end if
807 end do
809 end if
810 end if
811 ! diagnostics end ---------------------------------------------------------
814 ! do history file column-tendency fields
815 do l = 1, pcnstxx
816 lb = l + loffset
818 do jsrf = 1, 2
820 do jac = 1, 2
822 if (jac == 1) then
823 if (jsrf == jsrflx_gaexch) then
824 if ( .not. dotend(l) ) cycle
825 fieldname = trim(cnst_name(lb)) // '_sfgaex1'
826 else if (jsrf == jsrflx_rename) then
827 if ( .not. dotendrn(l) ) cycle
828 fieldname = trim(cnst_name(lb)) // '_sfgaex2'
829 else
830 cycle
831 end if
832 do i = 1, ncol
833 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf)*(adv_mass(l)/mwdry)
834 end do
835 call outfld( fieldname, qsrflx(:,l,jsrf), pcols, lchnk )
837 else
838 if (jsrf == jsrflx_gaexch) then
839 cycle
840 else if (jsrf == jsrflx_rename) then
841 if ( .not. dotendqqcwrn(l) ) cycle
842 fieldname = trim(cnst_name_cw(lb)) // '_sfgaex2'
843 else
844 cycle
845 end if
846 do i = 1, ncol
847 qqcwsrflx(i,l,jsrf) = qqcwsrflx(i,l,jsrf)*(adv_mass(l)/mwdry)
848 end do
849 call outfld( fieldname, qqcwsrflx(:,l,jsrf), pcols, lchnk )
850 end if
852 ! if (( masterproc ) .and. (nstep < 1)) &
853 ! write(*,'(2(a,2x),1p,e11.3)') &
854 ! 'modal_aero_newnuc_sub outfld', fieldname, adv_mass(l)
856 if (ldiag4 > 0) then
857 if (icol_diag > 0) then
858 i = icol_diag
859 if (jac == 1) then
860 tmp1 = qsrflx(i,l,jsrf)
861 else
862 tmp1 = qqcwsrflx(i,l,jsrf)
863 end if
864 write(*,'(a,4i5,2x,a,1p,2e12.4)') &
865 'gasaerexch nstep,lat,lon,l,fieldname,qsrflx,adv_mass', &
866 nstep, latndx(i), lonndx(i), l, fieldname, tmp1, adv_mass(l)
867 end if
868 end if
870 end do ! jac = ...
871 end do ! jsrf = ...
872 end do ! l = ...
875 return
876 !EOC
877 end subroutine modal_aero_gasaerexch_sub
880 !----------------------------------------------------------------------
881 !----------------------------------------------------------------------
882 subroutine gas_aer_uptkrates( ncol, loffset, &
883 q, t, pmid, &
884 dgncur_awet, uptkrate &
885 #ifdef WRF_PORT
886 , pcnstxx &
887 #endif
888 )
890 !
891 ! /
892 ! computes uptkrate = | dx dN/dx gas_conden_rate(Dp(x))
893 ! /
894 ! using Gauss-Hermite quadrature of order nghq=2
895 !
896 ! Dp = particle diameter (cm)
897 ! x = ln(Dp)
898 ! dN/dx = log-normal particle number density distribution
899 ! gas_conden_rate(Dp) = 2 * pi * gasdiffus * Dp * F(Kn,ac)
900 ! F(Kn,ac) = Fuchs-Sutugin correction factor
901 ! Kn = Knudsen number
902 ! ac = accomodation coefficient
903 !
905 !use modal_aero_data, only: ntot_amode, ntot_amode, nspec_amode, &
906 ! lspectype_amode, lmassptr_amode, &
907 ! sigmag_amode, &
908 ! specdens_amode, specmw_amode
909 use modal_aero_data, only: ntot_amode, ntot_amode, &
910 numptr_amode, &
911 sigmag_amode
912 #ifndef WRF_PORT
913 use ppgrid
914 use constituents, only: pcnst, cnst_name
915 #else
916 use module_cam_support, only: pcnst => pcnst_runtime, &
917 pcols, pver
918 use constituents, only: cnst_name
919 #endif
920 use physconst, only: mwdry, rair
922 implicit none
924 integer, intent(in) :: pcnstxx
925 integer, intent(in) :: ncol ! number of atmospheric column
926 integer, intent(in) :: loffset
927 real(r8), intent(in) :: q(ncol,pver,pcnstxx) ! Tracer array (mol,#/mol-air)
928 real(r8), intent(in) :: t(pcols,pver) ! Temperature in Kelvin
929 real(r8), intent(in) :: pmid(pcols,pver) ! Air pressure in Pa
930 real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)
932 real(r8), intent(out) :: uptkrate(ntot_amode,pcols,pver)
933 ! gas-to-aerosol mass transfer rates (1/s)
936 ! local
937 integer, parameter :: nghq = 2
938 integer :: i, iq, k, l1, l2, la, n
940 real(r8), parameter :: tworootpi = 3.5449077
941 real(r8), parameter :: root2 = 1.4142135
942 real(r8), parameter :: beta = 2.0
944 real(r8) :: aircon
945 real(r8) :: const
946 real(r8) :: dp, dum_m2v
947 real(r8) :: dryvol_a(pcols,pver)
948 real(r8) :: gasdiffus, gasspeed
949 real(r8) :: freepathx2, fuchs_sutugin
950 real(r8) :: knudsen
951 real(r8) :: lndp, lndpgn, lnsg
952 real(r8) :: num_a
953 real(r8) :: rhoair
954 real(r8) :: sumghq
955 real(r8), save :: xghq(nghq), wghq(nghq) ! quadrature abscissae and weights
957 data xghq / 0.70710678, -0.70710678 /
958 data wghq / 0.88622693, 0.88622693 /
961 ! outermost loop over all modes
962 do n = 1, ntot_amode
964 ! 22-aug-2007 rc easter - get number from q array rather
965 ! than computing a "bounded" number conc.
966 !! compute dry volume = sum_over_components{ component_mass / density }
967 !! (m3-AP/mol-air)
968 !! compute it for all i,k to improve accessing q array
969 ! dryvol_a(1:ncol,:) = 0.0_r8
970 ! do l1 = 1, nspec_amode(n)
971 ! l2 = lspectype_amode(l1,n)
972 !! dum_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
973 !! [m3-AP/kmol-AP]= [kg-AP/kmol-AP] / [kg-AP/m3-AP]
974 ! dum_m2v = specmw_amode(l2) / specdens_amode(l2)
975 ! la = lmassptr_amode(l1,n)
976 ! dryvol_a(1:ncol,:) = dryvol_a(1:ncol,:) &
977 ! + max(0.0_r8,q(1:ncol,:,la))*dum_m2v
978 ! end do
980 ! loops k and i
981 do k=1,pver
982 do i=1,ncol
984 rhoair = pmid(i,k)/(rair*t(i,k)) ! (kg-air/m3)
985 ! aircon = 1.0e3*rhoair/mwdry ! (mol-air/m3)
987 !! "bounded" number conc. (#/m3)
988 ! num_a = dryvol_a(i,k)*v2ncur_a(i,k,n)*aircon
990 ! number conc. (#/m3) -- note q(i,k,numptr) is (#/kmol-air)
991 ! so need aircon in (kmol-air/m3)
992 aircon = rhoair/mwdry ! (kmol-air/m3)
993 num_a = q(i,k,numptr_amode(n)-loffset)*aircon
995 ! gasdiffus = h2so4 gas diffusivity from mosaic code (m^2/s)
996 ! (pmid must be Pa)
997 gasdiffus = 0.557e-4 * (t(i,k)**1.75) / pmid(i,k)
998 ! gasspeed = h2so4 gas mean molecular speed from mosaic code (m/s)
999 gasspeed = 1.470e1 * sqrt(t(i,k))
1000 ! freepathx2 = 2 * (h2so4 mean free path) (m)
1001 freepathx2 = 6.0*gasdiffus/gasspeed
1003 lnsg = log( sigmag_amode(n) )
1004 lndpgn = log( dgncur_awet(i,k,n) ) ! (m)
1005 const = tworootpi * num_a * exp(beta*lndpgn + 0.5*(beta*lnsg)**2)
1007 ! sum over gauss-hermite quadrature points
1008 sumghq = 0.0
1009 do iq = 1, nghq
1010 lndp = lndpgn + beta*lnsg**2 + root2*lnsg*xghq(iq)
1011 dp = exp(lndp)
1013 ! knudsen number
1014 knudsen = freepathx2/dp
1015 ! following assumes accomodation coefficient = ac = 0.65
1016 ! (Adams & Seinfeld, 2002, JGR, and references therein)
1017 ! fuchs_sutugin = (0.75*ac*(1. + knudsen)) /
1018 ! (knudsen*(1.0 + knudsen + 0.283*ac) + 0.75*ac)
1019 fuchs_sutugin = (0.4875*(1. + knudsen)) / &
1020 (knudsen*(1.184 + knudsen) + 0.4875)
1022 sumghq = sumghq + wghq(iq)*dp*fuchs_sutugin/(dp**beta)
1023 end do
1024 uptkrate(n,i,k) = const * gasdiffus * sumghq
1026 end do ! "do i = 1, ncol"
1027 end do ! "do k = 1, pver"
1029 end do ! "do n = 1, ntot_soamode"
1032 return
1033 end subroutine gas_aer_uptkrates
1035 !----------------------------------------------------------------------
1038 subroutine modal_aero_soaexch( dtfull, temp, pres, &
1039 niter, niter_max, ntot_soamode, &
1040 g_soa_in, a_soa_in, a_poa_in, xferrate, &
1041 g_soa_tend, a_soa_tend )
1042 ! g_soa_tend, a_soa_tend, g0_soa, idiagss )
1044 !-----------------------------------------------------------------------
1045 !
1046 ! Purpose:
1047 !
1048 ! calculates condensation/evaporation of "soa gas"
1049 ! to/from multiple aerosol modes in 1 grid cell
1050 !
1051 ! key assumptions
1052 ! (1) ambient equilibrium vapor pressure of soa gas
1053 ! is given by p0_soa_298 and delh_vap_soa
1054 ! (2) equilibrium vapor pressure of soa gas at aerosol
1055 ! particle surface is given by raoults law in the form
1056 ! g_star = g0_soa*[a_soa/(a_soa + a_opoa)]
1057 ! (3) (oxidized poa)/(total poa) is equal to frac_opoa (constant)
1058 !
1059 !
1060 ! Author: R. Easter and R. Zaveri
1061 !
1062 !-----------------------------------------------------------------------
1063 implicit none
1065 real(r8), intent(in) :: dtfull ! full integration time step (s)
1066 real(r8), intent(in) :: temp ! air temperature (K)
1067 real(r8), intent(in) :: pres ! air pressure (Pa)
1068 integer, intent(out) :: niter ! number of iterations performed
1069 integer, intent(in) :: niter_max ! max allowed number of iterations
1070 integer, intent(in) :: ntot_soamode ! number of modes having soa
1071 real(r8), intent(in) :: g_soa_in ! initial soa gas mixrat (mol/mol)
1072 real(r8), intent(in) :: a_soa_in(ntot_soamode) ! initial soa aerosol mixrat (mol/mol)
1073 real(r8), intent(in) :: a_poa_in(ntot_soamode) ! initial poa aerosol mixrat (mol/mol)
1074 real(r8), intent(in) :: xferrate(ntot_soamode) ! gas-aerosol mass transfer rate (1/s)
1075 real(r8), intent(out) :: g_soa_tend ! soa gas mixrat tendency (mol/mol/s)
1076 real(r8), intent(out) :: a_soa_tend(ntot_soamode) ! soa aerosol mixrat tendency (mol/mol/s)
1077 ! real(r8), intent(out) :: g0_soa ! ambient soa gas equilib mixrat (mol/mol)
1078 ! integer, intent(in) :: idiagss
1080 integer :: luna=6
1081 integer :: m
1083 real(r8), parameter :: alpha = 0.05 ! parameter used in calc of time step
1084 real(r8), parameter :: g_min1 = 1.0e-20
1085 real(r8), parameter :: opoa_frac = 0.1 ! fraction of poa that is opoa
1086 real(r8), parameter :: delh_vap_soa = 156.0e3
1087 ! delh_vap_soa = heat of vaporization for gas soa (J/mol)
1088 real(r8), parameter :: p0_soa_298 = 1.0e-10
1089 ! p0_soa_298 = soa gas equilib vapor presssure (atm) at 298 k
1090 real(r8), parameter :: rgas = 8.3144 ! gas constant in J/K/mol
1092 real(r8) :: a_opoa(ntot_soamode) ! oxidized-poa aerosol mixrat (mol/mol)
1093 real(r8) :: a_soa(ntot_soamode) ! soa aerosol mixrat (mol/mol)
1094 real(r8) :: a_soa_tmp ! temporary soa aerosol mixrat (mol/mol)
1095 real(r8) :: beta(ntot_soamode) ! dtcur*xferrate
1096 real(r8) :: dtcur ! current time step (s)
1097 real(r8) :: dtmax ! = (dtfull-tcur)
1098 real(r8) :: g_soa ! soa gas mixrat (mol/mol)
1099 real(r8) :: g0_soa ! ambient soa gas equilib mixrat (mol/mol)
1100 real(r8) :: g_star(ntot_soamode) ! soa gas mixrat that is in equilib
1101 ! with each aerosol mode (mol/mol)
1102 real(r8) :: phi(ntot_soamode) ! "relative driving force"
1103 real(r8) :: p0_soa ! soa gas equilib vapor presssure (atm)
1104 real(r8) :: sat(ntot_soamode)
1105 real(r8) :: tcur ! current integration time (from 0 s)
1106 real(r8) :: tmpa, tmpb
1107 real(r8) :: tot_soa ! g_soa + sum( a_soa(:) )
1110 ! force things to be non-negative and calc tot_soa
1111 ! calc a_opoa (always slightly >0)
1112 g_soa = max( g_soa_in, 0.0_r8 )
1113 tot_soa = g_soa
1114 do m = 1, ntot_soamode
1115 a_soa(m) = max( a_soa_in(m), 0.0_r8 )
1116 tot_soa = tot_soa + a_soa(m)
1117 a_opoa(m) = opoa_frac*a_poa_in(m)
1118 a_opoa(m) = max( a_opoa(m), 1.0e-20_r8 ) ! force to small non-zero value
1119 end do
1121 ! calc ambient equilibrium soa gas
1122 p0_soa = p0_soa_298 * &
1123 exp( -(delh_vap_soa/rgas)*((1.0/temp)-(1.0/298.0)) )
1124 g0_soa = 1.01325e5*p0_soa/pres
1125 ! molecular weight adjustment
1126 ! the soa parameterization assumes that real gsoa and asoa have mw=150
1127 ! currently in cam3,
1128 ! mw=12 is used (this has to do with the mozart preprocessor)
1129 ! soag emission files (molec/cm2/s units) are set to give the desired
1130 ! mass emissions (kg/m2/s) and mass mixing ratios (kg/kg)
1131 ! when mw=12 is applied
1132 ! as a result, the molar mixing ratios for both gsoa and asoa
1133 ! are artificially scaled up by (150/12)
1134 ! and g0_soa must be similarly scaled up
1135 g0_soa = g0_soa*(150.0/12.0)
1136 ! g0_soa = 0.0 ! force irreversible uptake
1138 niter = 0
1139 tcur = 0.0
1140 dtcur = 0.0
1141 phi(:) = 0.0
1142 g_star(:) = 0.0
1144 ! if (idiagss > 0) then
1145 ! write(luna,'(a,1p,10e11.3)') 'p0, g0_soa', p0_soa, g0_soa
1146 ! write(luna,'(3a)') &
1147 ! 'niter, tcur, dtcur, phi(:), ', &
1148 ! 'g_star(:), ', &
1149 ! 'a_soa(:), g_soa'
1150 ! write(luna,'(3a)') &
1151 ! ' sat(:), ', &
1152 ! 'sat(:)*a_soa(:) ', &
1153 ! 'a_opoa(:)'
1154 ! write(luna,'(i3,1p,20e10.2)') niter, tcur, dtcur, &
1155 ! phi(:), g_star(:), a_soa(:), g_soa
1156 ! end if
1159 ! integration loop -- does multiple substeps to reach dtfull
1160 timeloop: do while (tcur < dtfull-1.0e-3_r8 )
1162 niter = niter + 1
1163 if (niter > niter_max) exit
1165 tmpa = 0.0
1166 do m = 1, ntot_soamode
1167 sat(m) = g0_soa/(a_soa(m) + a_opoa(m))
1168 g_star(m) = sat(m)*a_soa(m)
1169 phi(m) = (g_soa - g_star(m))/max(g_soa,g_star(m),g_min1)
1170 tmpa = tmpa + xferrate(m)*abs(phi(m))
1171 end do
1173 dtmax = dtfull-tcur
1174 if (dtmax*tmpa <= alpha) then
1175 ! here alpha/tmpa >= dtmax, so this is final substep
1176 dtcur = dtmax
1177 tcur = dtfull
1178 else
1179 dtcur = alpha/tmpa
1180 tcur = tcur + dtcur
1181 end if
1183 ! step 1 - for modes where soa is condensing, estimate "new" a_soa(m)
1184 ! using an explicit calculation with "old" g_soa
1185 ! and g_star(m) calculated using "old" a_soa(m)
1186 ! do this to get better estimate of "new" a_soa(m) and sat(m)
1187 do m = 1, ntot_soamode
1188 beta(m) = dtcur*xferrate(m)
1189 tmpa = g_soa - g_star(m)
1190 if (tmpa > 0.0_r8) then
1191 a_soa_tmp = a_soa(m) + beta(m)*tmpa
1192 sat(m) = g0_soa/(a_soa_tmp + a_opoa(m))
1193 g_star(m) = sat(m)*a_soa_tmp ! this just needed for diagnostics
1194 end if
1195 end do
1197 ! step 2 - implicit in g_soa and semi-implicit in a_soa,
1198 ! with g_star(m) calculated semi-implicitly
1199 tmpa = 0.0
1200 tmpb = 0.0
1201 do m = 1, ntot_soamode
1202 tmpa = tmpa + a_soa(m)/(1.0_r8 + beta(m)*sat(m))
1203 tmpb = tmpb + beta(m)/(1.0_r8 + beta(m)*sat(m))
1204 end do
1206 g_soa = (tot_soa - tmpa)/(1.0_r8 + tmpb)
1207 g_soa = max( 0.0_r8, g_soa )
1208 do m = 1, ntot_soamode
1209 a_soa(m) = (a_soa(m) + beta(m)*g_soa)/ &
1210 (1.0_r8 + beta(m)*sat(m))
1211 end do
1213 ! if (idiagss > 0) then
1214 ! write(luna,'(i3,1p,20e10.2)') niter, tcur, dtcur, &
1215 ! phi(:), g_star(:), a_soa(:), g_soa
1216 ! write(luna,'(23x,1p,20e10.2)') &
1217 ! sat(:), sat(:)*a_soa(:), a_opoa(:)
1218 ! end if
1220 ! if (niter > 9992000) then
1221 ! write(luna,'(a)') '*** to many iterations'
1222 ! exit
1223 ! end if
1225 end do timeloop
1228 g_soa_tend = (g_soa - g_soa_in)/dtfull
1229 do m = 1, ntot_soamode
1230 a_soa_tend(m) = (a_soa(m) - a_soa_in(m))/dtfull
1231 end do
1234 return
1235 end subroutine modal_aero_soaexch
1237 !----------------------------------------------------------------------
1240 subroutine modal_aero_gasaerexch_init
1242 !-----------------------------------------------------------------------
1243 !
1244 ! Purpose:
1245 ! set do_adjust and do_aitken flags
1246 ! create history fields for column tendencies associated with
1247 ! modal_aero_calcsize
1248 !
1249 ! Author: R. Easter
1250 !
1251 !-----------------------------------------------------------------------
1253 use modal_aero_data
1254 use modal_aero_rename
1255 #ifndef WRF_PORT
1256 use abortutils, only : endrun
1257 use cam_history, only : addfld, add_default, fieldname_len, phys_decomp
1258 use constituents, only : pcnst, cnst_get_ind, cnst_name
1259 use spmd_utils, only : masterproc
1260 use phys_control,only : phys_getopts
1261 #else
1262 use module_cam_support, only: pcnst => pcnst_runtime, &
1263 fieldname_len, &
1264 endrun, masterproc, phys_decomp, &
1265 add_default, addfld
1266 use constituents, only : cnst_name, cnst_get_ind
1267 #endif
1270 implicit none
1272 !-----------------------------------------------------------------------
1273 ! arguments
1275 !-----------------------------------------------------------------------
1276 ! local
1277 integer :: ipair, iq, iqfrm, iqfrm_aa, iqtoo, iqtoo_aa
1278 integer :: jac
1279 integer :: l, lsfrm, lstoo, lunout
1280 integer :: l_so4g, l_nh4g, l_msag, l_soag
1281 integer :: m, mfrm, mtoo
1282 integer :: nsamefrm, nsametoo, nspec
1283 integer :: n, nacc, nait
1285 logical :: do_msag, do_nh4g, do_soag
1286 logical :: dotend(pcnst), dotendqqcw(pcnst)
1288 real(r8) :: tmp1, tmp2
1290 character(len=fieldname_len) :: tmpnamea, tmpnameb
1291 character(len=fieldname_len+3) :: fieldname
1292 character(128) :: long_name
1293 character(8) :: unit
1295 logical :: history_aerosol ! Output the MAM aerosol tendencies
1297 !-----------------------------------------------------------------------
1298 #ifndef WRF_PORT
1299 call phys_getopts( history_aerosol_out = history_aerosol )
1300 #else
1301 history_aerosol = .FALSE.
1302 #endif
1304 lunout = 6
1305 !
1306 ! define "from mode" and "to mode" for primary carbon aging
1307 !
1308 ! skip (turn off) aging if either is absent,
1309 ! or if accum mode so4 is absent
1310 !
1311 modefrm_pcage = -999888777
1312 modetoo_pcage = -999888777
1313 if ((modeptr_pcarbon <= 0) .or. (modeptr_accum <= 0)) goto 15000
1314 l = lptr_so4_a_amode(modeptr_accum)
1315 if ((l < 1) .or. (l > pcnst)) goto 15000
1317 modefrm_pcage = modeptr_pcarbon
1318 modetoo_pcage = modeptr_accum
1320 !
1321 ! define species involved in each primary carbon aging pairing
1322 ! (include aerosol water)
1323 !
1324 !
1325 mfrm = modefrm_pcage
1326 mtoo = modetoo_pcage
1328 nspec = 0
1329 aa_iqfrm: do iqfrm = -1, nspec_amode(mfrm)
1331 if (iqfrm == -1) then
1332 lsfrm = numptr_amode(mfrm)
1333 lstoo = numptr_amode(mtoo)
1334 else if (iqfrm == 0) then
1335 ! bypass transfer of aerosol water due to primary-carbon aging
1336 cycle aa_iqfrm
1337 ! lsfrm = lwaterptr_amode(mfrm)
1338 ! lstoo = lwaterptr_amode(mtoo)
1339 else
1340 lsfrm = lmassptr_amode(iqfrm,mfrm)
1341 lstoo = 0
1342 end if
1343 if ((lsfrm < 1) .or. (lsfrm > pcnst)) cycle aa_iqfrm
1345 if (lsfrm>0 .and. iqfrm>0 ) then
1346 ! find "too" species having same lspectype_amode as the "frm" species
1347 do iqtoo = 1, nspec_amode(mtoo)
1348 if ( lspectype_amode(iqtoo,mtoo) .eq. &
1349 lspectype_amode(iqfrm,mfrm) ) then
1350 lstoo = lmassptr_amode(iqtoo,mtoo)
1351 exit
1352 end if
1353 end do
1354 end if
1356 if ((lstoo < 1) .or. (lstoo > pcnst)) lstoo = 0
1357 nspec = nspec + 1
1358 lspecfrm_pcage(nspec) = lsfrm
1359 lspectoo_pcage(nspec) = lstoo
1360 end do aa_iqfrm
1362 nspecfrm_pcage = nspec
1364 !
1365 ! output results
1366 !
1367 if ( masterproc ) then
1369 write(lunout,9310)
1371 mfrm = modefrm_pcage
1372 mtoo = modetoo_pcage
1373 write(lunout,9320) 1, mfrm, mtoo
1375 do iq = 1, nspecfrm_pcage
1376 lsfrm = lspecfrm_pcage(iq)
1377 lstoo = lspectoo_pcage(iq)
1378 if (lstoo .gt. 0) then
1379 write(lunout,9330) lsfrm, cnst_name(lsfrm), &
1380 lstoo, cnst_name(lstoo)
1381 else
1382 write(lunout,9340) lsfrm, cnst_name(lsfrm)
1383 end if
1384 end do
1386 write(lunout,*)
1388 end if ! ( masterproc )
1390 9310 format( / 'subr. modal_aero_gasaerexch_init - primary carbon aging pointers' )
1391 9320 format( 'pair', i3, 5x, 'mode', i3, ' ---> mode', i3 )
1392 9330 format( 5x, 'spec', i3, '=', a, ' ---> spec', i3, '=', a )
1393 9340 format( 5x, 'spec', i3, '=', a, ' ---> LOSS' )
1396 15000 continue
1398 ! set gas species indices
1399 call cnst_get_ind( 'H2SO4', l_so4g, .false. )
1400 call cnst_get_ind( 'NH3', l_nh4g, .false. )
1401 call cnst_get_ind( 'MSA', l_msag, .false. )
1402 call cnst_get_ind( 'SOAG', l_soag, .false. )
1403 if ((l_so4g <= 0) .or. (l_so4g > pcnst)) then
1404 write( *, '(/a/a,2i7)' ) &
1405 '*** modal_aero_gasaerexch_init -- cannot find H2SO4 species', &
1406 ' l_so4g=', l_so4g
1407 call endrun( 'modal_aero_gasaerexch_init error' )
1408 end if
1409 do_nh4g = .false.
1410 do_msag = .false.
1411 do_soag = .false.
1412 if ((l_nh4g > 0) .and. (l_nh4g <= pcnst)) do_nh4g = .true.
1413 if ((l_msag > 0) .and. (l_msag <= pcnst)) do_msag = .true.
1414 if ((l_soag > 0) .and. (l_soag <= pcnst)) do_soag = .true.
1416 ! set tendency flags
1417 dotend(:) = .false.
1418 dotend(l_so4g) = .true.
1419 if ( do_nh4g ) dotend(l_nh4g) = .true.
1420 if ( do_msag ) dotend(l_msag) = .true.
1421 if ( do_soag ) dotend(l_soag) = .true.
1422 do n = 1, ntot_amode
1423 l = lptr_so4_a_amode(n)
1424 if ((l > 0) .and. (l <= pcnst)) then
1425 dotend(l) = .true.
1426 if ( do_nh4g ) then
1427 l = lptr_nh4_a_amode(n)
1428 if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
1429 end if
1430 end if
1431 l = lptr_soa_a_amode(n)
1432 if ((l > 0) .and. (l <= pcnst)) then
1433 dotend(l) = .true.
1434 end if
1435 end do
1437 if (modefrm_pcage > 0) then
1438 do iq = 1, nspecfrm_pcage
1439 lsfrm = lspecfrm_pcage(iq)
1440 lstoo = lspectoo_pcage(iq)
1441 if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
1442 dotend(lsfrm) = .true.
1443 if ((lstoo > 0) .and. (lstoo <= pcnst)) then
1444 dotend(lstoo) = .true.
1445 end if
1446 end if
1447 end do
1448 end if
1451 ! define history fields for basic gas-aer exchange
1452 ! and primary carbon aging from that
1453 do l = 1, pcnst
1454 if ( .not. dotend(l) ) cycle
1456 tmpnamea = cnst_name(l)
1457 fieldname = trim(tmpnamea) // '_sfgaex1'
1458 long_name = trim(tmpnamea) // ' gas-aerosol-exchange primary column tendency'
1459 unit = 'kg/m2/s'
1460 call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
1461 if ( history_aerosol ) then
1462 call add_default( fieldname, 1, ' ' )
1463 endif
1464 if ( masterproc ) write(*,'(3(a,3x))') 'gasaerexch addfld', fieldname, unit
1466 end do ! l = ...
1469 ! define history fields for aitken-->accum renaming
1470 dotend(:) = .false.
1471 dotendqqcw(:) = .false.
1472 do ipair = 1, npair_renamexf
1473 do iq = 1, nspecfrm_renamexf(ipair)
1474 lsfrm = lspecfrma_renamexf(iq,ipair)
1475 lstoo = lspectooa_renamexf(iq,ipair)
1476 if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
1477 dotend(lsfrm) = .true.
1478 if ((lstoo > 0) .and. (lstoo <= pcnst)) then
1479 dotend(lstoo) = .true.
1480 end if
1481 end if
1483 lsfrm = lspecfrmc_renamexf(iq,ipair)
1484 lstoo = lspectooc_renamexf(iq,ipair)
1485 if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
1486 dotendqqcw(lsfrm) = .true.
1487 if ((lstoo > 0) .and. (lstoo <= pcnst)) then
1488 dotendqqcw(lstoo) = .true.
1489 end if
1490 end if
1491 end do ! iq = ...
1492 end do ! ipair = ...
1494 do l = 1, pcnst
1495 do jac = 1, 2
1496 if (jac == 1) then
1497 if ( .not. dotend(l) ) cycle
1498 tmpnamea = cnst_name(l)
1499 else
1500 if ( .not. dotendqqcw(l) ) cycle
1501 tmpnamea = cnst_name_cw(l)
1502 end if
1504 fieldname = trim(tmpnamea) // '_sfgaex2'
1505 long_name = trim(tmpnamea) // ' gas-aerosol-exchange renaming column tendency'
1506 unit = 'kg/m2/s'
1507 if ((tmpnamea(1:3) == 'num') .or. &
1508 (tmpnamea(1:3) == 'NUM')) unit = '#/m2/s'
1509 call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
1510 if ( history_aerosol ) then
1511 call add_default( fieldname, 1, ' ' )
1512 endif
1513 if ( masterproc ) write(*,'(3(a,3x))') 'gasaerexch addfld', fieldname, unit
1514 end do ! jac = ...
1515 end do ! l = ...
1518 ! calculate soa_equivso4_factor
1519 ! if do_soag == .false., then set it to zero as a safety measure
1520 soa_equivso4_factor = 0.0
1521 if ( do_soag ) then
1522 tmp1 = -1.0 ; tmp2 = -1.0
1523 do l = 1, ntot_aspectype
1524 if (specname_amode(l) == 's-organic') tmp1 = spechygro(l)
1525 if (specname_amode(l) == 'sulfate' ) tmp2 = spechygro(l)
1526 end do
1527 if ((tmp1 > 0.0_r8) .and. (tmp2 > 0.0_r8)) then
1528 soa_equivso4_factor = tmp1/tmp2
1529 else
1530 write(*,'(a/a,1p,2e10.2)') &
1531 '*** subr modal_aero_gasaerexch_init', &
1532 ' cannot find hygros - tmp1/2 =', tmp1, tmp2
1533 call endrun()
1534 end if
1535 end if
1537 return
1538 end subroutine modal_aero_gasaerexch_init
1541 !----------------------------------------------------------------------
1543 end module modal_aero_gasaerexch