Line data Source code
1 : ! modal_aero_gasaerexch.F90
2 :
3 :
4 : !----------------------------------------------------------------------
5 : !----------------------------------------------------------------------
6 : !BOP
7 : !
8 : ! !MODULE: modal_aero_gasaerexch --- does modal aerosol gas-aerosol exchange
9 : !
10 : ! !INTERFACE:
11 : module modal_aero_gasaerexch
12 :
13 : ! !USES:
14 : use shr_kind_mod, only: r8 => shr_kind_r8
15 : use chem_mods, only: gas_pcnst
16 : use modal_aero_data, only: nspec_max, nsoa, npoa, soa_multi_species
17 : use ref_pres, only: top_lev => clim_modal_aero_top_lev
18 : use ppgrid, only: pcols, pver
19 : use modal_aero_data, only: ntot_amode, numptr_amode, sigmag_amode
20 : use modal_aero_data, only: lptr2_soa_g_amode, lptr2_soa_a_amode, lptr2_pom_a_amode
21 :
22 : implicit none
23 : private
24 : save
25 :
26 : ! !PUBLIC MEMBER FUNCTIONS:
27 : public modal_aero_gasaerexch_sub, modal_aero_gasaerexch_init
28 :
29 : ! !PUBLIC DATA MEMBERS:
30 : integer, parameter :: pcnstxx = gas_pcnst
31 : integer, protected, public :: maxspec_pcage != nspec_max
32 :
33 : integer, protected, public :: modefrm_pcage
34 : integer, protected, public :: nspecfrm_pcage
35 : integer :: modetoo_pcage
36 :
37 : integer, protected, allocatable, public :: lspecfrm_pcage(:)
38 : integer, protected, allocatable, public :: lspectoo_pcage(:)
39 :
40 : real(r8), parameter, public :: n_so4_monolayers_pcage = 8.0_r8
41 :
42 : ! number of so4(+nh4) monolayers needed to "age" a carbon particle
43 :
44 : real(r8), parameter, public :: &
45 : dr_so4_monolayers_pcage = n_so4_monolayers_pcage * 4.76e-10_r8
46 : ! thickness of the so4 monolayers (m)
47 : ! for so4(+nh4), use bi-sulfate mw and 1.77 g/cm3,
48 : ! --> 1 mol so4(+nh4) = 65 cm^3 --> 1 molecule = (4.76e-10 m)^3
49 : ! aging criterion is approximate so do not try to distinguish
50 : ! sulfuric acid, bisulfate, ammonium sulfate
51 :
52 : real(r8), protected, allocatable, public :: soa_equivso4_factor(:)
53 : ! this factor converts an soa volume to a volume of so4(+nh4)
54 : ! having same hygroscopicity as the soa
55 :
56 : real (r8) :: fac_m2v_nh4, fac_m2v_so4
57 : real (r8), allocatable :: fac_m2v_soa(:)
58 :
59 : real (r8), allocatable :: fac_m2v_pcarbon(:)
60 :
61 : ! !DESCRIPTION: This module implements ...
62 : !
63 : ! !REVISION HISTORY:
64 : !
65 : ! RCE 07.04.13: Adapted from MIRAGE2 code
66 : !
67 : !EOP
68 : !----------------------------------------------------------------------
69 : !BOC
70 :
71 : ! list private module data here
72 :
73 : !EOC
74 : !----------------------------------------------------------------------
75 :
76 :
77 : contains
78 :
79 :
80 : !----------------------------------------------------------------------
81 : !----------------------------------------------------------------------
82 : !BOP
83 : ! !ROUTINE: modal_aero_gasaerexch_sub --- ...
84 : !
85 : ! !INTERFACE:
86 1490712 : subroutine modal_aero_gasaerexch_sub( &
87 : lchnk, ncol, nstep, &
88 : loffset, deltat, &
89 : t, pmid, pdel, &
90 : qh2o, troplev, &
91 1489176 : q, qqcw, &
92 1489176 : dqdt_other, dqqcwdt_other, &
93 1489176 : dgncur_a, dgncur_awet, &
94 : sulfeq )
95 :
96 : ! !USES:
97 : use modal_aero_data, only: alnsg_amode,lmassptr_amode,cnst_name_cw
98 : use modal_aero_data, only: lptr_so4_a_amode,lptr_nh4_a_amode
99 : use modal_aero_data, only: modeptr_pcarbon,nspec_amode,specmw_amode,specdens_amode
100 : use modal_aero_rename, only: modal_aero_rename_sub
101 : use rad_constituents, only: rad_cnst_get_info
102 : use constituents, only: pcnst, cnst_mw
103 :
104 : use cam_history, only: outfld, fieldname_len
105 : use chem_mods, only: adv_mass
106 : use constituents, only: pcnst, cnst_name, cnst_get_ind
107 : use mo_tracname, only: solsym
108 : use physconst, only: gravit, mwdry, rair
109 : use cam_abortutils, only: endrun
110 : use spmd_utils, only: iam, masterproc
111 : use phys_control, only: cam_chempkg_is
112 :
113 : implicit none
114 :
115 : ! !PARAMETERS:
116 : integer, intent(in) :: lchnk ! chunk identifier
117 : integer, intent(in) :: ncol ! number of atmospheric column
118 : integer, intent(in) :: nstep ! model time-step number
119 : integer, intent(in) :: loffset ! offset applied to modal aero "ptrs"
120 : integer, intent(in) :: troplev(pcols) ! tropopause vertical index
121 : real(r8), intent(in) :: deltat ! time step (s)
122 :
123 : real(r8), intent(inout) :: q(ncol,pver,pcnstxx) ! tracer mixing ratio (TMR) array
124 : ! *** MUST BE #/kmol-air for number
125 : ! *** MUST BE mol/mol-air for mass
126 : ! *** NOTE ncol dimension
127 : real(r8), intent(inout) :: qqcw(ncol,pver,pcnstxx)
128 : ! like q but for cloud-borner tracers
129 : real(r8), intent(in) :: dqdt_other(ncol,pver,pcnstxx)
130 : ! TMR tendency from other continuous
131 : ! growth processes (aqchem, soa??)
132 : ! *** NOTE ncol dimension
133 : real(r8), intent(in) :: dqqcwdt_other(ncol,pver,pcnstxx)
134 : ! like dqdt_other but for cloud-borner tracers
135 : real(r8), intent(in) :: t(pcols,pver) ! temperature at model levels (K)
136 : real(r8), intent(in) :: pmid(pcols,pver) ! pressure at model levels (Pa)
137 : real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickness of levels (Pa)
138 : real(r8), intent(in) :: qh2o(pcols,pver) ! water vapor mixing ratio (kg/kg)
139 : real(r8), intent(in) :: dgncur_a(pcols,pver,ntot_amode)
140 : real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)
141 : real(r8), pointer :: sulfeq(:,:,:)
142 :
143 : ! dry & wet geo. mean dia. (m) of number distrib.
144 :
145 : ! !DESCRIPTION:
146 : ! computes TMR (tracer mixing ratio) tendencies for gas condensation
147 : ! onto aerosol particles
148 : !
149 : ! this version does condensation of H2SO4, NH3, and MSA, both treated as
150 : ! completely non-volatile (gas --> aerosol, but no aerosol --> gas)
151 : ! gas H2SO4 goes to aerosol SO4
152 : ! gas MSA (if present) goes to aerosol SO4
153 : ! aerosol MSA is not distinguished from aerosol SO4
154 : ! gas NH3 (if present) goes to aerosol NH4
155 : ! if gas NH3 is not present, then ????
156 : !
157 : !
158 : ! !REVISION HISTORY:
159 : ! RCE 07.04.13: Adapted from MIRAGE2 code
160 : !
161 : !EOP
162 : !----------------------------------------------------------------------
163 : !BOC
164 :
165 : ! local variables
166 : integer, parameter :: jsrflx_gaexch = 1
167 : integer, parameter :: jsrflx_rename = 2
168 : integer, parameter :: ldiag1=-1, ldiag2=-1, ldiag3=-1, ldiag4=-1
169 : integer, parameter :: method_soa = 2
170 : ! method_soa=0 is no uptake
171 : ! method_soa=1 is irreversible uptake done like h2so4 uptake
172 : ! method_soa=2 is reversible uptake using subr modal_aero_soaexch
173 :
174 : integer :: i, iq, itmpa
175 : integer :: idiagss
176 2978352 : integer :: ido_so4a(ntot_amode), ido_nh4a(ntot_amode)
177 2978352 : integer :: ido_soaa(ntot_amode,nsoa)
178 : integer :: j, jac, jsrf, jsoa
179 : integer :: k,p
180 : integer :: l, l2, lb, lsfrm, lstoo
181 : integer :: l_so4g, l_nh4g, l_msag
182 2978352 : integer :: l_soag(nsoa)
183 : integer :: n, nn, niter, niter_max, ntot_soamode
184 :
185 2978352 : logical :: is_dorename_atik, dorename_atik(ncol,pver)
186 :
187 : character(len=fieldname_len+3) :: fieldname
188 : character(len=100) :: msg ! string for endrun calls
189 : character(len=32) :: spec_type
190 :
191 2978352 : real (r8) :: avg_uprt_nh4, avg_uprt_so4, avg_uprt_soa(nsoa)
192 : real (r8) :: deltatxx
193 2978352 : real (r8) :: dqdt_nh4(ntot_amode), dqdt_so4(ntot_amode)
194 2978352 : real (r8) :: dqdt_soa(ntot_amode,nsoa)
195 2978352 : real (r8) :: dqdt_soag(nsoa)
196 : real (r8) :: fac_volsfc_pcarbon
197 2978352 : real (r8) :: fgain_nh4(ntot_amode), fgain_so4(ntot_amode)
198 2978352 : real (r8) :: fgain_soa(ntot_amode,nsoa)
199 : real (r8) :: g0_soa(nsoa)
200 2978352 : real(r8) :: mw_poa_host(npoa) ! molec wght of poa used in host code
201 2978352 : real(r8) :: mw_soa_host(nsoa) ! molec wght of poa used in host code
202 : real (r8) :: pdel_fac
203 : real (r8) :: qmax_nh4, qnew_nh4, qnew_so4
204 2978352 : real (r8) :: qold_nh4(ntot_amode), qold_so4(ntot_amode)
205 2978352 : real (r8) :: qold_poa(ntot_amode,npoa)
206 2978352 : real (r8) :: qold_soa(ntot_amode,nsoa)
207 2978352 : real (r8) :: qold_soag(nsoa)
208 : real (r8) :: sum_dqdt_msa, sum_dqdt_so4
209 2978352 : real (r8) :: sum_dqdt_soa(nsoa)
210 : real (r8) :: sum_dqdt_nh4, sum_dqdt_nh4_b
211 : real (r8) :: sum_uprt_msa, sum_uprt_nh4, sum_uprt_so4
212 2978352 : real (r8) :: sum_uprt_soa(nsoa)
213 : real (r8) :: tmp1, tmp2, tmpa
214 : real (r8) :: tmp_kxt, tmp_pxt
215 : real (r8) :: tmp_so4a_bgn, tmp_so4a_end
216 : real (r8) :: tmp_so4g_avg, tmp_so4g_bgn, tmp_so4g_equ
217 2978352 : real (r8) :: uptkrate(ntot_amode,pcols,pver)
218 2978352 : real (r8) :: uptkratebb(ntot_amode)
219 2978352 : real (r8) :: uptkrate_soa(ntot_amode,nsoa)
220 : ! gas-to-aerosol mass transfer rates (1/s)
221 : real (r8) :: vol_core, vol_shell
222 : real (r8) :: xferfrac_pcage, xferfrac_max
223 : real (r8) :: xferrate
224 :
225 : logical :: do_msag ! true if msa gas is a species
226 : logical :: do_nh4g ! true if nh3 gas is a species
227 : logical :: do_soag_any ! true if soa gas is a species
228 2978352 : logical :: do_soag(nsoa) ! true if soa gas is a species
229 :
230 : logical :: dotend(pcnstxx) ! identifies species directly involved in
231 : ! gas-aerosol exchange (gas condensation)
232 : logical :: dotendqqcw(pcnstxx) ! like dotend but for cloud-borner tracers
233 : logical :: dotendrn(pcnstxx), dotendqqcwrn(pcnstxx)
234 : ! identifies species involved in renaming
235 : ! after "continuous growth"
236 : ! (gas-aerosol exchange and aqchem)
237 :
238 : integer, parameter :: nsrflx = 2 ! last dimension of qsrflx
239 2978352 : real(r8) :: dqdt(ncol,pver,pcnstxx) ! TMR "delta q" array - NOTE dims
240 2978352 : real(r8) :: dqqcwdt(ncol,pver,pcnstxx) ! like dqdt but for cloud-borner tracers
241 : real(r8) :: qsrflx(pcols,pcnstxx,nsrflx)
242 : ! process-specific column tracer tendencies
243 : ! (1=renaming, 2=gas condensation)
244 : real(r8) :: qconff(pcols,pver),qevapff(pcols,pver)
245 : real(r8) :: qconbb(pcols,pver),qevapbb(pcols,pver)
246 : real(r8) :: qconbg(pcols,pver),qevapbg(pcols,pver)
247 : real(r8) :: qcon(pcols,pver),qevap(pcols,pver)
248 :
249 : real(r8) :: qqcwsrflx(pcols,pcnstxx,nsrflx)
250 :
251 : ! following only needed for diagnostics
252 2978352 : real(r8) :: qold(ncol,pver,pcnstxx) ! NOTE dims
253 : real(r8) :: qnew(ncol,pver,pcnstxx) ! NOTE dims
254 : real(r8) :: qdel(ncol,pver,pcnstxx) ! NOTE dims
255 : real(r8) :: dumavec(1000), dumbvec(1000), dumcvec(1000)
256 2978352 : real(r8) :: qqcwold(ncol,pver,pcnstxx)
257 2978352 : real(r8) :: dqdtsv1(ncol,pver,pcnstxx)
258 2978352 : real(r8) :: dqqcwdtsv1(ncol,pver,pcnstxx)
259 :
260 :
261 : !----------------------------------------------------------------------
262 :
263 : ! set gas species indices
264 1489176 : call cnst_get_ind( 'H2SO4', l_so4g, .false. )
265 1489176 : call cnst_get_ind( 'NH3', l_nh4g, .false. )
266 1489176 : if ( .not. cam_chempkg_is('geoschem_mam4') ) then
267 1489176 : call cnst_get_ind( 'MSA', l_msag, .false. )
268 : else
269 0 : l_msag = 0
270 : endif
271 1489176 : l_so4g = l_so4g - loffset
272 1489176 : l_nh4g = l_nh4g - loffset
273 1489176 : l_msag = l_msag - loffset
274 1489176 : if ((l_so4g <= 0) .or. (l_so4g > pcnstxx)) then
275 : write( *, '(/a/a,2i7)' ) &
276 0 : '*** modal_aero_gasaerexch_sub -- cannot find H2SO4 species', &
277 0 : ' l_so4g, loffset =', l_so4g, loffset
278 0 : call endrun( 'modal_aero_gasaerexch_sub error' )
279 : end if
280 1489176 : do_nh4g = .false.
281 1489176 : do_msag = .false.
282 1489176 : if ((l_nh4g > 0) .and. (l_nh4g <= pcnstxx)) do_nh4g = .true.
283 1489176 : if ((l_msag > 0) .and. (l_msag <= pcnstxx)) do_msag = .true.
284 :
285 1489176 : do_soag_any = .false.
286 2978352 : do_soag(:) = .false.
287 2978352 : do jsoa = 1, nsoa
288 1489176 : l_soag(jsoa) = lptr2_soa_g_amode(jsoa) - loffset
289 1489176 : if ((method_soa == 1) .or. (method_soa == 2)) then
290 1489176 : if ((l_soag(jsoa) > 0) .and. (l_soag(jsoa) <= pcnstxx)) then
291 1489176 : do_soag_any = .true.
292 1489176 : do_soag(jsoa) = .true.
293 : end if
294 : else if (method_soa /= 0) then
295 : write(*,'(/a,1x,i10)') '*** modal_aero_gasaerexch_sub - bad method_soa =', method_soa
296 : call endrun( 'modal_aero_gasaerexch_sub error' )
297 : end if
298 : end do ! jsoa
299 :
300 : ! set tendency flags
301 1489176 : dotend(:) = .false.
302 1489176 : dotendqqcw(:) = .false.
303 7445880 : ido_so4a(:) = 0
304 7445880 : ido_nh4a(:) = 0
305 8935056 : ido_soaa(:,:) = 0
306 :
307 1489176 : dotend(l_so4g) = .true.
308 1489176 : if ( do_nh4g ) dotend(l_nh4g) = .true.
309 1489176 : if ( do_msag ) dotend(l_msag) = .true.
310 2978352 : do jsoa = 1, nsoa
311 2978352 : if ( do_soag(jsoa) ) dotend(l_soag(jsoa)) = .true.
312 : end do
313 :
314 1489176 : ntot_soamode = 0
315 7445880 : do n = 1, ntot_amode
316 5956704 : l = lptr_so4_a_amode(n)-loffset
317 5956704 : if ((l > 0) .and. (l <= pcnstxx)) then
318 4467528 : dotend(l) = .true.
319 4467528 : ido_so4a(n) = 1
320 4467528 : if ( do_nh4g ) then
321 0 : l = lptr_nh4_a_amode(n)-loffset
322 0 : if ((l > 0) .and. (l <= pcnstxx)) then
323 0 : dotend(l) = .true.
324 0 : ido_nh4a(n) = 1
325 : end if
326 : end if
327 : end if
328 :
329 13402584 : do jsoa = 1, nsoa
330 11913408 : if ( do_soag(jsoa) ) then
331 5956704 : l = lptr2_soa_a_amode(n,jsoa)-loffset
332 5956704 : if ((l > 0) .and. (l <= pcnstxx)) then
333 2978352 : dotend(l) = .true.
334 2978352 : ido_soaa(n,jsoa) = 1
335 2978352 : ntot_soamode = n
336 : end if
337 : end if
338 : end do ! jsoa
339 : end do ! n
340 :
341 :
342 1489176 : if ( do_soag_any ) ntot_soamode = max( ntot_soamode, modefrm_pcage )
343 :
344 1489176 : if (modefrm_pcage > 0) then
345 1489176 : ido_so4a(modefrm_pcage) = 2
346 1489176 : if (ido_nh4a(modetoo_pcage) == 1) ido_nh4a(modefrm_pcage) = 2
347 2978352 : do jsoa = 1, nsoa
348 2978352 : if (ido_soaa(modetoo_pcage,jsoa) == 1) ido_soaa(modefrm_pcage,jsoa) = 2
349 : end do
350 5956704 : do iq = 1, nspecfrm_pcage
351 4467528 : lsfrm = lspecfrm_pcage(iq)-loffset
352 4467528 : lstoo = lspectoo_pcage(iq)-loffset
353 5956704 : if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
354 4467528 : dotend(lsfrm) = .true.
355 4467528 : if ((lstoo > 0) .and. (lstoo <= pcnst)) then
356 4467528 : dotend(lstoo) = .true.
357 : end if
358 : end if
359 : end do
360 :
361 :
362 1489176 : n = modeptr_pcarbon
363 1489176 : fac_volsfc_pcarbon = exp( 2.5_r8*(alnsg_amode(n)**2) )
364 1489176 : xferfrac_max = 1.0_r8 - 10.0_r8*epsilon(1.0_r8) ! 1-eps
365 : end if
366 :
367 :
368 : ! zero out tendencies and other
369 71735685840 : dqdt(:,:,:) = 0.0_r8
370 71735685840 : dqqcwdt(:,:,:) = 0.0_r8
371 1489176 : qsrflx(:,:,:) = 0.0_r8
372 1489176 : qqcwsrflx(:,:,:) = 0.0_r8
373 :
374 : !-------Initialize evap/cond diagnostics (ncols x pver)-----------
375 1489176 : qconff(:,:) = 0.0_r8
376 1489176 : qevapff(:,:) = 0.0_r8
377 1489176 : qconbb(:,:) = 0.0_r8
378 1489176 : qevapbb(:,:) = 0.0_r8
379 1489176 : qconbg(:,:) = 0.0_r8
380 1489176 : qevapbg(:,:) = 0.0_r8
381 1489176 : qcon(:,:) = 0.0_r8
382 1489176 : qevap(:,:) = 0.0_r8
383 : !---------------------------------------------------
384 :
385 : ! compute gas-to-aerosol mass transfer rates
386 : call gas_aer_uptkrates( ncol, loffset, &
387 : q, t, pmid, &
388 1489176 : dgncur_awet, uptkrate )
389 :
390 :
391 : ! use this for tendency calcs to avoid generating very small negative values
392 1489176 : deltatxx = deltat * (1.0_r8 + 1.0e-15_r8)
393 :
394 :
395 1489176 : jsrf = jsrflx_gaexch
396 139982544 : do k=top_lev,pver
397 2314006344 : do i=1,ncol
398 :
399 : ! fgain_so4(n) = fraction of total h2so4 uptake going to mode n
400 : ! fgain_nh4(n) = fraction of total nh3 uptake going to mode n
401 4348047600 : sum_uprt_so4 = 0.0_r8
402 4348047600 : sum_uprt_nh4 = 0.0_r8
403 4348047600 : sum_uprt_soa = 0.0_r8
404 10870119000 : do n = 1, ntot_amode
405 8696095200 : uptkratebb(n) = uptkrate(n,i,k)
406 8696095200 : if (ido_so4a(n) > 0) then
407 8696095200 : fgain_so4(n) = uptkratebb(n)
408 8696095200 : sum_uprt_so4 = sum_uprt_so4 + fgain_so4(n)
409 8696095200 : if (ido_so4a(n) == 1) then
410 6522071400 : qold_so4(n) = q(i,k,lptr_so4_a_amode(n)-loffset)
411 : else
412 2174023800 : qold_so4(n) = 0.0_r8
413 : end if
414 : else
415 0 : fgain_so4(n) = 0.0_r8
416 0 : qold_so4(n) = 0.0_r8
417 : end if
418 :
419 8696095200 : if (ido_nh4a(n) > 0) then
420 : ! 2.08 factor is for gas diffusivity (nh3/h2so4)
421 : ! differences in fuch-sutugin and accom coef ignored
422 0 : fgain_nh4(n) = uptkratebb(n)*2.08_r8
423 0 : sum_uprt_nh4 = sum_uprt_nh4 + fgain_nh4(n)
424 0 : if (ido_nh4a(n) == 1) then
425 0 : qold_nh4(n) = q(i,k,lptr_nh4_a_amode(n)-loffset)
426 : else
427 0 : qold_nh4(n) = 0.0_r8
428 : end if
429 : else
430 8696095200 : fgain_nh4(n) = 0.0_r8
431 8696095200 : qold_nh4(n) = 0.0_r8
432 : end if
433 :
434 17392190400 : do j = 1, npoa
435 8696095200 : l = lptr2_pom_a_amode(n,j)-loffset
436 17392190400 : if (l > 0) then
437 4348047600 : qold_poa(n,j) = q(i,k,l)
438 : else
439 4348047600 : qold_poa(n,j) = 0.0_r8
440 : end if
441 : end do
442 :
443 8696095200 : itmpa = 0
444 17392190400 : do jsoa = 1, nsoa
445 8696095200 : if (ido_soaa(n,jsoa) > 0) then
446 : ! 0.81 factor is for gas diffusivity (soa/h2so4)
447 : ! (differences in fuch-sutugin and accom coef ignored)
448 6522071400 : fgain_soa(n,jsoa) = uptkratebb(n)*0.81_r8
449 6522071400 : sum_uprt_soa(jsoa) = sum_uprt_soa(jsoa) + fgain_soa(n,jsoa)
450 6522071400 : if (ido_soaa(n,jsoa) == 1) then
451 4348047600 : l = lptr2_soa_a_amode(n,jsoa)-loffset
452 4348047600 : qold_soa(n,jsoa) = q(i,k,l)
453 4348047600 : itmpa = itmpa + 1
454 : else
455 2174023800 : qold_soa(n,jsoa) = 0.0_r8
456 : end if
457 : else
458 2174023800 : fgain_soa(n,jsoa) = 0.0_r8
459 2174023800 : qold_soa(n,jsoa) = 0.0_r8
460 : end if
461 17392190400 : uptkrate_soa(n,jsoa) = fgain_soa(n,jsoa)
462 : end do ! jsoa
463 : ! in previous code versions with nsoa=1,
464 : ! qold_poa was non-zero (i.e., loaded from q) only when ido_soaa(n)=1
465 : ! thus qold_poa=0 for the primary carbon mode which has ido_soaa=2
466 : ! this is probably not how it should be
467 15218166600 : if (itmpa == 0) qold_poa(n,:) = 0.0_r8
468 :
469 : end do ! n
470 :
471 2174023800 : if (sum_uprt_so4 > 0.0_r8) then
472 10870119000 : do n = 1, ntot_amode
473 10870119000 : fgain_so4(n) = fgain_so4(n) / sum_uprt_so4
474 : end do
475 : end if
476 : ! at this point (sum_uprt_so4 <= 0.0) only when all the fgain_so4 are zero
477 2174023800 : if (sum_uprt_nh4 > 0.0_r8) then
478 0 : do n = 1, ntot_amode
479 0 : fgain_nh4(n) = fgain_nh4(n) / sum_uprt_nh4
480 : end do
481 : end if
482 :
483 4348047600 : do jsoa = 1, nsoa
484 4348047600 : if (sum_uprt_soa(jsoa) > 0.0_r8) then
485 10870119000 : do n = 1, ntot_amode
486 10870119000 : fgain_soa(n,jsoa) = fgain_soa(n,jsoa) / sum_uprt_soa(jsoa)
487 : end do
488 : end if
489 : end do
490 :
491 : ! uptake amount (fraction of gas uptaken) over deltat
492 2174023800 : avg_uprt_so4 = (1.0_r8 - exp(-deltatxx*sum_uprt_so4))/deltatxx
493 2174023800 : avg_uprt_nh4 = (1.0_r8 - exp(-deltatxx*sum_uprt_nh4))/deltatxx
494 :
495 4348047600 : do jsoa = 1, nsoa
496 4348047600 : avg_uprt_soa(jsoa) = (1.0_r8 - exp(-deltatxx*sum_uprt_soa(jsoa)))/deltatxx
497 : end do
498 :
499 : ! sum_dqdt_so4 = so4_a tendency from h2so4 gas uptake (mol/mol/s)
500 : ! sum_dqdt_msa = msa_a tendency from msa gas uptake (mol/mol/s)
501 : ! sum_dqdt_nh4 = nh4_a tendency from nh3 gas uptake (mol/mol/s)
502 : ! sum_dqdt_soa = soa_a tendency from soa gas uptake (mol/mol/s)
503 2174023800 : sum_dqdt_so4 = q(i,k,l_so4g) * avg_uprt_so4
504 2174023800 : if ( do_msag ) then
505 0 : sum_dqdt_msa = q(i,k,l_msag) * avg_uprt_so4
506 : else
507 : sum_dqdt_msa = 0.0_r8
508 : end if
509 2174023800 : if ( do_nh4g ) then
510 0 : sum_dqdt_nh4 = q(i,k,l_nh4g) * avg_uprt_nh4
511 : else
512 : sum_dqdt_nh4 = 0.0_r8
513 : end if
514 :
515 4348047600 : do jsoa = 1, nsoa
516 4348047600 : if ( do_soag(jsoa) ) then
517 2174023800 : sum_dqdt_soa(jsoa) = q(i,k,l_soag(jsoa)) * avg_uprt_soa(jsoa)
518 : else
519 0 : sum_dqdt_soa(jsoa) = 0.0_r8
520 : end if
521 : end do
522 :
523 2174023800 : if ( associated(sulfeq) .and. (k <= troplev(i)) ) then
524 : ! compute TMR tendencies for so4 interstial aerosol due to reversible gas uptake
525 : ! only above the tropopause
526 :
527 0 : tmp_kxt = deltatxx*sum_uprt_so4 ! sum over modes of uptake_rate*deltat
528 0 : tmp_pxt = 0.0_r8
529 0 : do n = 1, ntot_amode
530 0 : if (ido_so4a(n) <= 0) cycle
531 0 : tmp_pxt = tmp_pxt + uptkratebb(n)*sulfeq(i,k,n)
532 : end do
533 0 : tmp_pxt = max( 0.0_r8, tmp_pxt*deltatxx ) ! sum over modes of uptake_rate*sulfeq*deltat
534 0 : tmp_so4g_bgn = q(i,k,l_so4g)
535 : ! calc avg h2so4(g) over deltat
536 0 : if (tmp_kxt >= 1.0e-5_r8) then
537 : ! exponential decay towards equilibrium value solution
538 0 : tmp_so4g_equ = tmp_pxt/tmp_kxt
539 0 : tmp_so4g_avg = tmp_so4g_equ + (tmp_so4g_bgn-tmp_so4g_equ)*(1.0_r8-exp(-tmp_kxt))/tmp_kxt
540 : else
541 : ! first order approx for tmp_kxt small
542 0 : tmp_so4g_avg = tmp_so4g_bgn*(1.0_r8-0.5_r8*tmp_kxt) + 0.5_r8*tmp_pxt
543 : end if
544 0 : sum_dqdt_so4 = 0.0_r8
545 0 : do n = 1, ntot_amode
546 0 : if (ido_so4a(n) <= 0) cycle
547 : ! calc change to so4(a) in mode n
548 0 : if (ido_so4a(n) == 1) then
549 0 : l = lptr_so4_a_amode(n)-loffset
550 0 : tmp_so4a_bgn = q(i,k,l)
551 : else
552 : tmp_so4a_bgn = 0.0_r8
553 : end if
554 0 : tmp_so4a_end = tmp_so4a_bgn + deltatxx*uptkratebb(n)*(tmp_so4g_avg-sulfeq(i,k,n))
555 0 : tmp_so4a_end = max( 0.0_r8, tmp_so4a_end )
556 0 : dqdt_so4(n) = (tmp_so4a_end - tmp_so4a_bgn)/deltatxx
557 0 : sum_dqdt_so4 = sum_dqdt_so4 + dqdt_so4(n)
558 : end do
559 : ! do not allow msa condensation in stratosphere
560 : ! ( Note that the code for msa has never been used.
561 : ! The plan was to simulate msa(g), treat it as non-volatile (like h2so4(g)),
562 : ! and treat condensed msa as sulfate, so just one additional tracer. )
563 0 : if ( do_msag ) sum_dqdt_msa = 0.0_r8
564 :
565 : else
566 : ! compute TMR tendencies for so4 interstial aerosol due to simple gas uptake
567 10870119000 : do n = 1, ntot_amode
568 10870119000 : dqdt_so4(n) = fgain_so4(n)*(sum_dqdt_so4 + sum_dqdt_msa)
569 : end do
570 : end if
571 :
572 : ! compute TMR tendencies for nh4 interstial aerosol due to simple gas uptake
573 : ! but force nh4/so4 molar ratio <= 2
574 2174023800 : sum_dqdt_nh4_b = 0.0_r8
575 10870119000 : dqdt_nh4(:) = 0._r8
576 2174023800 : if ( do_nh4g ) then
577 0 : do n = 1, ntot_amode
578 0 : dqdt_nh4(n) = fgain_nh4(n)*sum_dqdt_nh4
579 0 : qnew_nh4 = qold_nh4(n) + dqdt_nh4(n)*deltat
580 0 : qnew_so4 = qold_so4(n) + dqdt_so4(n)*deltat
581 0 : qmax_nh4 = 2.0_r8*qnew_so4
582 0 : if (qnew_nh4 > qmax_nh4) then
583 0 : dqdt_nh4(n) = (qmax_nh4 - qold_nh4(n))/deltatxx
584 : end if
585 0 : sum_dqdt_nh4_b = sum_dqdt_nh4_b + dqdt_nh4(n)
586 : end do
587 : end if
588 :
589 2174023800 : if (( do_soag_any ) .and. (method_soa > 1)) then
590 : ! compute TMR tendencies for soag and soa interstial aerosol
591 : ! using soa parameterization
592 2174023800 : niter_max = 1000
593 13044142800 : dqdt_soa(:,:) = 0.0_r8
594 4348047600 : dqdt_soag(:) = 0.0_r8
595 4348047600 : do jsoa = 1, nsoa
596 4348047600 : qold_soag(jsoa) = q(i,k,l_soag(jsoa))
597 : end do
598 :
599 : ! get molecular weight from the host model
600 10870119000 : do n = 1, ntot_amode
601 43480476000 : do l = 1, nspec_amode(n)
602 32610357000 : call rad_cnst_get_info(0, n, l, spec_type=spec_type )
603 8696095200 : select case( spec_type )
604 : case('s-organic')
605 8696095200 : mw_soa_host(:) = specmw_amode(l,n)
606 : case('p-organic')
607 36958404600 : mw_poa_host(:) = specmw_amode(l,n)
608 : end select
609 : end do
610 : end do
611 :
612 : call modal_aero_soaexch( deltat, t(i,k), pmid(i,k), &
613 : niter, niter_max, ntot_amode, ntot_soamode, npoa, nsoa, &
614 : mw_poa_host, mw_soa_host, &
615 : qold_soag, qold_soa, qold_poa, uptkrate_soa, &
616 2174023800 : dqdt_soag, dqdt_soa )
617 4348047600 : sum_dqdt_soa(:) = -dqdt_soag(:)
618 :
619 : else if ( do_soag_any ) then
620 : ! compute TMR tendencies for soa interstial aerosol
621 : ! due to simple gas uptake
622 :
623 : do jsoa = 1, nsoa
624 : do n = 1, ntot_amode
625 : dqdt_soa(n,jsoa) = fgain_soa(n,jsoa)*sum_dqdt_soa(jsoa)
626 : end do
627 : end do
628 : else ! method_soa is neither 1 nor 2, no uptake
629 0 : dqdt_soa(:,:) = 0.0_r8
630 : end if
631 :
632 2174023800 : pdel_fac = pdel(i,k)/gravit
633 10870119000 : do n = 1, ntot_amode
634 8696095200 : if (ido_so4a(n) == 1) then
635 6522071400 : l = lptr_so4_a_amode(n)-loffset
636 6522071400 : dqdt(i,k,l) = dqdt_so4(n)
637 6522071400 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_so4(n)*pdel_fac
638 : end if
639 :
640 8696095200 : if ( do_nh4g ) then
641 0 : if (ido_nh4a(n) == 1) then
642 0 : l = lptr_nh4_a_amode(n)-loffset
643 0 : dqdt(i,k,l) = dqdt_nh4(n)
644 0 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_nh4(n)*pdel_fac
645 : end if
646 : end if
647 :
648 19566214200 : do jsoa = 1, nsoa
649 17392190400 : if ( do_soag(jsoa) ) then
650 8696095200 : if (ido_soaa(n,jsoa) == 1) then
651 4348047600 : l = lptr2_soa_a_amode(n,jsoa)-loffset
652 4348047600 : dqdt(i,k,l) = dqdt_soa(n,jsoa) !calculated by modal_aero_soaexch for method_soa=2
653 4348047600 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_soa(n,jsoa)*pdel_fac
654 : !------- Add code for condensation/evaporation diagnostics---
655 4348047600 : if (nsoa.eq.15) then !check for current SOA package
656 0 : if(jsoa.ge.1.and.jsoa.le.5) then ! Fossil SOA species
657 0 : if (dqdt_soa(n,jsoa).ge.0.0_r8) then
658 0 : qconff(i,k)=qconff(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
659 0 : elseif(dqdt_soa(n,jsoa).lt.0.0_r8) then
660 0 : qevapff(i,k)=qevapff(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
661 : endif
662 :
663 0 : elseif(jsoa.ge.6.and.jsoa.le.10) then ! Biomass SOA species
664 0 : if (dqdt_soa(n,jsoa).ge.0.0_r8) then
665 0 : qconbb(i,k)=qconbb(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
666 0 : elseif(dqdt_soa(n,jsoa).lt.0.0_r8) then
667 0 : qevapbb(i,k)=qevapbb(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
668 : endif
669 :
670 0 : elseif(jsoa.ge.11.and.jsoa.le.15) then ! Biomass SOA species
671 0 : if (dqdt_soa(n,jsoa).ge.0.0_r8) then
672 0 : qconbg(i,k)=qconbg(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
673 0 : elseif(dqdt_soa(n,jsoa).lt.0.0_r8) then
674 0 : qevapbg(i,k)=qevapbg(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
675 : endif
676 :
677 : endif ! jsoa
678 : endif !nsoa
679 4348047600 : if (nsoa.eq.5) then !check for current SOA package
680 0 : if (dqdt_soa(n,jsoa).ge.0.0_r8) then
681 0 : qcon(i,k)=qcon(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
682 0 : elseif(dqdt_soa(n,jsoa).lt.0.0_r8) then
683 0 : qevap(i,k)=qevap(i,k)+dqdt_soa(n,jsoa)*(adv_mass(l)/mwdry)
684 : endif
685 : endif !nsoa
686 : !---------------------------------------------------------------------------------------------------------------------
687 : end if
688 : end if
689 : end do
690 : end do ! n
691 :
692 : ! compute TMR tendencies for h2so4, nh3, and msa gas
693 : ! due to simple gas uptake
694 2174023800 : l = l_so4g
695 2174023800 : dqdt(i,k,l) = -sum_dqdt_so4
696 2174023800 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
697 :
698 2174023800 : if ( do_msag ) then
699 0 : l = l_msag
700 0 : dqdt(i,k,l) = -sum_dqdt_msa
701 0 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
702 : end if
703 :
704 2174023800 : if ( do_nh4g ) then
705 0 : l = l_nh4g
706 0 : dqdt(i,k,l) = -sum_dqdt_nh4_b
707 0 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
708 : end if
709 :
710 4348047600 : do jsoa = 1, nsoa
711 4348047600 : if ( do_soag(jsoa) ) then
712 2174023800 : l = l_soag(jsoa)
713 2174023800 : dqdt(i,k,l) = -sum_dqdt_soa(jsoa)
714 : ! dqdt for gas is negative of the sum of dqdt for aerosol soa species in each mode: Manish
715 2174023800 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
716 : end if
717 : end do
718 :
719 : ! compute TMR tendencies associated with primary carbon aging
720 2312517168 : if (modefrm_pcage > 0) then
721 2174023800 : n = modeptr_pcarbon
722 2174023800 : tmpa = 0.0_r8
723 4348047600 : do jsoa = 1, nsoa
724 4348047600 : tmpa = tmpa + dqdt_soa(n,jsoa)*fac_m2v_soa(jsoa)*soa_equivso4_factor(jsoa)
725 : end do
726 : vol_shell = deltat * &
727 2174023800 : ( dqdt_so4(n)*fac_m2v_so4 + dqdt_nh4(n)*fac_m2v_nh4 + tmpa )
728 2174023800 : vol_core = 0.0_r8
729 6522071400 : do l = 1, nspec_amode(n)
730 : vol_core = vol_core + &
731 6522071400 : q(i,k,lmassptr_amode(l,n)-loffset)*fac_m2v_pcarbon(l)
732 : end do
733 : ! ratio1 = vol_shell/vol_core =
734 : ! actual hygroscopic-shell-volume/carbon-core-volume after gas uptake
735 : ! ratio2 = 6.0_r8*dr_so4_monolayers_pcage/(dgncur_a*fac_volsfc_pcarbon)
736 : ! = (shell-volume corresponding to n_so4_monolayers_pcage)/core-volume
737 : ! The 6.0/(dgncur_a*fac_volsfc_pcarbon) = (mode-surface-area/mode-volume)
738 : ! Note that vol_shell includes both so4+nh4 AND soa as "equivalent so4",
739 : ! The soa_equivso4_factor accounts for the lower hygroscopicity of soa.
740 : !
741 : ! Define xferfrac_pcage = min( 1.0, ratio1/ratio2)
742 : ! But ratio1/ratio2 == tmp1/tmp2, and coding below avoids possible overflow
743 : !
744 2174023800 : tmp1 = vol_shell*dgncur_a(i,k,n)*fac_volsfc_pcarbon
745 2174023800 : tmp2 = max( 6.0_r8*dr_so4_monolayers_pcage*vol_core, 0.0_r8 )
746 2174023800 : if (tmp1 >= tmp2) then
747 : xferfrac_pcage = xferfrac_max
748 : else
749 2173297671 : xferfrac_pcage = min( tmp1/tmp2, xferfrac_max )
750 : end if
751 :
752 2174023800 : if (xferfrac_pcage > 0.0_r8) then
753 8696095200 : do iq = 1, nspecfrm_pcage
754 6522071400 : lsfrm = lspecfrm_pcage(iq)-loffset
755 6522071400 : lstoo = lspectoo_pcage(iq)-loffset
756 6522071400 : xferrate = (xferfrac_pcage/deltat)*q(i,k,lsfrm)
757 6522071400 : dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferrate
758 6522071400 : qsrflx(i,lsfrm,jsrf) = qsrflx(i,lsfrm,jsrf) - xferrate*pdel_fac
759 8696095200 : if ((lstoo > 0) .and. (lstoo <= pcnst)) then
760 6522071400 : dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferrate
761 6522071400 : qsrflx(i,lstoo,jsrf) = qsrflx(i,lstoo,jsrf) + xferrate*pdel_fac
762 : end if
763 : end do
764 :
765 2174023800 : if (ido_so4a(modetoo_pcage) > 0) then
766 2174023800 : l = lptr_so4_a_amode(modetoo_pcage)-loffset
767 2174023800 : dqdt(i,k,l) = dqdt(i,k,l) + dqdt_so4(modefrm_pcage)
768 2174023800 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_so4(modefrm_pcage)*pdel_fac
769 : end if
770 :
771 2174023800 : if (ido_nh4a(modetoo_pcage) > 0) then
772 0 : l = lptr_nh4_a_amode(modetoo_pcage)-loffset
773 0 : dqdt(i,k,l) = dqdt(i,k,l) + dqdt_nh4(modefrm_pcage)
774 0 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_nh4(modefrm_pcage)*pdel_fac
775 : end if
776 :
777 4348047600 : do jsoa = 1, nsoa
778 4348047600 : if (ido_soaa(modetoo_pcage,jsoa) > 0) then
779 2174023800 : l = lptr2_soa_a_amode(modetoo_pcage,jsoa)-loffset
780 2174023800 : dqdt(i,k,l) = dqdt(i,k,l) + dqdt_soa(modefrm_pcage,jsoa)
781 2174023800 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_soa(modefrm_pcage,jsoa)*pdel_fac
782 : end if
783 : end do
784 :
785 : end if
786 :
787 : end if
788 :
789 : end do ! "i = 1, ncol"
790 : end do ! "k = top_lev, pver"
791 :
792 : ! set "temporary testing arrays"
793 71735685840 : qold(:,:,:) = q(:,:,:)
794 71735685840 : qqcwold(:,:,:) = qqcw(:,:,:)
795 71735685840 : dqdtsv1(:,:,:) = dqdt(:,:,:)
796 71735685840 : dqqcwdtsv1(:,:,:) = dqqcwdt(:,:,:)
797 :
798 :
799 : !
800 : ! do renaming calcs
801 : !
802 1489176 : dotendrn(:) = .false.
803 1489176 : dotendqqcwrn(:) = .false.
804 2314006344 : dorename_atik(1:ncol,:) = .true.
805 1489176 : is_dorename_atik = .true.
806 : call modal_aero_rename_sub( &
807 : 'modal_aero_gasaerexch_sub', &
808 : lchnk, ncol, nstep, &
809 : loffset, deltat, &
810 : pdel, troplev, &
811 : dotendrn, q, &
812 : dqdt, dqdt_other, &
813 : dotendqqcwrn, qqcw, &
814 : dqqcwdt, dqqcwdt_other, &
815 : is_dorename_atik, dorename_atik, &
816 : jsrflx_rename, nsrflx, &
817 1489176 : qsrflx, qqcwsrflx )
818 :
819 :
820 : ! This applies dqdt tendencies for all species
821 : ! apply the dqdt to update q (and same for qqcw)
822 : !
823 47653632 : do l = 1, pcnstxx
824 46164456 : if ( dotend(l) .or. dotendrn(l) ) then
825 2939633424 : do k = top_lev, pver
826 48594133224 : do i = 1, ncol
827 48562860528 : q(i,k,l) = q(i,k,l) + dqdt(i,k,l)*deltat
828 : end do
829 : end do
830 : end if
831 47653632 : if ( dotendqqcw(l) .or. dotendqqcwrn(l) ) then
832 1959755616 : do k = top_lev, pver
833 32396088816 : do i = 1, ncol
834 32375240352 : qqcw(i,k,l) = qqcw(i,k,l) + dqqcwdt(i,k,l)*deltat
835 : end do
836 : end do
837 : end if
838 : end do
839 :
840 : ! diagnostics start -------------------------------------------------------
841 : !!$ if (ldiag3 > 0) then
842 : !!$ if (icol_diag > 0) then
843 : !!$ i = icol_diag
844 : !!$ write(*,'(a,3i5)') 'gasaerexch ppp nstep,lat,lon', nstep, latndx(i), lonndx(i)
845 : !!$ write(*,'(2i5,3(2x,a))') 0, 0, 'ppp', 'pdel for all k'
846 : !!$ write(*,'(1p,7e12.4)') (pdel(i,k), k=top_lev,pver)
847 : !!$
848 : !!$ write(*,'(a,3i5)') 'gasaerexch ddd nstep,lat,lon', nstep, latndx(i), lonndx(i)
849 : !!$ do l = 1, pcnstxx
850 : !!$ lb = l + loffset
851 : !!$
852 : !!$ if ( dotend(l) .or. dotendrn(l) ) then
853 : !!$ write(*,'(2i5,3(2x,a))') 1, l, 'ddd1', cnst_name(lb), 'qold for all k'
854 : !!$ write(*,'(1p,7e12.4)') (qold(i,k,l), k=top_lev,pver)
855 : !!$ write(*,'(2i5,3(2x,a))') 1, l, 'ddd2', cnst_name(lb), 'qnew for all k'
856 : !!$ write(*,'(1p,7e12.4)') (q(i,k,l), k=top_lev,pver)
857 : !!$ write(*,'(2i5,3(2x,a))') 1, l, 'ddd3', cnst_name(lb), 'dqdt from conden for all k'
858 : !!$ write(*,'(1p,7e12.4)') (dqdtsv1(i,k,l), k=top_lev,pver)
859 : !!$ write(*,'(2i5,3(2x,a))') 1, l, 'ddd4', cnst_name(lb), 'dqdt from rename for all k'
860 : !!$ write(*,'(1p,7e12.4)') ((dqdt(i,k,l)-dqdtsv1(i,k,l)), k=top_lev,pver)
861 : !!$ write(*,'(2i5,3(2x,a))') 1, l, 'ddd5', cnst_name(lb), 'dqdt other for all k'
862 : !!$ write(*,'(1p,7e12.4)') (dqdt_other(i,k,l), k=top_lev,pver)
863 : !!$ end if
864 : !!$
865 : !!$ if ( dotendqqcw(l) .or. dotendqqcwrn(l) ) then
866 : !!$ write(*,'(2i5,3(2x,a))') 2, l, 'ddd1', cnst_name_cw(lb), 'qold for all k'
867 : !!$ write(*,'(1p,7e12.4)') (qqcwold(i,k,l), k=top_lev,pver)
868 : !!$ write(*,'(2i5,3(2x,a))') 2, l, 'ddd2', cnst_name_cw(lb), 'qnew for all k'
869 : !!$ write(*,'(1p,7e12.4)') (qqcw(i,k,l), k=top_lev,pver)
870 : !!$ write(*,'(2i5,3(2x,a))') 2, l, 'ddd3', cnst_name_cw(lb), 'dqdt from conden for all k'
871 : !!$ write(*,'(1p,7e12.4)') (dqqcwdtsv1(i,k,l), k=top_lev,pver)
872 : !!$ write(*,'(2i5,3(2x,a))') 2, l, 'ddd4', cnst_name_cw(lb), 'dqdt from rename for all k'
873 : !!$ write(*,'(1p,7e12.4)') ((dqqcwdt(i,k,l)-dqqcwdtsv1(i,k,l)), k=top_lev,pver)
874 : !!$ write(*,'(2i5,3(2x,a))') 2, l, 'ddd5', cnst_name_cw(lb), 'dqdt other for all k'
875 : !!$ write(*,'(1p,7e12.4)') (dqqcwdt_other(i,k,l), k=top_lev,pver)
876 : !!$ end if
877 : !!$
878 : !!$ end do
879 : !!$
880 : !!$ write(*,'(a,3i5)') 'gasaerexch fff nstep,lat,lon', nstep, latndx(i), lonndx(i)
881 : !!$ do l = 1, pcnstxx
882 : !!$ lb = l + loffset
883 : !!$ if ( dotend(l) .or. dotendrn(l) .or. dotendqqcw(l) .or. dotendqqcwrn(l) ) then
884 : !!$ write(*,'(i5,2(2x,a,2l3))') l, &
885 : !!$ cnst_name(lb), dotend(l), dotendrn(l), &
886 : !!$ cnst_name_cw(lb), dotendqqcw(l), dotendqqcwrn(l)
887 : !!$ end if
888 : !!$ end do
889 : !!$
890 : !!$ end if
891 : !!$ end if
892 : ! diagnostics end ---------------------------------------------------------
893 :
894 : !-----Outfld for condensation/evaporation------------------------------
895 1489176 : if (nsoa.eq.5) then !check for current SOA package
896 0 : call outfld(trim('qcon_gaex'), qcon(:,:), pcols, lchnk )
897 0 : call outfld(trim('qevap_gaex'), qevap(:,:), pcols, lchnk )
898 : endif
899 : !-----------------------------------------------------------------------
900 1489176 : if (nsoa.eq.15) then !check for current SOA package
901 0 : call outfld(trim('qconff_gaex'), qconff(:,:), pcols, lchnk )
902 0 : call outfld(trim('qevapff_gaex'), qevapff(:,:), pcols, lchnk )
903 0 : call outfld(trim('qconbb_gaex'), qconbb(:,:), pcols, lchnk )
904 0 : call outfld(trim('qevapbb_gaex'), qevapbb(:,:), pcols, lchnk )
905 0 : call outfld(trim('qconbg_gaex'), qconbg(:,:), pcols, lchnk )
906 0 : call outfld(trim('qevapbg_gaex'), qevapbg(:,:), pcols, lchnk )
907 : endif
908 : !-----------------------------------------------------------------------
909 :
910 : ! do history file column-tendency fields
911 47653632 : do l = 1, pcnstxx
912 46164456 : lb = l + loffset
913 139982544 : do jsrf = 1, 2
914 323151192 : do jac = 1, 2
915 276986736 : if (jac == 1) then
916 92328912 : if (jsrf == jsrflx_gaexch) then
917 46164456 : if ( .not. dotend(l) ) cycle
918 19359288 : fieldname = trim(cnst_name(lb)) // '_sfgaex1'
919 46164456 : else if (jsrf == jsrflx_rename) then
920 46164456 : if ( .not. dotendrn(l) ) cycle
921 20848464 : fieldname = trim(cnst_name(lb)) // '_sfgaex2'
922 : else
923 : cycle
924 : end if
925 671375952 : do i = 1, ncol
926 671375952 : qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf)*(adv_mass(l)/mwdry)
927 : end do
928 40207752 : call outfld( fieldname, qsrflx(:,l,jsrf), pcols, lchnk )
929 : else
930 92328912 : if (jsrf == jsrflx_gaexch) then
931 : cycle
932 46164456 : else if (jsrf == jsrflx_rename) then
933 46164456 : if ( .not. dotendqqcwrn(l) ) cycle
934 20848464 : fieldname = trim(cnst_name_cw(lb)) // '_sfgaex2'
935 : else
936 : cycle
937 : end if
938 348120864 : do i = 1, ncol
939 348120864 : qqcwsrflx(i,l,jsrf) = qqcwsrflx(i,l,jsrf)*(adv_mass(l)/mwdry)
940 : end do
941 20848464 : call outfld( fieldname, qqcwsrflx(:,l,jsrf), pcols, lchnk )
942 : end if
943 : end do ! jac = ...
944 : end do ! jsrf = ...
945 : end do ! l = ...
946 :
947 1489176 : return
948 1489176 : end subroutine modal_aero_gasaerexch_sub
949 :
950 :
951 : !----------------------------------------------------------------------
952 : !----------------------------------------------------------------------
953 1489176 : subroutine gas_aer_uptkrates( ncol, loffset, &
954 1489176 : q, t, pmid, &
955 1489176 : dgncur_awet, uptkrate )
956 :
957 : !
958 : ! /
959 : ! computes uptkrate = | dx dN/dx gas_conden_rate(Dp(x))
960 : ! /
961 : ! using Gauss-Hermite quadrature of order nghq=2
962 : !
963 : ! Dp = particle diameter (cm)
964 : ! x = ln(Dp)
965 : ! dN/dx = log-normal particle number density distribution
966 : ! gas_conden_rate(Dp) = 2 * pi * gasdiffus * Dp * F(Kn,ac)
967 : ! F(Kn,ac) = Fuchs-Sutugin correction factor
968 : ! Kn = Knudsen number
969 : ! ac = accomodation coefficient
970 : !
971 :
972 1489176 : use physconst, only: mwdry, rair
973 :
974 : implicit none
975 :
976 :
977 : integer, intent(in) :: ncol ! number of atmospheric column
978 : integer, intent(in) :: loffset
979 : real(r8), intent(in) :: q(ncol,pver,pcnstxx) ! Tracer array (mol,#/mol-air)
980 : real(r8), intent(in) :: t(pcols,pver) ! Temperature in Kelvin
981 : real(r8), intent(in) :: pmid(pcols,pver) ! Air pressure in Pa
982 : real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)
983 :
984 : real(r8), intent(out) :: uptkrate(ntot_amode,pcols,pver)
985 : ! gas-to-aerosol mass transfer rates (1/s)
986 :
987 :
988 : ! local
989 : integer, parameter :: nghq = 2
990 : integer :: i, iq, k, l1, l2, la, n
991 :
992 : ! Can use sqrt here once Lahey is gone.
993 : real(r8), parameter :: tworootpi = 3.5449077_r8
994 : real(r8), parameter :: root2 = 1.4142135_r8
995 : real(r8), parameter :: beta = 2.0_r8
996 :
997 : real(r8) :: aircon
998 : real(r8) :: const
999 : real(r8) :: dp, dum_m2v
1000 : real(r8) :: dryvol_a(pcols,pver)
1001 : real(r8) :: gasdiffus, gasspeed
1002 : real(r8) :: freepathx2, fuchs_sutugin
1003 : real(r8) :: knudsen
1004 : real(r8) :: lndp, lndpgn, lnsg
1005 : real(r8) :: num_a
1006 : real(r8) :: rhoair
1007 : real(r8) :: sumghq
1008 : real(r8), save :: xghq(nghq), wghq(nghq) ! quadrature abscissae and weights
1009 :
1010 : data xghq / 0.70710678_r8, -0.70710678_r8 /
1011 : data wghq / 0.88622693_r8, 0.88622693_r8 /
1012 :
1013 :
1014 : ! outermost loop over all modes
1015 7445880 : do n = 1, ntot_amode
1016 :
1017 : ! 22-aug-2007 rc easter - get number from q array rather
1018 : ! than computing a "bounded" number conc.
1019 : !! compute dry volume = sum_over_components{ component_mass / density }
1020 : !! (m3-AP/mol-air)
1021 : !! compute it for all i,k to improve accessing q array
1022 : ! dryvol_a(1:ncol,:) = 0.0_r8
1023 : ! do l1 = 1, nspec_amode(n)
1024 : ! l2 = lspectype_amode(l1,n)
1025 : !! dum_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
1026 : !! [m3-AP/kmol-AP]= [kg-AP/kmol-AP] / [kg-AP/m3-AP]
1027 : ! dum_m2v = specmw_amode(l2) / specdens_amode(l2)
1028 : ! la = lmassptr_amode(l1,n)
1029 : ! dryvol_a(1:ncol,:) = dryvol_a(1:ncol,:) &
1030 : ! + max(0.0_r8,q(1:ncol,:,la))*dum_m2v
1031 : ! end do
1032 :
1033 : ! loops k and i
1034 561419352 : do k=top_lev,pver
1035 9256025376 : do i=1,ncol
1036 :
1037 8696095200 : rhoair = pmid(i,k)/(rair*t(i,k)) ! (kg-air/m3)
1038 : ! aircon = 1.0e3*rhoair/mwdry ! (mol-air/m3)
1039 :
1040 : !! "bounded" number conc. (#/m3)
1041 : ! num_a = dryvol_a(i,k)*v2ncur_a(i,k,n)*aircon
1042 :
1043 : ! number conc. (#/m3) -- note q(i,k,numptr) is (#/kmol-air)
1044 : ! so need aircon in (kmol-air/m3)
1045 8696095200 : aircon = rhoair/mwdry ! (kmol-air/m3)
1046 8696095200 : num_a = q(i,k,numptr_amode(n)-loffset)*aircon
1047 :
1048 : ! gasdiffus = h2so4 gas diffusivity from mosaic code (m^2/s)
1049 : ! (pmid must be Pa)
1050 8696095200 : gasdiffus = 0.557e-4_r8 * (t(i,k)**1.75_r8) / pmid(i,k)
1051 : ! gasspeed = h2so4 gas mean molecular speed from mosaic code (m/s)
1052 8696095200 : gasspeed = 1.470e1_r8 * sqrt(t(i,k))
1053 : ! freepathx2 = 2 * (h2so4 mean free path) (m)
1054 8696095200 : freepathx2 = 6.0_r8*gasdiffus/gasspeed
1055 :
1056 8696095200 : lnsg = log( sigmag_amode(n) )
1057 8696095200 : lndpgn = log( dgncur_awet(i,k,n) ) ! (m)
1058 8696095200 : const = tworootpi * num_a * exp(beta*lndpgn + 0.5_r8*(beta*lnsg)**2)
1059 :
1060 : ! sum over gauss-hermite quadrature points
1061 8696095200 : sumghq = 0.0_r8
1062 26088285600 : do iq = 1, nghq
1063 17392190400 : lndp = lndpgn + beta*lnsg**2 + root2*lnsg*xghq(iq)
1064 17392190400 : dp = exp(lndp)
1065 :
1066 : ! knudsen number
1067 17392190400 : knudsen = freepathx2/dp
1068 : ! Changed by Manish Shrivastava on 7/17/2013 to use accom=1; because we do not know better
1069 : ! following assumes accomodation coefficient = ac = 1. instead 0.65 ! answer change needs to be tested
1070 : ! (Adams & Seinfeld, 2002, JGR, and references therein)
1071 : ! fuchs_sutugin = (0.75*ac*(1. + knudsen)) /
1072 : ! (knudsen*(1.0 + knudsen + 0.283*ac) + 0.75*ac)
1073 : fuchs_sutugin = (0.4875_r8*(1._r8 + knudsen)) / &
1074 17392190400 : (knudsen*(1.184_r8 + knudsen) + 0.4875_r8)
1075 26088285600 : sumghq = sumghq + wghq(iq)*dp*fuchs_sutugin/(dp**beta)
1076 : end do
1077 9250068672 : uptkrate(n,i,k) = const * gasdiffus * sumghq
1078 :
1079 : end do ! "do i = 1, ncol"
1080 : end do ! "do k = 1, pver"
1081 :
1082 : end do ! "do n = 1, ntot_soamode"
1083 :
1084 :
1085 1489176 : return
1086 : end subroutine gas_aer_uptkrates
1087 :
1088 : !----------------------------------------------------------------------
1089 :
1090 2174023800 : subroutine modal_aero_soaexch( dtfull, temp, pres, &
1091 : niter, niter_max, ntot_amode, ntot_soamode, ntot_poaspec, ntot_soaspec, &
1092 2174023800 : mw_poa_host, mw_soa_host, &
1093 2174023800 : g_soa_in, a_soa_in, a_poa_in, xferrate_in, &
1094 2174023800 : g_soa_tend, a_soa_tend )
1095 : ! g_soa_tend, a_soa_tend, g0_soa, idiagss )
1096 :
1097 : !-----------------------------------------------------------------------
1098 : !
1099 : ! Purpose:
1100 : !
1101 : ! calculates condensation/evaporation of "soa gas"
1102 : ! to/from multiple aerosol modes in 1 grid cell
1103 : !
1104 : ! key assumptions
1105 : ! (1) ambient equilibrium vapor pressure of soa gas
1106 : ! is given by p0_soa_298 and delh_vap_soa
1107 : ! (2) equilibrium vapor pressure of soa gas at aerosol
1108 : ! particle surface is given by raoults law in the form
1109 : ! g_star = g0_soa*[a_soa/(a_soa + a_opoa)]
1110 : ! (3) (oxidized poa)/(total poa) is equal to frac_opoa (constant)
1111 : !
1112 : !
1113 : ! Author: R. Easter and R. Zaveri
1114 : ! Additions to run with multiple BC, SOA and POM's: Shrivastava et al., 2015
1115 : !-----------------------------------------------------------------------
1116 :
1117 : use mo_constants, only: rgas ! Gas constant (J/K/mol)
1118 :
1119 : implicit none
1120 :
1121 : real(r8), intent(in) :: dtfull ! full integration time step (s)
1122 : real(r8), intent(in) :: temp ! air temperature (K)
1123 : real(r8), intent(in) :: pres ! air pressure (Pa)
1124 : integer, intent(out) :: niter ! number of iterations performed
1125 : integer, intent(in) :: niter_max ! max allowed number of iterations
1126 : integer, intent(in) :: ntot_amode ! number of modes
1127 : integer, intent(in) :: ntot_soamode ! number of modes having soa
1128 : integer, intent(in) :: ntot_poaspec ! number of poa species
1129 : integer, intent(in) :: ntot_soaspec ! number of soa species
1130 : real(r8), intent(in) :: mw_poa_host(ntot_poaspec) ! molec wght of poa used in host code
1131 : real(r8), intent(in) :: mw_soa_host(ntot_soaspec) ! molec wght of poa used in host code
1132 : real(r8), intent(in) :: g_soa_in(ntot_soaspec) ! initial soa gas mixrat (mol/mol at host mw)
1133 : real(r8), intent(in) :: a_soa_in(ntot_amode,ntot_soaspec) ! initial soa aerosol mixrat (mol/mol at host mw)
1134 : real(r8), intent(in) :: a_poa_in(ntot_amode,ntot_poaspec) ! initial poa aerosol mixrat (mol/mol at host mw)
1135 : real(r8), intent(in) :: xferrate_in(ntot_amode,ntot_soaspec) ! gas-aerosol mass transfer rate (1/s)
1136 : real(r8), intent(out) :: g_soa_tend(ntot_soaspec) ! soa gas mixrat tendency (mol/mol/s at host mw)
1137 : real(r8), intent(out) :: a_soa_tend(ntot_amode,ntot_soaspec) ! soa aerosol mixrat tendency (mol/mol/s at host mw)
1138 : ! integer, intent(in) :: idiagss
1139 :
1140 : integer :: ll
1141 : integer :: m,k
1142 :
1143 4348047600 : logical :: skip_soamode(ntot_amode) ! true if this mode does not have soa
1144 :
1145 : real(r8), parameter :: a_min1 = 1.0e-20_r8
1146 : real(r8), parameter :: g_min1 = 1.0e-20_r8
1147 : real(r8), parameter :: alpha = 0.05_r8 ! parameter used in calc of time step
1148 : real(r8), parameter :: dtsub_fixed = -1.0_r8 ! fixed sub-step for time integration (s)
1149 :
1150 4348047600 : real(r8) :: a_ooa_sum_tmp(ntot_soamode) ! total ooa (=soa+opoa) in a mode
1151 4348047600 : real(r8) :: a_opoa(ntot_soamode) ! oxidized-poa aerosol mixrat (mol/mol at actual mw)
1152 4348047600 : real(r8) :: a_soa(ntot_soamode,ntot_soaspec) ! soa aerosol mixrat (mol/mol at actual mw)
1153 4348047600 : real(r8) :: a_soa_tmp(ntot_soamode,ntot_soaspec) ! temporary soa aerosol mixrat (mol/mol)
1154 4348047600 : real(r8) :: beta(ntot_soamode,ntot_soaspec) ! dtcur*xferrate
1155 4348047600 : real(r8) :: delh_vap_soa(ntot_soaspec) ! delh_vap_soa = heat of vaporization for gas soa (J/mol)
1156 4348047600 : real(r8) :: del_g_soa_tmp(ntot_soaspec)
1157 : real(r8) :: dtcur ! current time step (s)
1158 : real(r8) :: dtmax ! = (dtfull-tcur)
1159 4348047600 : real(r8) :: g0_soa(ntot_soaspec) ! ambient soa gas equilib mixrat (mol/mol at actual mw)
1160 4348047600 : real(r8) :: g_soa(ntot_soaspec) ! soa gas mixrat (mol/mol at actual mw)
1161 4348047600 : real(r8) :: g_star(ntot_soamode,ntot_soaspec) ! soa gas mixrat that is in equilib
1162 : ! with each aerosol mode (mol/mol)
1163 4348047600 : real(r8) :: mw_poa(ntot_poaspec) ! actual molec wght of poa
1164 4348047600 : real(r8) :: mw_soa(ntot_soaspec) ! actual molec wght of soa
1165 4348047600 : real(r8) :: opoa_frac(ntot_poaspec) ! fraction of poa that is opoa
1166 4348047600 : real(r8) :: phi(ntot_soamode,ntot_soaspec) ! "relative driving force"
1167 4348047600 : real(r8) :: p0_soa(ntot_soaspec) ! soa gas equilib vapor presssure (atm)
1168 4348047600 : real(r8) :: p0_soa_298(ntot_soaspec) ! p0_soa_298 = soa gas equilib vapor presssure (atm) at 298 k
1169 4348047600 : real(r8) :: sat(ntot_soamode,ntot_soaspec) ! sat(m,ll) = g0_soa(ll)/a_ooa_sum_tmp(m) = g_star(m,ll)/a_soa(m,ll)
1170 : ! used by the numerical integration scheme -- it is not a saturation rato!
1171 : real(r8) :: tcur ! current integration time (from 0 s)
1172 : real(r8) :: tmpa, tmpb, tmpf
1173 4348047600 : real(r8) :: tot_soa(ntot_soaspec) ! g_soa + sum( a_soa(:) )
1174 2174023800 : real(r8) :: xferrate(ntot_amode,ntot_soaspec) ! gas-aerosol mass transfer rate (1/s)
1175 :
1176 : ! Changed by Manish Shrivastava
1177 4348047600 : opoa_frac(:) = 0.0_r8 !POA does not form solution with SOA for all runs; set opoa_frac=0.0_r8 by Manish Shrivastava
1178 4348047600 : mw_poa(:) = 250.0_r8
1179 4348047600 : mw_soa(:) = 250.0_r8
1180 :
1181 : ! New SOA properties added by Manish Shrivastava on 09/27/2012
1182 2174023800 : if (ntot_soaspec ==1) then
1183 4348047600 : p0_soa_298(:) = 9.7831E-11_r8
1184 4348047600 : delh_vap_soa(:) = 131.0e3_r8
1185 4348047600 : opoa_frac(:) = 0.0_r8
1186 0 : elseif (ntot_soaspec ==2) then
1187 : ! same for anthropogenic and biomass burning species
1188 0 : p0_soa_298 (1) = 1.0e-10_r8
1189 0 : p0_soa_298 (2) = 1.0e-10_r8
1190 0 : delh_vap_soa(:) = 156.0e3_r8
1191 0 : elseif(ntot_soaspec ==5) then
1192 : ! 5 volatility bins for each of the a combined SOA classes ( including biomass burning, fossil fuel, biogenic)
1193 0 : p0_soa_298 (1) = 9.7831E-13_r8 !soaff0 C*=0.01ug/m3
1194 0 : p0_soa_298 (2) = 9.7831E-12_r8 !soaff1 C*=0.10ug/m3
1195 0 : p0_soa_298 (3) = 9.7831E-11_r8 !soaff2 C*=1.0ug/m3
1196 0 : p0_soa_298 (4) = 9.7831E-10_r8 !soaff3 C*=10.0ug/m3
1197 0 : p0_soa_298 (5) = 9.7831E-9_r8 !soaff4 C*=100.0ug/m3
1198 :
1199 0 : delh_vap_soa(1) = 153.0e3_r8
1200 0 : delh_vap_soa(2) = 142.0e3_r8
1201 0 : delh_vap_soa(3) = 131.0e3_r8
1202 0 : delh_vap_soa(4) = 120.0e3_r8
1203 0 : delh_vap_soa(5) = 109.0e3_r8
1204 0 : elseif(ntot_soaspec ==15) then
1205 : !
1206 : ! 5 volatility bins for each of the 3 SOA classes ( biomass burning, fossil fuel, biogenic)
1207 : ! SOA species 1-5 are for anthropogenic while 6-10 are for biomass burning SOA
1208 : ! SOA species 11-15 are for biogenic SOA, based on Cappa et al., Reference needs to be updated
1209 : ! For MW=250.0
1210 0 : p0_soa_298 (1) = 9.7831E-13_r8 !soaff0 C*=0.01ug/m3
1211 0 : p0_soa_298 (2) = 9.7831E-12_r8 !soaff1 C*=0.10ug/m3
1212 0 : p0_soa_298 (3) = 9.7831E-11_r8 !soaff2 C*=1.0ug/m3
1213 0 : p0_soa_298 (4) = 9.7831E-10_r8 !soaff3 C*=10.0ug/m3
1214 0 : p0_soa_298 (5) = 9.7831E-9_r8 !soaff4 C*=100.0ug/m3
1215 0 : p0_soa_298 (6) = 9.7831E-13_r8 !soabb0 C*=0.01ug/m3
1216 0 : p0_soa_298 (7) = 9.7831E-12_r8 !soabb1 C*=0.10ug/m3
1217 0 : p0_soa_298 (8) = 9.7831E-11_r8 !soabb2 C*=1.0ug/m3
1218 0 : p0_soa_298 (9) = 9.7831E-10_r8 !soabb3 C*=10.0ug/m3
1219 0 : p0_soa_298 (10) = 9.7831E-9_r8 !soabb4 C*=100.0ug/m3
1220 0 : p0_soa_298 (11) = 9.7831E-13_r8 !soabg0 C*=0.01ug/m3
1221 0 : p0_soa_298 (12) = 9.7831E-12_r8 !soabg1 C*=0.1ug/m3
1222 0 : p0_soa_298 (13) = 9.7831E-11_r8 !soabg2 C*=1.0ug/m3
1223 0 : p0_soa_298 (14) = 9.7831E-10_r8 !soabg3 C*=10.0ug/m3
1224 0 : p0_soa_298 (15) = 9.7831E-9_r8 !soabg4 C*=100.0ug/m3
1225 :
1226 : !
1227 : ! have to be adjusted to 15 species, following the numbers by Epstein et al., 2012
1228 : !
1229 0 : delh_vap_soa(1) = 153.0e3_r8
1230 0 : delh_vap_soa(2) = 142.0e3_r8
1231 0 : delh_vap_soa(3) = 131.0e3_r8
1232 0 : delh_vap_soa(4) = 120.0e3_r8
1233 0 : delh_vap_soa(5) = 109.0e3_r8
1234 0 : delh_vap_soa(6) = 153.0e3_r8
1235 0 : delh_vap_soa(7) = 142.0e3_r8
1236 0 : delh_vap_soa(8) = 131.0e3_r8
1237 0 : delh_vap_soa(9) = 120.0e3_r8
1238 0 : delh_vap_soa(10) = 109.0e3_r8
1239 0 : delh_vap_soa(11) = 153.0e3_r8
1240 0 : delh_vap_soa(12) = 142.0e3_r8
1241 0 : delh_vap_soa(13) = 131.0e3_r8
1242 0 : delh_vap_soa(14) = 120.0e3_r8
1243 0 : delh_vap_soa(15) = 109.0e3_r8
1244 : endif
1245 :
1246 : !BSINGH - Initialized g_soa_tend and a_soa_tend to circumvent the undefined behavior (04/16/12)
1247 4348047600 : g_soa_tend(:) = 0.0_r8
1248 13044142800 : a_soa_tend(:,:) = 0.0_r8
1249 :
1250 : ! determine which modes have non-zero transfer rates
1251 : ! and are involved in the soa gas-aerosol transfer
1252 : ! for diameter = 1 nm and number = 1 #/cm3, xferrate ~= 1e-9 s-1
1253 10870119000 : do m = 1, ntot_soamode
1254 8696095200 : skip_soamode(m) = .true.
1255 19566214200 : do ll = 1, ntot_soaspec
1256 8696095200 : xferrate(m,ll) = xferrate_in(m,ll)
1257 17392190400 : skip_soamode(m) = .false.
1258 : end do
1259 : end do
1260 :
1261 : ! convert incoming mixing ratios from mol/mol at the "host-code" molec. weight (12.0 in cam5)
1262 : ! to mol/mol at the "actual" molec. weight (currently assumed to be 250.0)
1263 : ! also
1264 : ! force things to be non-negative
1265 : ! calc tot_soa(ll)
1266 : ! calc a_opoa (always slightly >0)
1267 4348047600 : do ll = 1, ntot_soaspec
1268 2174023800 : tmpf = mw_soa_host(ll)/mw_soa(ll)
1269 2174023800 : g_soa(ll) = max( g_soa_in(ll), 0.0_r8 ) * tmpf
1270 2174023800 : tot_soa(ll) = g_soa(ll)
1271 13044142800 : do m = 1, ntot_soamode
1272 8696095200 : if ( skip_soamode(m) ) cycle
1273 8696095200 : a_soa(m,ll) = max( a_soa_in(m,ll), 0.0_r8 ) * tmpf
1274 10870119000 : tot_soa(ll) = tot_soa(ll) + a_soa(m,ll)
1275 : end do
1276 : end do
1277 :
1278 10870119000 : do m = 1, ntot_soamode
1279 8696095200 : if ( skip_soamode(m) ) cycle
1280 8696095200 : a_opoa(m) = 0.0_r8
1281 19566214200 : do ll = 1, ntot_poaspec
1282 17392190400 : a_opoa(m) = opoa_frac(ll)*a_poa_in(m,ll)
1283 : end do
1284 : end do
1285 :
1286 : ! calc ambient equilibrium soa gas
1287 4348047600 : do ll = 1, ntot_soaspec
1288 2174023800 : p0_soa(ll) = p0_soa_298(ll) * &
1289 2174023800 : exp( -(delh_vap_soa(ll)/rgas)*((1.0_r8/temp)-(1.0_r8/298.0_r8)) )
1290 4348047600 : g0_soa(ll) = 1.01325e5_r8*p0_soa(ll)/pres
1291 : end do
1292 :
1293 2174023800 : niter = 0
1294 2174023800 : tcur = 0.0_r8
1295 2174023800 : dtcur = 0.0_r8
1296 13044142800 : phi(:,:) = 0.0_r8
1297 13044142800 : g_star(:,:) = 0.0_r8
1298 :
1299 : ! if (idiagss > 0) then
1300 : ! write(luna,'(a,1p,10e11.3)') 'p0, g0_soa', p0_soa, g0_soa
1301 : ! write(luna,'(3a)') &
1302 : ! 'niter, tcur, dtcur, phi(:), ', &
1303 : ! 'g_star(:), ', &
1304 : ! 'a_soa(:), g_soa'
1305 : ! write(luna,'(3a)') &
1306 : ! ' sat(:), ', &
1307 : ! 'sat(:)*a_soa(:) ', &
1308 : ! 'a_opoa(:)'
1309 : ! write(luna,'(i3,1p,20e10.2)') niter, tcur, dtcur, &
1310 : ! phi(:), g_star(:), a_soa(:), g_soa
1311 : ! end if
1312 :
1313 :
1314 : ! integration loop -- does multiple substeps to reach dtfull
1315 : time_loop: &
1316 8438889498 : do while (tcur < dtfull-1.0e-3_r8 )
1317 :
1318 6264865698 : niter = niter + 1
1319 6264865698 : if (niter > niter_max) exit
1320 :
1321 31324328490 : tmpa = 0.0_r8 ! time integration parameter for all soa species
1322 31324328490 : do m = 1, ntot_soamode
1323 25059462792 : if ( skip_soamode(m) ) cycle
1324 56383791282 : a_ooa_sum_tmp(m) = a_opoa(m) + sum( a_soa(m,1:ntot_soaspec) )
1325 : end do
1326 12529731396 : do ll = 1, ntot_soaspec
1327 : tmpb = 0.0_r8 ! time integration parameter for a single soa species
1328 31324328490 : do m = 1, ntot_soamode
1329 25059462792 : if ( skip_soamode(m) ) cycle
1330 25059462792 : sat(m,ll) = g0_soa(ll)/max( a_ooa_sum_tmp(m), a_min1 )
1331 25059462792 : g_star(m,ll) = sat(m,ll)*a_soa(m,ll)
1332 25059462792 : phi(m,ll) = (g_soa(ll) - g_star(m,ll))/max( g_soa(ll), g_star(m,ll), g_min1 )
1333 31324328490 : tmpb = tmpb + xferrate(m,ll)*abs(phi(m,ll))
1334 : end do
1335 12529731396 : tmpa = max( tmpa, tmpb )
1336 : end do
1337 :
1338 : if (dtsub_fixed > 0.0_r8) then
1339 : dtcur = dtsub_fixed
1340 : tcur = tcur + dtcur
1341 : else
1342 6264865698 : dtmax = dtfull-tcur
1343 6264865698 : if (dtmax*tmpa <= alpha) then
1344 : ! here alpha/tmpa >= dtmax, so this is final substep
1345 : dtcur = dtmax
1346 : tcur = dtfull
1347 : else
1348 4090843651 : dtcur = alpha/tmpa
1349 4090843651 : tcur = tcur + dtcur
1350 : end if
1351 : end if
1352 :
1353 : ! step 1 - for modes where soa is condensing, estimate "new" a_soa(m,ll)
1354 : ! using an explicit calculation with "old" g_soa
1355 : ! and g_star(m,ll) calculated using "old" a_soa(m,ll)
1356 : ! do this to get better estimate of "new" a_soa(m,ll) and sat(m,ll)
1357 31324328490 : do m = 1, ntot_soamode
1358 25059462792 : if ( skip_soamode(m) ) cycle
1359 50118925584 : do ll = 1, ntot_soaspec
1360 : ! first ll loop calcs a_soa_tmp(m,ll) & a_ooa_sum_tmp
1361 25059462792 : a_soa_tmp(m,ll) = a_soa(m,ll)
1362 25059462792 : beta(m,ll) = dtcur*xferrate(m,ll)
1363 25059462792 : del_g_soa_tmp(ll) = g_soa(ll) - g_star(m,ll)
1364 50118925584 : if (del_g_soa_tmp(ll) > 0.0_r8) then
1365 19816107746 : a_soa_tmp(m,ll) = a_soa(m,ll) + beta(m,ll)*del_g_soa_tmp(ll)
1366 : end if
1367 : end do
1368 50118925584 : a_ooa_sum_tmp(m) = a_opoa(m) + sum( a_soa_tmp(m,1:ntot_soaspec) )
1369 56383791282 : do ll = 1, ntot_soaspec
1370 : ! second ll loop calcs sat & g_star
1371 50118925584 : if (del_g_soa_tmp(ll) > 0.0_r8) then
1372 19816107746 : sat(m,ll) = g0_soa(ll)/max( a_ooa_sum_tmp(m), a_min1 )
1373 19816107746 : g_star(m,ll) = sat(m,ll)*a_soa_tmp(m,ll) ! this just needed for diagnostics
1374 : end if
1375 : end do
1376 : end do
1377 :
1378 : ! step 2 - implicit in g_soa and semi-implicit in a_soa,
1379 : ! with g_star(m,ll) calculated semi-implicitly
1380 14703755196 : do ll = 1, ntot_soaspec
1381 : tmpa = 0.0_r8
1382 : tmpb = 0.0_r8
1383 31324328490 : do m = 1, ntot_soamode
1384 25059462792 : if ( skip_soamode(m) ) cycle
1385 25059462792 : tmpa = tmpa + a_soa(m,ll)/(1.0_r8 + beta(m,ll)*sat(m,ll))
1386 31324328490 : tmpb = tmpb + beta(m,ll)/(1.0_r8 + beta(m,ll)*sat(m,ll))
1387 : end do
1388 :
1389 6264865698 : g_soa(ll) = (tot_soa(ll) - tmpa)/(1.0_r8 + tmpb)
1390 6264865698 : g_soa(ll) = max( 0.0_r8, g_soa(ll) )
1391 37589194188 : do m = 1, ntot_soamode
1392 25059462792 : if ( skip_soamode(m) ) cycle
1393 25059462792 : a_soa(m,ll) = (a_soa(m,ll) + beta(m,ll)*g_soa(ll))/ &
1394 31324328490 : (1.0_r8 + beta(m,ll)*sat(m,ll))
1395 : end do
1396 : end do
1397 :
1398 : ! if (idiagss > 0) then
1399 : ! write(luna,'(i3,1p,20e10.2)') niter, tcur, dtcur, &
1400 : ! phi(:), g_star(:), a_soa(:), g_soa
1401 : ! write(luna,'(23x,1p,20e10.2)') &
1402 : ! sat(:), sat(:)*a_soa(:), a_opoa(:)
1403 : ! end if
1404 :
1405 : ! if (niter > 9992000) then
1406 : ! write(luna,'(a)') '*** to many iterations'
1407 : ! exit
1408 : ! end if
1409 :
1410 : end do time_loop
1411 :
1412 :
1413 : ! calculate outgoing tendencies (at the host-code molec. weight)
1414 : ! (a_soa & g_soa are at actual mw, but a_soa_in & g_soa_in are at host-code mw)
1415 4348047600 : do ll = 1, ntot_soaspec
1416 2174023800 : tmpf = mw_soa(ll)/mw_soa_host(ll)
1417 2174023800 : g_soa_tend(ll) = (g_soa(ll)*tmpf - g_soa_in(ll))/dtfull
1418 13044142800 : do m = 1, ntot_soamode
1419 8696095200 : if ( skip_soamode(m) ) cycle
1420 10870119000 : a_soa_tend(m,ll) = (a_soa(m,ll)*tmpf - a_soa_in(m,ll))/dtfull
1421 : end do
1422 : end do
1423 :
1424 :
1425 2174023800 : return
1426 :
1427 : end subroutine modal_aero_soaexch
1428 :
1429 : !----------------------------------------------------------------------
1430 :
1431 : !----------------------------------------------------------------------
1432 :
1433 1536 : subroutine modal_aero_gasaerexch_init
1434 :
1435 : !-----------------------------------------------------------------------
1436 : !
1437 : ! Purpose:
1438 : ! set do_adjust and do_aitken flags
1439 : ! create history fields for column tendencies associated with
1440 : ! modal_aero_calcsize
1441 : !
1442 : ! Author: R. Easter
1443 : !
1444 : !-----------------------------------------------------------------------
1445 :
1446 : use modal_aero_data
1447 : use modal_aero_rename
1448 :
1449 : use cam_abortutils, only: endrun
1450 : use cam_history, only: addfld, add_default, fieldname_len, horiz_only
1451 : use constituents, only: pcnst, cnst_get_ind, cnst_name
1452 : use spmd_utils, only: masterproc
1453 : use phys_control, only: phys_getopts
1454 :
1455 : implicit none
1456 :
1457 : !-----------------------------------------------------------------------
1458 : ! arguments
1459 :
1460 : !-----------------------------------------------------------------------
1461 : ! local
1462 : integer :: ipair, iq, iqfrm, iqfrm_aa, iqtoo, iqtoo_aa
1463 : integer :: jac,jsoa,p
1464 : integer :: l, l1, l2, lsfrm, lstoo, lunout
1465 : integer :: l_so4g, l_nh4g, l_msag
1466 : integer :: m, mfrm, mtoo
1467 : integer :: n, nacc, nait
1468 : integer :: nchfrm, nchfrmskip, nchtoo, nchtooskip, nspec
1469 :
1470 : logical :: do_msag, do_nh4g
1471 3072 : logical :: do_soag_any, do_soag(nsoa)
1472 : logical :: dotend(pcnst), dotendqqcw(pcnst)
1473 :
1474 : real(r8) :: tmp1, tmp2
1475 :
1476 : character(len=fieldname_len) :: tmpnamea, tmpnameb
1477 : character(len=fieldname_len+3) :: fieldname
1478 : character(128) :: long_name
1479 : character(128) :: msg
1480 : character(8) :: unit
1481 :
1482 : logical :: history_aerosol ! Output the MAM aerosol tendencies
1483 : logical :: history_aerocom ! Output the aerocom history
1484 : !-----------------------------------------------------------------------
1485 :
1486 1536 : call phys_getopts( history_aerosol_out = history_aerosol )
1487 :
1488 1536 : maxspec_pcage = nspec_max
1489 4608 : allocate(lspecfrm_pcage(maxspec_pcage))
1490 3072 : allocate(lspectoo_pcage(maxspec_pcage))
1491 4608 : allocate(soa_equivso4_factor(nsoa))
1492 3072 : allocate(fac_m2v_soa(nsoa))
1493 4608 : allocate(fac_m2v_pcarbon(nspec_max))
1494 1536 : lunout = 6
1495 : !
1496 : ! define "from mode" and "to mode" for primary carbon aging
1497 : !
1498 : ! skip (turn off) aging if either is absent,
1499 : ! or if accum mode so4 is absent
1500 : !
1501 1536 : modefrm_pcage = -999888777
1502 1536 : modetoo_pcage = -999888777
1503 1536 : if ((modeptr_pcarbon <= 0) .or. (modeptr_accum <= 0)) goto 15000
1504 1536 : l = lptr_so4_a_amode(modeptr_accum)
1505 1536 : if ((l < 1) .or. (l > pcnst)) goto 15000
1506 :
1507 1536 : modefrm_pcage = modeptr_pcarbon
1508 1536 : modetoo_pcage = modeptr_accum
1509 :
1510 : !
1511 : ! define species involved in each primary carbon aging pairing
1512 : ! (include aerosol water)
1513 : !
1514 : !
1515 1536 : mfrm = modefrm_pcage
1516 1536 : mtoo = modetoo_pcage
1517 :
1518 1536 : if (mfrm < 10) then
1519 : nchfrmskip = 1
1520 0 : else if (mfrm < 100) then
1521 : nchfrmskip = 2
1522 : else
1523 0 : nchfrmskip = 3
1524 : end if
1525 1536 : if (mtoo < 10) then
1526 : nchtooskip = 1
1527 0 : else if (mtoo < 100) then
1528 : nchtooskip = 2
1529 : else
1530 0 : nchtooskip = 3
1531 : end if
1532 1536 : nspec = 0
1533 :
1534 7680 : aa_iqfrm: do iqfrm = -1, nspec_amode(mfrm)
1535 :
1536 6144 : if (iqfrm == -1) then
1537 1536 : lsfrm = numptr_amode(mfrm)
1538 1536 : lstoo = numptr_amode(mtoo)
1539 4608 : else if (iqfrm == 0) then
1540 : ! bypass transfer of aerosol water due to primary-carbon aging
1541 : cycle aa_iqfrm
1542 : ! lsfrm = lwaterptr_amode(mfrm)
1543 : ! lstoo = lwaterptr_amode(mtoo)
1544 : else
1545 3072 : lsfrm = lmassptr_amode(iqfrm,mfrm)
1546 3072 : lstoo = 0
1547 : end if
1548 4608 : if ((lsfrm < 1) .or. (lsfrm > pcnst)) cycle aa_iqfrm
1549 :
1550 4608 : if (lsfrm>0 .and. iqfrm>0 ) then
1551 3072 : nchfrm = len( trim( cnst_name(lsfrm) ) ) - nchfrmskip
1552 :
1553 : ! find "too" species having same lspectype_amode as the "frm" species
1554 : ! AND same cnst_name (except for last 1/2/3 characters which are the mode index)
1555 9216 : do iqtoo = 1, nspec_amode(mtoo)
1556 : ! if ( lspectype_amode(iqtoo,mtoo) .eq. &
1557 : ! lspectype_amode(iqfrm,mfrm) ) then
1558 9216 : lstoo = lmassptr_amode(iqtoo,mtoo)
1559 9216 : nchtoo = len( trim( cnst_name(lstoo) ) ) - nchtooskip
1560 9216 : if (cnst_name(lsfrm)(1:nchfrm) == cnst_name(lstoo)(1:nchtoo)) then
1561 : exit
1562 : else
1563 6144 : lstoo = 0
1564 : end if
1565 : ! end if
1566 : end do
1567 : end if
1568 :
1569 4608 : if ((lstoo < 1) .or. (lstoo > pcnst)) lstoo = 0
1570 4608 : nspec = nspec + 1
1571 4608 : lspecfrm_pcage(nspec) = lsfrm
1572 7680 : lspectoo_pcage(nspec) = lstoo
1573 : end do aa_iqfrm
1574 :
1575 1536 : nspecfrm_pcage = nspec
1576 :
1577 : !
1578 : ! output results
1579 : !
1580 1536 : if ( masterproc ) then
1581 :
1582 2 : write(lunout,9310)
1583 :
1584 2 : mfrm = modefrm_pcage
1585 2 : mtoo = modetoo_pcage
1586 2 : write(lunout,9320) 1, mfrm, mtoo
1587 :
1588 8 : do iq = 1, nspecfrm_pcage
1589 6 : lsfrm = lspecfrm_pcage(iq)
1590 6 : lstoo = lspectoo_pcage(iq)
1591 8 : if (lstoo .gt. 0) then
1592 6 : write(lunout,9330) lsfrm, cnst_name(lsfrm), &
1593 12 : lstoo, cnst_name(lstoo)
1594 : else
1595 0 : write(lunout,9340) lsfrm, cnst_name(lsfrm)
1596 : end if
1597 : end do
1598 :
1599 2 : write(lunout,*)
1600 :
1601 : end if ! ( masterproc )
1602 :
1603 : 9310 format( / 'subr. modal_aero_gasaerexch_init - primary carbon aging pointers' )
1604 : 9320 format( 'pair', i3, 5x, 'mode', i3, ' ---> mode', i3 )
1605 : 9330 format( 5x, 'spec', i3, '=', a, ' ---> spec', i3, '=', a )
1606 : 9340 format( 5x, 'spec', i3, '=', a, ' ---> LOSS' )
1607 :
1608 :
1609 : 15000 continue
1610 :
1611 : ! set tendency flags and gas species indices and flags
1612 1536 : dotend(:) = .false.
1613 :
1614 1536 : call cnst_get_ind( 'H2SO4', l_so4g, .false. )
1615 1536 : if ((l_so4g <= 0) .or. (l_so4g > pcnst)) then
1616 : write( *, '(/a/a,2i7)' ) &
1617 0 : '*** modal_aero_gasaerexch_init -- cannot find H2SO4 species', &
1618 0 : ' l_so4g=', l_so4g
1619 0 : call endrun( 'modal_aero_gasaerexch_init error' )
1620 : end if
1621 1536 : dotend(l_so4g) = .true.
1622 :
1623 1536 : call cnst_get_ind( 'NH3', l_nh4g, .false. )
1624 1536 : do_nh4g = .false.
1625 1536 : if ((l_nh4g > 0) .and. (l_nh4g <= pcnst)) then
1626 0 : do_nh4g = .true.
1627 0 : dotend(l_nh4g) = .true.
1628 : end if
1629 :
1630 1536 : call cnst_get_ind( 'MSA', l_msag, .false. )
1631 1536 : do_msag = .false.
1632 1536 : if ((l_msag > 0) .and. (l_msag <= pcnst)) then
1633 0 : do_msag = .true.
1634 0 : dotend(l_msag) = .true.
1635 : end if
1636 :
1637 1536 : do_soag_any = .false.
1638 3072 : do_soag(:) = .false.
1639 3072 : do jsoa = 1, nsoa
1640 1536 : l = lptr2_soa_g_amode(jsoa)
1641 3072 : if ((l > 0) .and. (l <= pcnst)) then
1642 1536 : do_soag_any = .true.
1643 1536 : do_soag(jsoa) = .true.
1644 1536 : dotend(l) = .true.
1645 : end if
1646 : end do
1647 :
1648 :
1649 7680 : do n = 1, ntot_amode
1650 6144 : l = lptr_so4_a_amode(n)
1651 6144 : if ((l > 0) .and. (l <= pcnst)) then
1652 4608 : dotend(l) = .true.
1653 4608 : if ( do_nh4g ) then
1654 0 : l = lptr_nh4_a_amode(n)
1655 0 : if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
1656 : end if
1657 : end if
1658 13824 : do jsoa = 1, nsoa
1659 12288 : if ( do_soag(jsoa) ) then
1660 6144 : l = lptr2_soa_a_amode(n,jsoa)
1661 6144 : if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
1662 : end if
1663 : end do
1664 : end do
1665 :
1666 1536 : if (modefrm_pcage > 0) then
1667 6144 : do iq = 1, nspecfrm_pcage
1668 4608 : lsfrm = lspecfrm_pcage(iq)
1669 4608 : lstoo = lspectoo_pcage(iq)
1670 6144 : if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
1671 4608 : dotend(lsfrm) = .true.
1672 4608 : if ((lstoo > 0) .and. (lstoo <= pcnst)) then
1673 4608 : dotend(lstoo) = .true.
1674 : end if
1675 : end if
1676 : end do
1677 : end if
1678 :
1679 : !---------define history fields for new cond/evap diagnostics----------------------------------------
1680 1536 : fieldname=trim('qconff_gaex')
1681 1536 : long_name = trim('3D fields for Fossil SOA condensation')
1682 1536 : unit = 'kg/kg/s'
1683 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
1684 1536 : if ( history_aerosol ) then
1685 0 : call add_default( fieldname, 1, ' ' )
1686 : endif
1687 1536 : if ( masterproc ) write(*,'(3(a,3x))') 'qconff addfld', fieldname, unit
1688 :
1689 1536 : fieldname=trim('qevapff_gaex')
1690 1536 : long_name = trim('3D fields for Fossil SOA evaporation')
1691 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
1692 1536 : if ( history_aerosol ) then
1693 0 : call add_default( fieldname, 1, ' ' )
1694 : endif
1695 1536 : if ( masterproc ) write(*,'(3(a,3x))') 'qevapff addfld', fieldname, unit
1696 :
1697 1536 : fieldname=trim('qconbb_gaex')
1698 1536 : long_name = trim('3D fields for Biomass SOA condensation')
1699 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
1700 1536 : if ( history_aerosol ) then
1701 0 : call add_default( fieldname, 1, ' ' )
1702 : endif
1703 1536 : if ( masterproc ) write(*,'(3(a,3x))') 'qconbb addfld', fieldname, unit
1704 :
1705 1536 : fieldname=trim('qevapbb_gaex')
1706 1536 : long_name = trim('3D fields for Biomass SOA evaporation')
1707 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
1708 1536 : if ( history_aerosol ) then
1709 0 : call add_default( fieldname, 1, ' ' )
1710 : endif
1711 1536 : if ( masterproc ) write(*,'(3(a,3x))') 'qevapbb addfld', fieldname, unit
1712 :
1713 1536 : fieldname=trim('qconbg_gaex')
1714 1536 : long_name = trim('3D fields for Biogenic SOA condensation')
1715 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
1716 1536 : if ( history_aerosol ) then
1717 0 : call add_default( fieldname, 1, ' ' )
1718 : endif
1719 1536 : if ( masterproc ) write(*,'(3(a,3x))') 'qconbg addfld', fieldname, unit
1720 :
1721 1536 : fieldname=trim('qevapbg_gaex')
1722 1536 : long_name = trim('3D fields for Biogenic SOA evaporation')
1723 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
1724 1536 : if ( history_aerosol ) then
1725 0 : call add_default( fieldname, 1, ' ' )
1726 : endif
1727 1536 : if ( masterproc ) write(*,'(3(a,3x))') 'qevapbg addfld', fieldname, unit
1728 :
1729 1536 : fieldname=trim('qcon_gaex')
1730 1536 : long_name = trim('3D fields for SOA condensation')
1731 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
1732 1536 : if ( history_aerosol ) then
1733 0 : call add_default( fieldname, 1, ' ' )
1734 : endif
1735 1536 : if ( masterproc ) write(*,'(3(a,3x))') 'qcon addfld', fieldname, unit
1736 :
1737 1536 : fieldname=trim('qevap_gaex')
1738 1536 : long_name = trim('3D fields for Biogenic SOA evaporation')
1739 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
1740 1536 : if ( history_aerosol ) then
1741 0 : call add_default( fieldname, 1, ' ' )
1742 : endif
1743 1536 : if ( masterproc ) write(*,'(3(a,3x))') 'qevap addfld', fieldname, unit
1744 : !------------------------------------------------------------------------------
1745 :
1746 : ! define history fields for basic gas-aer exchange
1747 : ! and primary carbon aging from that
1748 64512 : do l = 1, pcnst
1749 62976 : if ( .not. dotend(l) ) cycle
1750 :
1751 19968 : tmpnamea = cnst_name(l)
1752 19968 : fieldname = trim(tmpnamea) // '_sfgaex1'
1753 19968 : long_name = trim(tmpnamea) // ' gas-aerosol-exchange primary column tendency'
1754 19968 : unit = 'kg/m2/s'
1755 19968 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
1756 19968 : if ( history_aerosol ) then
1757 0 : call add_default( fieldname, 1, ' ' )
1758 : endif
1759 21504 : if ( masterproc ) write(*,'(3(a,3x))') 'gasaerexch addfld', fieldname, unit
1760 :
1761 : end do ! l = ...
1762 : ! define history fields for aitken-->accum renaming
1763 1536 : dotend(:) = .false.
1764 1536 : dotendqqcw(:) = .false.
1765 6144 : do ipair = 1, npair_renamexf
1766 26112 : do iq = 1, nspecfrm_renamexf(ipair)
1767 19968 : lsfrm = lspecfrma_renamexf(iq,ipair)
1768 19968 : lstoo = lspectooa_renamexf(iq,ipair)
1769 19968 : if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
1770 19968 : dotend(lsfrm) = .true.
1771 19968 : if ((lstoo > 0) .and. (lstoo <= pcnst)) then
1772 19968 : dotend(lstoo) = .true.
1773 : end if
1774 : end if
1775 :
1776 19968 : lsfrm = lspecfrmc_renamexf(iq,ipair)
1777 19968 : lstoo = lspectooc_renamexf(iq,ipair)
1778 24576 : if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
1779 19968 : dotendqqcw(lsfrm) = .true.
1780 19968 : if ((lstoo > 0) .and. (lstoo <= pcnst)) then
1781 19968 : dotendqqcw(lstoo) = .true.
1782 : end if
1783 : end if
1784 : end do ! iq = ...
1785 : end do ! ipair = ...
1786 :
1787 64512 : do l = 1, pcnst
1788 190464 : do jac = 1, 2
1789 125952 : if (jac == 1) then
1790 62976 : if ( .not. dotend(l) ) cycle
1791 21504 : tmpnamea = cnst_name(l)
1792 : else
1793 62976 : if ( .not. dotendqqcw(l) ) cycle
1794 21504 : tmpnamea = cnst_name_cw(l)
1795 : end if
1796 :
1797 43008 : fieldname = trim(tmpnamea) // '_sfgaex2'
1798 43008 : long_name = trim(tmpnamea) // ' gas-aerosol-exchange renaming column tendency'
1799 43008 : unit = 'kg/m2/s'
1800 43008 : if ((tmpnamea(1:3) == 'num') .or. &
1801 9216 : (tmpnamea(1:3) == 'NUM')) unit = '#/m2/s'
1802 43008 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
1803 43008 : if ( history_aerosol ) then
1804 0 : call add_default( fieldname, 1, ' ' )
1805 : endif
1806 105984 : if ( masterproc ) write(*,'(3(a,3x))') 'gasaerexch addfld', fieldname, unit
1807 : end do ! jac = ...
1808 : end do ! l = ...
1809 :
1810 :
1811 : ! set for used in aging calcs:
1812 : ! fac_m2v_so4, fac_m2v_nh4, fac_m2v_soa(:)
1813 : ! soa_equivso4_factor(:)
1814 3072 : soa_equivso4_factor = 0.0_r8
1815 1536 : if (modefrm_pcage > 0) then
1816 1536 : n = modeptr_accum
1817 1536 : l = lptr_so4_a_amode(n) ; l2 = -1
1818 1536 : if (l <= 0) call endrun( 'modal_aero_gasaerexch_init error a001 finding accum. so4' )
1819 10752 : do l1 = 1, nspec_amode(n)
1820 10752 : if (lmassptr_amode(l1,n) == l) then
1821 : ! l2 = lspectype_amode(l1,n)
1822 1536 : l2 = l1
1823 : ! fac_m2v_so4 = specmw_amode(l2) / specdens_amode(l2)
1824 1536 : fac_m2v_so4 = specmw_amode(l1,n) / specdens_amode(l1,n)
1825 : ! tmp2 = spechygro(l2)
1826 1536 : tmp2 = spechygro(l1,n)
1827 :
1828 : end if
1829 : end do
1830 1536 : if (l2 <= 0) call endrun( 'modal_aero_gasaerexch_init error a002 finding accum. so4' )
1831 :
1832 1536 : l = lptr_nh4_a_amode(n) ; l2 = -1
1833 1536 : if (l > 0) then
1834 0 : do l1 = 1, nspec_amode(n)
1835 0 : if (lmassptr_amode(l1,n) == l) then
1836 : ! l2 = lspectype_amode(l1,n)
1837 0 : l2 = l1
1838 : ! fac_m2v_nh4 = specmw_amode(l2) / specdens_amode(l2)
1839 0 : fac_m2v_nh4 = specmw_amode(l1,n) / specdens_amode(l1,n)
1840 :
1841 : end if
1842 : end do
1843 0 : if (l2 <= 0) call endrun( 'modal_aero_gasaerexch_init error a002 finding accum. nh4' )
1844 : else
1845 1536 : fac_m2v_nh4 = fac_m2v_so4
1846 : end if
1847 :
1848 3072 : do jsoa = 1, nsoa
1849 1536 : l = lptr2_soa_a_amode(n,jsoa) ; l2 = -1
1850 1536 : if (l <= 0) then
1851 0 : write( msg, '(a,i4)') 'modal_aero_gasaerexch_init error a001 finding accum. jsoa =', jsoa
1852 0 : call endrun( msg )
1853 : end if
1854 10752 : do l1 = 1, nspec_amode(n)
1855 10752 : if (lmassptr_amode(l1,n) == l) then
1856 : ! l2 = lspectype_amode(l1,n)
1857 1536 : l2 = l1
1858 : ! fac_m2v_soa(jsoa) = specmw_amode(l2) / specdens_amode(l2)
1859 1536 : fac_m2v_soa(jsoa) = specmw_amode(l1,n) / specdens_amode(l1,n)
1860 : ! soa_equivso4_factor(jsoa) = spechygro(l2)/tmp2
1861 1536 : soa_equivso4_factor(jsoa) = spechygro(l1,n)/tmp2
1862 : end if
1863 : end do
1864 3072 : if (l2 <= 0) then
1865 0 : write( msg, '(a,i4)') 'modal_aero_gasaerexch_init error a002 finding accum. jsoa =', jsoa
1866 0 : call endrun( msg )
1867 : end if
1868 : end do
1869 :
1870 10752 : fac_m2v_pcarbon(:) = 0.0_r8
1871 1536 : n = modeptr_pcarbon
1872 4608 : do l = 1, nspec_amode(n)
1873 : ! l2 = lspectype_amode(l,n)
1874 : ! fac_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
1875 : ! [m3-AP/kmol-AP] = [kg-AP/kmol-AP] / [kg-AP/m3-AP]
1876 : ! fac_m2v_pcarbon(l) = specmw_amode(l2) / specdens_amode(l2)
1877 4608 : fac_m2v_pcarbon(l) = specmw_amode(l,n) / specdens_amode(l,n)
1878 : end do
1879 : end if
1880 :
1881 :
1882 1536 : return
1883 :
1884 3072 : end subroutine modal_aero_gasaerexch_init
1885 :
1886 :
1887 : !----------------------------------------------------------------------
1888 :
1889 : end module modal_aero_gasaerexch
|
