Line data Source code
1 : !----------------------------------------------------------------------
2 : !----------------------------------------------------------------------
3 : !BOP
4 : !
5 : ! !MODULE: carma_aero_gasaerexch --- does carma aerosol gas-aerosol exchange for SOA
6 : !
7 : ! !INTERFACE:
8 : module carma_aero_gasaerexch
9 :
10 : ! !USES:
11 : use shr_kind_mod, only: r8 => shr_kind_r8
12 : use chem_mods, only: gas_pcnst
13 : use ref_pres, only: top_lev => clim_modal_aero_top_lev
14 : use ppgrid, only: pcols, pver
15 : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_info_by_bin, rad_cnst_get_bin_props_by_idx, &
16 : rad_cnst_get_info_by_bin_spec
17 : use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field
18 :
19 : implicit none
20 : private
21 : public :: carma_aero_gasaerexch_sub
22 : public :: carma_aero_gasaerexch_init
23 :
24 : !PUBLIC DATA MEMBERS:
25 :
26 : ! description of bin aerosols
27 : integer, public, protected :: nspec_max = 0
28 : integer, public, protected :: nbins = 0
29 : integer, public, protected :: nsoa_vbs = 0
30 : integer, public, protected :: nsoa = 0
31 : integer, public, protected :: npoa = 0
32 : integer, public, protected, allocatable :: nspec(:)
33 :
34 : ! Misc private data
35 : character(len=32), allocatable :: fldname(:) ! names for interstitial output fields
36 : character(len=32), allocatable :: fldname_cw(:) ! names for cloud_borne output fields
37 :
38 : ! local indexing for bins
39 : integer, allocatable :: bin_idx(:,:) ! table for local indexing of modal aero number and mmr
40 : integer :: ncnst_tot ! total number of mode number conc + mode species
41 :
42 : real(r8) :: mw_soa = 250._r8
43 : integer :: fracvbs_idx = -1
44 : integer, allocatable :: dqdtsoa_idx(:,:)
45 : integer, allocatable :: cnsoa(:) ! true if soa gas is a species and carma soa in bin
46 : integer, allocatable :: cnpoa(:) ! true if soa gas is a species and carma soa in bin
47 : integer, allocatable :: l_soag(:) ! true if soa gas is a species and carma soa in bin
48 :
49 : logical, allocatable :: do_soag_any(:) ! true if soa gas is a species and carma soa in bin
50 : ! !DESCRIPTION: This module implements ...
51 : !
52 : ! !REVISION HISTORY:
53 : !
54 : ! RCE 07.04.13: Adapted from MIRAGE2 code
55 : !
56 : !EOP
57 : !----------------------------------------------------------------------
58 : !BOC
59 :
60 : ! list private module data here
61 :
62 : !EOC
63 : !----------------------------------------------------------------------
64 :
65 :
66 : contains
67 :
68 : !----------------------------------------------------------------------
69 :
70 1536 : subroutine carma_aero_gasaerexch_init
71 :
72 : !-----------------------------------------------------------------------
73 : !
74 : ! Purpose:
75 : ! gas-aerosol exchange SOAG <-> soa
76 : !
77 : ! Author: Simone Tilmes
78 : !
79 : !-----------------------------------------------------------------------
80 :
81 : use cam_history, only: addfld, add_default, fieldname_len, horiz_only
82 : use constituents, only: pcnst, cnst_name
83 : use phys_control, only: phys_getopts
84 : use mo_chem_utls, only: get_spc_ndx
85 :
86 : !-----------------------------------------------------------------------
87 : ! arguments
88 :
89 : !-----------------------------------------------------------------------
90 : ! local
91 : integer :: j
92 : integer :: i, ii
93 : integer :: l
94 : integer :: m
95 : integer :: ns
96 : character(len=fieldname_len+3) :: fieldname
97 : character(len=32) :: spectype
98 : character(len=32) :: spec_name
99 : character(128) :: long_name
100 : character(8) :: unit
101 : character(len=2) :: outsoa
102 :
103 : logical :: history_aerosol ! Output the MAM aerosol tendencies
104 : !-----------------------------------------------------------------------
105 :
106 1536 : call phys_getopts( history_aerosol_out = history_aerosol )
107 :
108 : !
109 : ! get info about the bin aerosols
110 : ! get nbins
111 :
112 1536 : call rad_cnst_get_info( 0, nbins=nbins)
113 :
114 4608 : allocate( nspec(nbins) )
115 3072 : allocate( cnsoa(nbins) )
116 3072 : allocate( cnpoa(nbins) )
117 :
118 62976 : do m = 1, nbins
119 62976 : call rad_cnst_get_info_by_bin(0, m, nspec=nspec(m))
120 : end do
121 :
122 62976 : nspec_max = maxval(nspec)
123 :
124 1536 : ncnst_tot = nspec(1)
125 61440 : do m = 2, nbins
126 61440 : ncnst_tot = ncnst_tot + nspec(m)
127 : end do
128 :
129 0 : allocate( bin_idx(nbins,nspec_max), &
130 0 : do_soag_any(nbins), &
131 0 : fldname_cw(ncnst_tot), &
132 13824 : fldname(ncnst_tot) )
133 :
134 : ! Local indexing compresses the mode and number/mass indicies into one index.
135 : ! This indexing is used by the pointer arrays used to reference state and pbuf
136 : ! fields.
137 : ! for CARMA we add number = 0, total mass = 1, and mass from each constituence into mm.
138 1536 : ii = 0
139 62976 : do m = 1, nbins
140 278016 : do l = 1, nspec(m) ! do through nspec
141 215040 : ii = ii + 1
142 276480 : bin_idx(m,l) = ii
143 : end do
144 : end do
145 :
146 : ! SOAG / SOA / POM information
147 : ! Define number of VBS bins (nsoa) based on number of SOAG chemistry species
148 :
149 1536 : nsoa_vbs = 0
150 336384 : do i = 1, pcnst
151 336384 : if (cnst_name(i)(:4) == 'SOAG') then
152 1536 : nsoa_vbs = nsoa_vbs + 1
153 : end if
154 : end do
155 4608 : allocate( l_soag(nsoa_vbs) )
156 1536 : nsoa_vbs = 0
157 336384 : do i = 1, pcnst
158 336384 : if (cnst_name(i)(:4) == 'SOAG') then
159 1536 : nsoa_vbs = nsoa_vbs + 1
160 1536 : l_soag(nsoa_vbs) = get_spc_ndx(cnst_name(i))
161 : end if
162 : end do
163 :
164 1536 : fracvbs_idx = pbuf_get_index('FRACVBS')
165 :
166 : ! identify number of SOA and POA in CARMA code (CARMA number cn)
167 62976 : do m = 1, nbins
168 61440 : cnsoa(m) = 0
169 61440 : cnpoa(m) = 0
170 278016 : do l = 1, nspec(m)
171 215040 : call rad_cnst_get_bin_props_by_idx(0, m, l,spectype=spectype)
172 215040 : if (trim(spectype) == 's-organic') then
173 30720 : cnsoa(m) = cnsoa(m) + 1
174 : end if
175 276480 : if (trim(spectype) == 'p-organic') then
176 30720 : cnpoa(m) = cnpoa(m) + 1
177 : end if
178 : end do
179 : end do
180 : ! some bins don't contain soa or poa
181 62976 : nsoa= maxval(cnsoa)
182 62976 : npoa= maxval(cnpoa)
183 :
184 6144 : allocate( dqdtsoa_idx(nbins,nsoa) )
185 62976 : do m = 1, nbins
186 61440 : ns = 0
187 278016 : do l = 1, nspec(m)
188 215040 : call rad_cnst_get_bin_props_by_idx(0, m, l,spectype=spectype)
189 276480 : if (trim(spectype) == 's-organic') then
190 30720 : call rad_cnst_get_info_by_bin_spec(0, m, l, spec_name=spec_name)
191 30720 : ns = ns + 1
192 30720 : dqdtsoa_idx(m,ns) = pbuf_get_index('DQDT_'//trim(spec_name))
193 : end if
194 : end do
195 : end do
196 :
197 62976 : do m = 1, nbins
198 62976 : do_soag_any(m) = cnsoa(m)>0
199 : end do
200 :
201 : !---------define history fields for new cond/evap diagnostics----------------------------------------
202 :
203 1536 : fieldname=trim('qcon_gaex')
204 1536 : long_name = trim('3D fields for SOA condensation')
205 1536 : unit = 'kg/kg/s'
206 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
207 1536 : if ( history_aerosol ) then
208 0 : call add_default( fieldname, 1, ' ' )
209 : endif
210 :
211 3072 : do j = 1, nsoa_vbs
212 1536 : write (outsoa, "(I2.2)") j
213 1536 : fieldname = 'qcon_gaex'//outsoa
214 1536 : long_name = '3D fields for SOA condensation for VBS bin'//outsoa
215 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
216 1536 : if ( history_aerosol ) then
217 0 : call add_default( fieldname, 1, ' ' )
218 : endif
219 1536 : fieldname = 'qevap_gaex'//outsoa
220 1536 : long_name = '3D fields for SOA evaporation for VBS bin'//outsoa
221 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
222 3072 : if ( history_aerosol ) then
223 0 : call add_default( fieldname, 1, ' ' )
224 : endif
225 : end do
226 :
227 1536 : fieldname=trim('qevap_gaex')
228 1536 : long_name = trim('3D fields for SOA evaporation')
229 3072 : call addfld(fieldname, (/'lev'/), 'A', unit, long_name )
230 1536 : if ( history_aerosol ) then
231 0 : call add_default( fieldname, 1, ' ' )
232 : endif
233 :
234 : !------------------------------------------------------------------------------
235 :
236 : ! define history fields for basic gas-aer exchange
237 62976 : do m = 1, nbins
238 276480 : do l = 1, nspec(m) ! do through nspec
239 215040 : ii = bin_idx(m,l)
240 276480 : if (l <= nspec(m) ) then ! species
241 215040 : call rad_cnst_get_info_by_bin_spec(0, m, l, spec_name=fldname(ii) )
242 : ! only write out SOA exchange here
243 215040 : call rad_cnst_get_bin_props_by_idx(0, m, l,spectype=spectype)
244 215040 : if (trim(spectype) == 's-organic') then
245 30720 : fieldname= trim(fldname(ii)) // '_sfgaex1'
246 30720 : long_name = trim(fldname(ii)) // ' gas-aerosol-exchange primary column tendency'
247 30720 : unit = 'kg/m2/s'
248 30720 : call addfld( fieldname, horiz_only, 'A', unit, long_name )
249 30720 : if ( history_aerosol ) then
250 0 : call add_default( fieldname, 1, ' ' )
251 : endif
252 : end if
253 : end if
254 : end do
255 :
256 61440 : write(fieldname,'("WETRAD_bin",I2.2)') m
257 61440 : write(long_name,'("bin ",I2.2," wet radius in carma_aero_gasaerexch")') m
258 :
259 122880 : call addfld(fieldname, (/'lev'/), 'A', 'cm', long_name )
260 61440 : if ( history_aerosol ) then
261 0 : call add_default( fieldname, 1, ' ' )
262 : endif
263 :
264 61440 : write(fieldname,'("UPTKRATE_bin",I2.2)') m
265 61440 : write(long_name,'("bin ",I2.2," up take rate in carma_aero_gasaerexch")') m
266 :
267 122880 : call addfld(fieldname, (/'lev'/), 'A', 'sec-1', long_name )
268 61440 : if ( history_aerosol ) then
269 0 : call add_default( fieldname, 1, ' ' )
270 : endif
271 :
272 61440 : write(fieldname,'("NUMDENS_bin",I2.2)') m
273 61440 : write(long_name,'("bin ",I2.2," number density carma_aero_gasaerexch")') m
274 :
275 122880 : call addfld(fieldname, (/'lev'/), 'A', 'm-3', long_name )
276 62976 : if ( history_aerosol ) then
277 0 : call add_default( fieldname, 1, ' ' )
278 : endif
279 :
280 : end do
281 :
282 1536 : fieldname=trim('UPTKRATE')
283 1536 : long_name = trim('total uptake rate in carma_aero_gasaerexch')
284 3072 : call addfld(fieldname, (/'lev'/), 'A', 'sec-1', long_name )
285 1536 : if ( history_aerosol ) then
286 0 : call add_default( fieldname, 1, ' ' )
287 : endif
288 :
289 :
290 3072 : end subroutine carma_aero_gasaerexch_init
291 :
292 :
293 : !----------------------------------------------------------------------
294 :
295 : !----------------------------------------------------------------------
296 : !----------------------------------------------------------------------
297 : !BOP
298 : ! !ROUTINE: carma_aero_gasaerexch_sub --- ...
299 : !
300 : ! !INTERFACE:
301 72960 : subroutine carma_aero_gasaerexch_sub( state, &
302 : pbuf, lchnk, ncol, nstep, &
303 72960 : loffset, deltat, mbar, &
304 : t, pmid, pdel, &
305 : qh2o, troplev, &
306 72960 : q, raervmr, &
307 72960 : wetr_n )
308 :
309 : ! !USES:
310 1536 : use cam_history, only: outfld, fieldname_len
311 : use physconst, only: gravit, mwdry
312 : use cam_abortutils, only: endrun
313 : use time_manager, only: is_first_step
314 : use carma_aerosol_state_mod, only: carma_aerosol_state
315 : use physics_types, only: physics_state
316 : use physconst, only: mwdry, rair
317 :
318 : ! !PARAMETERS:
319 : type(physics_state), target, intent(in) :: state ! Physics state variables
320 : type(physics_buffer_desc), pointer :: pbuf(:)
321 :
322 : integer, intent(in) :: lchnk ! chunk identifier
323 : integer, intent(in) :: ncol ! number of atmospheric column
324 : integer, intent(in) :: nstep ! model time-step number
325 : integer, intent(in) :: loffset ! offset applied to modal aero "ptrs"
326 : integer, intent(in) :: troplev(pcols) ! tropopause vertical index
327 : real(r8), intent(in) :: deltat ! time step (s)
328 : real(r8), intent(in) :: mbar(ncol,pver) ! mean wet atmospheric mass ( amu )
329 :
330 : real(r8), intent(inout) :: q(ncol,pver,gas_pcnst) ! tracer mixing ratio (TMR) array
331 : ! *** MUST BE #/kmol-air for number
332 : ! *** MUST BE mol/mol-air for mass
333 : ! *** NOTE ncol dimension
334 : real(r8), intent(in) :: raervmr (ncol,pver,ncnst_tot) ! aerosol mixing rations (vmr)
335 : real(r8), intent(in) :: t(pcols,pver) ! temperature at model levels (K)
336 : real(r8), intent(in) :: pmid(pcols,pver) ! pressure at model levels (Pa)
337 : real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickness of levels (Pa)
338 : real(r8), intent(in) :: qh2o(pcols,pver) ! water vapor mixing ratio (kg/kg)
339 : real(r8), intent(in) :: wetr_n(pcols,pver,nbins) !wet geo. mean dia. (cm) of number distrib.
340 :
341 : ! !DESCRIPTION:
342 : ! this version does only do condensation for SOA for CARMA
343 : ! method_soa=0 is no uptake
344 : ! method_soa=1 is irreversible uptake done like h2so4 uptake
345 : ! method_soa=2 is reversible uptake using subr carma_aero_soaexch
346 : !
347 : ! !REVISION HISTORY:
348 : ! RCE 07.04.13: Adapted from MIRAGE2 code
349 : !
350 : !EOP
351 : !----------------------------------------------------------------------
352 : !BOC
353 :
354 : ! local variables
355 : integer, parameter :: ldiag1=-1, ldiag2=-1, ldiag3=-1, ldiag4=-1
356 : integer, parameter :: method_soa = 2
357 :
358 : real (r8), parameter :: mw_poa_host = 12.0_r8 ! molec wght of poa used in host code
359 : real (r8), parameter :: mw_soa_host = 250.0_r8 ! molec wght of soa used in host code
360 :
361 : integer :: i
362 : integer :: j, jsoa
363 : integer :: k
364 : integer :: l
365 : integer :: mm, m, n, nn, niter, niter_max
366 :
367 : character(len=fieldname_len+3) :: fieldname
368 : character(len=32) :: spectype
369 : character(len=2) :: outsoa
370 :
371 145920 : real (r8) :: avg_uprt_soa(nsoa_vbs)
372 : real (r8) :: deltatxx
373 145920 : real (r8) :: dqdt_soa_vbs(nbins,nsoa_vbs)
374 145920 : real (r8) :: dqdt_soa_all(pcols,pver,nbins,nsoa)
375 145920 : real (r8) :: dqdt_soag(nsoa_vbs)
376 145920 : real (r8) :: fgain_soa(nbins,nsoa_vbs)
377 : real (r8) :: pdel_fac
378 145920 : real (r8) :: num_bin(pcols,pver,nbins)
379 145920 : real (r8) :: soa_vbs(pcols,pver,nbins,nsoa_vbs)
380 145920 : real (r8) :: soa_c(pcols,pver,nbins,nsoa) ! SOA from CARMA
381 145920 : real (r8) :: poa_c(pcols,pver,nbins,npoa) ! POA from CARMA
382 145920 : real (r8) :: qold_poa(nbins,npoa) ! POA from CARMA old
383 145920 : real (r8) :: qold_soa(nbins,nsoa_vbs) ! SOA on VBS bins old
384 145920 : real (r8) :: qnew_soa_vbs(nbins,nsoa_vbs) ! SOA on VBS bins new
385 145920 : real (r8) :: qnew_soa(nbins) ! SOA new for combined VBS bin new for combined VBS binss
386 145920 : real (r8) :: qold_soag(nsoa_vbs)
387 145920 : real (r8) :: sum_dqdt_soa(nsoa_vbs) ! sum_dqdt_soa = soa tendency from soa gas uptake (mol/mol/s)
388 145920 : real (r8) :: sum_uprt_soa(nsoa_vbs) ! total soa uptake rate over all bin, for each soa vbs bin
389 145920 : real (r8) :: uptkrate(pcols,pver,nbins)
390 : real (r8) :: uptkrate_all(pcols,pver)
391 145920 : real (r8) :: uptkratebb(nbins)
392 145920 : real (r8) :: uptkrate_soa(nbins,nsoa_vbs)
393 : ! gas-to-aerosol mass transfer rates (1/s)
394 :
395 : integer, parameter :: nsrflx = 1 ! only one dimension of qsrflx, no renaming or changes in size for CARMA currently
396 145920 : real(r8) :: dqdt(ncol,pver,gas_pcnst) ! TMR "delta q" array - NOTE dims
397 145920 : real(r8) :: qsrflx(pcols,nbins,nsoa)
398 : ! process-specific column tracer tendencies
399 : ! (1=gas condensation)
400 145920 : real(r8) :: qcon_vbs(pcols,pver,nsoa_vbs)
401 145920 : real(r8) :: qevap_vbs(pcols,pver,nsoa_vbs)
402 : real(r8) :: qcon(pcols,pver)
403 : real(r8) :: qevap(pcols,pver)
404 : real(r8) :: total_soag
405 145920 : real(r8) :: soag(nsoa_vbs)
406 :
407 72960 : real(r8), pointer :: frac_vbs(:,:,:,:) ! fraction of vbs SOA bins to total SOA
408 72960 : real(r8), pointer :: dqdt_soa(:,:)
409 :
410 : real(r8) :: rhoair(pcols,pver)
411 72960 : real(r8), pointer :: nmr(:,:)
412 : type(carma_aerosol_state), pointer :: aero_state
413 :
414 : !----------------------------------------------------------------------
415 145920 : aero_state => carma_aerosol_state(state, pbuf)
416 :
417 : ! map CARMA soa to working soa(nbins,nsoa)
418 :
419 72960 : call pbuf_get_field(pbuf, fracvbs_idx, frac_vbs)
420 :
421 3475887360 : num_bin(:,:,:) = 0._r8
422 3475960320 : soa_c(:,:,:,:) = 0._r8
423 3475960320 : poa_c(:,:,:,:) = 0._r8
424 :
425 78723840 : rhoair(:ncol,:) = pmid(:ncol,:)/(rair*t(:ncol,:)) ! (kg-air/m3)
426 :
427 2991360 : do m = 1, nbins ! main loop over aerosol bins
428 2991360 : if (do_soag_any(m)) then ! only bins that contain soa
429 1459200 : n = 0
430 1459200 : nn = 0
431 10214400 : do l = 1, nspec(m)
432 8755200 : mm = bin_idx(m, l)
433 8755200 : call rad_cnst_get_bin_props_by_idx(0, m, l, spectype=spectype)
434 8755200 : if (trim(spectype) == 's-organic') then
435 1459200 : n = n + 1
436 1574476800 : soa_c(:ncol,:,m,n) = raervmr(:ncol,:,mm)
437 : end if
438 10214400 : if (trim(spectype) == 'p-organic') then
439 1459200 : nn = nn + 1
440 1574476800 : poa_c(:ncol,:,m,nn) = raervmr(:ncol,:,mm)
441 : end if
442 : end do
443 1459200 : if (npoa .gt. 1) then
444 0 : call endrun( 'carma_aero_gasaerexch_sub error: CARMA currently only supports 1 POA element' )
445 : end if
446 :
447 1459200 : if (nsoa_vbs.eq.nsoa) then
448 62981990400 : soa_vbs(:ncol,:,:,:) = soa_c(:ncol,:,:,:)
449 : else
450 0 : if (nsoa.eq.1) then
451 0 : if (is_first_step()) then
452 : !first time step initialization only
453 0 : do k=top_lev,pver
454 0 : do i=1,ncol
455 0 : total_soag = 0.0_r8
456 0 : do j = 1, nsoa_vbs
457 0 : soag(j) = q(i,k,l_soag(j))
458 0 : total_soag = total_soag + soag(j)
459 : end do
460 0 : if (total_soag .gt. 0.0_r8) then
461 0 : do j= 1, nsoa_vbs
462 0 : frac_vbs(i,k,m,j) = soag(j)/total_soag
463 : end do
464 : end if
465 : end do
466 : end do
467 : end if
468 : ! end first time step, after that use fraction from previous time step
469 0 : do k=top_lev,pver
470 0 : do i=1,ncol
471 0 : do j= 1, nsoa_vbs
472 0 : soa_vbs(i,k,m,j) = frac_vbs(i,k,m,j)*soa_c(i,k,m,nsoa)
473 : end do
474 : end do
475 : end do
476 : else
477 : ! error message this code only works if SOAG and SOA CARMA have the same number of species,
478 : ! or if SOA CARMA has only one species.
479 0 : call endrun( 'carma_aero_gasaerexch_sub error in number of SOA species' )
480 : end if
481 :
482 : end if
483 :
484 : ! get bin number densities for all aerosols
485 1459200 : call aero_state%get_ambient_num(m,nmr) ! #/kg
486 1574476800 : num_bin(:ncol,:,m) = nmr(:ncol,:)*rhoair(:ncol,:) ! #/m3
487 :
488 : end if
489 : end do
490 :
491 :
492 : ! SOA will be updated in CARMA
493 :
494 : ! zero out tendencies and other
495 6376704000 : dqdt(:,:,:) = 0.0_r8
496 49758720 : qsrflx(:,:,:) = 0.0_r8
497 :
498 : !-------Initialize evap/cond diagnostics (ncols x pver)-----------
499 86968320 : qcon_vbs(:,:,:) = 0.0_r8
500 86968320 : qevap_vbs(:,:,:) = 0.0_r8
501 72960 : qcon(:,:) = 0.0_r8
502 72960 : qevap(:,:) = 0.0_r8
503 : !---------------------------------------------------
504 : ! compute gas-to-aerosol mass transfer rates
505 : ! check if only number is needed for this calculatuion!
506 : call gas_aer_uptkrates( ncol, loffset, &
507 : num_bin, t, pmid, &
508 72960 : wetr_n, uptkrate )
509 :
510 2991360 : do m = 1, nbins
511 :
512 2918400 : write(fieldname,'("NUMDENS_bin",I2.2)') m
513 3148953600 : call outfld(fieldname, num_bin(:ncol,:,m), ncol, lchnk )
514 :
515 2918400 : write(fieldname,'("WETRAD_bin",I2.2)') m
516 3148953600 : call outfld(fieldname, wetr_n(:ncol,:,m), ncol, lchnk )
517 :
518 2918400 : write(fieldname,'("UPTKRATE_bin",I2.2)') m
519 3148953600 : call outfld(fieldname, uptkrate(:ncol,:,m), ncol, lchnk )
520 :
521 3149026560 : uptkrate_all(:ncol,:) = uptkrate_all(:ncol,:) + uptkrate(:ncol,:,m)
522 : end do
523 :
524 72960 : fieldname = trim('UPTKRATE')
525 78723840 : call outfld(fieldname, uptkrate_all(:ncol,:), ncol, lchnk )
526 :
527 : ! use this for tendency calcs to avoid generating very small negative values
528 72960 : deltatxx = deltat * (1.0_r8 + 1.0e-15_r8)
529 :
530 3475960320 : dqdt_soa_all(:,:,:,:) = 0.0_r8
531 5180160 : do k=top_lev,pver
532 78723840 : do i=1,ncol
533 147087360 : sum_uprt_soa(:) = 0.0_r8
534 3088834560 : uptkrate_soa(:,:) = 0.0_r8
535 3015290880 : do n = 1, nbins
536 3015290880 : if (do_soag_any(n)) then ! only bins that contain soa
537 1470873600 : uptkratebb(n) = uptkrate(i,k,n)
538 1470873600 : if (npoa .gt. 0) then
539 2941747200 : do j = 1, npoa
540 2941747200 : qold_poa(n,j) = poa_c(i,k,n,j)
541 : end do
542 : else
543 0 : qold_poa(n,j) = 0.0_r8
544 : end if
545 2941747200 : do jsoa = 1, nsoa_vbs
546 : ! 0.81 factor is for gas diffusivity (soa/h2so4)
547 : ! (differences in fuch-sutugin and accom coef ignored)
548 1470873600 : fgain_soa(n,jsoa) = uptkratebb(n)*0.81_r8
549 1470873600 : uptkrate_soa(n,jsoa) = fgain_soa(n,jsoa)
550 1470873600 : sum_uprt_soa(jsoa) = sum_uprt_soa(jsoa) + fgain_soa(n,jsoa)
551 2941747200 : qold_soa(n,jsoa) = soa_vbs(i,k,n,jsoa)
552 : end do
553 : else
554 2941747200 : qold_poa(n,:) = 0.0_r8
555 2941747200 : qold_soa(n,:) = 0.0_r8
556 2941747200 : fgain_soa(n,:) = 0.0_r8
557 : end if
558 : end do ! n
559 :
560 147087360 : do jsoa = 1, nsoa_vbs
561 147087360 : if (sum_uprt_soa(jsoa) > 0.0_r8) then
562 3015290880 : do n = 1, nbins
563 3015290880 : if (do_soag_any(n)) then ! only bins that contain soa
564 1470873600 : fgain_soa(n,jsoa) = fgain_soa(n,jsoa) / sum_uprt_soa(jsoa)
565 : end if
566 : end do
567 : end if
568 : end do
569 :
570 : ! uptake amount (fraction of gas uptaken) over deltat
571 147087360 : do jsoa = 1, nsoa_vbs
572 147087360 : avg_uprt_soa(jsoa) = (1.0_r8 - exp(-deltatxx*sum_uprt_soa(jsoa)))/deltatxx
573 : end do
574 :
575 : ! sum_dqdt_soa = soa_a tendency from soa gas uptake (mol/mol/s)
576 :
577 147087360 : do jsoa = 1, nsoa_vbs
578 147087360 : sum_dqdt_soa(jsoa) = q(i,k,l_soag(jsoa)) * avg_uprt_soa(jsoa)
579 : end do
580 :
581 : if (method_soa > 1) then
582 : ! compute TMR tendencies for soag and soa interstial aerosol
583 : ! using soa parameterization
584 73543680 : niter_max = 1000
585 3088834560 : dqdt_soa_vbs(:,:) = 0.0_r8
586 147087360 : dqdt_soag(:) = 0.0_r8
587 147087360 : do jsoa = 1, nsoa_vbs
588 147087360 : qold_soag(jsoa) = q(i,k,l_soag(jsoa))
589 : end do
590 :
591 147087360 : call carma_aero_soaexch( deltat, t(i,k), pmid(i,k), &
592 : niter, niter_max, nbins, nsoa_vbs, npoa, &
593 : mw_poa_host, mw_soa_host, &
594 : qold_soag, qold_soa, qold_poa, uptkrate_soa, &
595 220631040 : dqdt_soag, dqdt_soa_vbs )
596 :
597 147087360 : sum_dqdt_soa(:) = -dqdt_soag(:)
598 :
599 : else if ( method_soa .eq. 1) then
600 : ! compute TMR tendencies for soa interstial aerosol
601 : ! due to simple gas uptake
602 :
603 : do n = 1, nbins
604 : if (do_soag_any(n) ) then
605 : do jsoa = 1, nsoa_vbs
606 : dqdt_soa_vbs(n,jsoa) = fgain_soa(n,jsoa)*sum_dqdt_soa(jsoa)
607 : end do
608 : end if
609 : end do
610 :
611 : end if
612 :
613 : ! update soa to calcuate fractions (state variables and pbuf is not updated for SOA, will be done in CARMA)
614 73543680 : pdel_fac = pdel(i,k)/gravit
615 3015290880 : qnew_soa(:) =0.0_r8
616 3088834560 : qnew_soa_vbs(:,:) =0.0_r8
617 :
618 3015290880 : do n = 1, nbins
619 3015290880 : if ( do_soag_any(n) ) then
620 1470873600 : if (nsoa.eq.nsoa_vbs) then
621 2941747200 : do jsoa = 1, nsoa_vbs
622 1470873600 : qsrflx(i,n,jsoa) = qsrflx(i,n,jsoa) + dqdt_soa_vbs(n,jsoa)*pdel_fac
623 2941747200 : dqdt_soa_all(i,k,n,jsoa) = dqdt_soa_vbs(n,jsoa) ! sum up for different volatility bins
624 : end do
625 0 : else if (nsoa.eq.1) then
626 0 : do jsoa = 1, nsoa_vbs
627 : ! sum up for different volatility bins
628 0 : dqdt_soa_all(i,k,n,nsoa) = dqdt_soa_all(i,k,n,nsoa) + dqdt_soa_vbs(n,jsoa)
629 : end do
630 0 : do jsoa = 1, nsoa_vbs
631 0 : qsrflx(i,n,nsoa) = qsrflx(i,n,nsoa) + dqdt_soa_vbs(n,jsoa)*pdel_fac
632 0 : qnew_soa_vbs(n,jsoa) = qold_soa(n,jsoa) + dqdt_soa_vbs(n,jsoa)*deltat
633 0 : qnew_soa(n) = qnew_soa(n) + qnew_soa_vbs(n,jsoa) ! derive new fraction of SOA bin contributions
634 : end do
635 0 : do jsoa = 1, nsoa_vbs
636 0 : if (qnew_soa(n) .gt. 0.0_r8) then
637 0 : frac_vbs(i,k,n,jsoa) = qnew_soa_vbs(n,jsoa) / qnew_soa(n)
638 : end if
639 : end do
640 : else
641 0 : call endrun( 'carma_aero_gasaerexch_sub error' )
642 : end if
643 :
644 : !------- Add code for condensation/evaporation diagnostics sum of all bin---
645 2941747200 : do jsoa = 1, nsoa_vbs
646 2941747200 : if (dqdt_soa_vbs(n,jsoa).ge.0.0_r8) then
647 1300721000 : qcon_vbs(i,k,jsoa)=qcon_vbs(i,k,jsoa) + dqdt_soa_vbs(n,jsoa)*(mw_soa/mwdry)
648 1300721000 : qcon(i,k)=qcon(i,k)+dqdt_soa_vbs(n,jsoa)*(mw_soa/mwdry)
649 170152600 : else if (dqdt_soa_vbs(n,jsoa).lt.0.0_r8) then
650 170152600 : qevap_vbs(i,k,jsoa)=qevap_vbs(i,k,jsoa) + dqdt_soa_vbs(n,jsoa)*(mw_soa/mwdry)
651 170152600 : qevap(i,k)=qevap(i,k)+dqdt_soa_vbs(n,jsoa)*(mw_soa/mwdry)
652 : endif
653 : end do
654 : !---------------------------------------------------------------------------------------------------------------------
655 : end if
656 : end do ! n
657 :
658 : ! compute TMR tendencies for SAOG gas
659 : ! due to simple gas uptake
660 152194560 : do jsoa = 1, nsoa
661 147087360 : dqdt(i,k,l_soag(jsoa)) = -sum_dqdt_soa(jsoa)
662 : end do
663 :
664 : end do ! "i = 1, ncol"
665 : end do ! "k = top_lev, pver"
666 :
667 : ! This applies dqdt tendencies for SOAG only , soa is done in CARMA
668 : ! apply the dqdt to update q
669 : !
670 145920 : do jsoa = 1, nsoa_vbs
671 5253120 : do k = top_lev, pver
672 78723840 : do i = 1, ncol
673 78650880 : q(i,k,l_soag(jsoa)) = max (q(i,k,l_soag(jsoa)) + dqdt(i,k,l_soag(jsoa))*deltat, 1.0e-40_r8)
674 : end do
675 : end do
676 : end do
677 :
678 :
679 : !-----Outfld for condensation/evaporation------------------------------
680 72960 : call outfld(trim('qcon_gaex'), qcon(:,:), pcols, lchnk )
681 72960 : call outfld(trim('qevap_gaex'), qevap(:,:), pcols, lchnk )
682 145920 : do jsoa = 1, nsoa_vbs
683 72960 : write (outsoa, "(I2.2)") jsoa
684 72960 : call outfld(trim('qcon_gaex')//outsoa, qcon_vbs(:,:,jsoa), pcols, lchnk )
685 145920 : call outfld(trim('qevap_gaex')//outsoa, qevap_vbs(:,:,jsoa), pcols, lchnk )
686 : end do
687 : !-----------------------------------------------------------------------
688 : ! do history file of column-tendency fields over SOA fields (as defined in CARMA) and set pointer
689 2991360 : do m = 1, nbins
690 2991360 : if (do_soag_any(m)) then
691 1459200 : j = 0
692 10214400 : do l = 1, nspec(m)
693 8755200 : mm = bin_idx(m,l)
694 8755200 : call rad_cnst_get_bin_props_by_idx(0, m, l,spectype=spectype)
695 10214400 : if (trim(spectype) == 's-organic') then
696 1459200 : j = j + 1
697 1459200 : fieldname= trim(fldname(mm)) // '_sfgaex1'
698 22471680 : do i = 1, ncol
699 22471680 : qsrflx(i,m,j) = qsrflx(i,m,j)*(mw_soa/mwdry)
700 : end do
701 1459200 : call outfld( fieldname, qsrflx(:,m,j), pcols, lchnk )
702 :
703 : !set pointer field
704 1459200 : call pbuf_get_field(pbuf, dqdtsoa_idx(m,j), dqdt_soa )
705 :
706 1574476800 : dqdt_soa(:ncol,:) = dqdt_soa_all(:ncol,:,m,j) *(mw_soa/mbar(:ncol,:))
707 : end if
708 : end do ! l = ...
709 : end if
710 : end do ! m = ...
711 :
712 145920 : end subroutine carma_aero_gasaerexch_sub
713 :
714 :
715 : !----------------------------------------------------------------------
716 : !----------------------------------------------------------------------
717 72960 : subroutine gas_aer_uptkrates( ncol, loffset, &
718 72960 : num_bin, t, pmid, &
719 72960 : wetr, uptkrate )
720 :
721 : !
722 : ! /
723 : ! computes uptkrate = | dx dN/dx gas_conden_rate(Dp(x))
724 : ! /
725 : ! using Gauss-Hermite quadrature of order nghq=2
726 : !
727 : ! Dp = particle diameter (cm)
728 : ! x = ln(Dp)
729 : ! dN/dx = log-normal particle number density distribution
730 : ! gas_conden_rate(Dp) = 2 * pi * gasdiffus * Dp * F(Kn,ac)
731 : ! F(Kn,ac) = Fuchs-Sutugin correction factor
732 : ! Kn = Knudsen number
733 : ! ac = accomodation coefficient
734 : !
735 :
736 : integer, intent(in) :: ncol ! number of atmospheric column
737 : integer, intent(in) :: loffset
738 : real(r8), intent(in) :: t(pcols,pver) ! Temperature in Kelvin
739 : real(r8), intent(in) :: pmid(pcols,pver) ! Air pressure in Pa
740 : real(r8), intent(in) :: wetr(pcols,pver,nbins)
741 : real(r8), intent(in) :: num_bin(pcols,pver,nbins)
742 :
743 : real(r8), intent(out) :: uptkrate(pcols,pver,nbins)
744 : ! gas-to-aerosol mass transfer rates (1/s)
745 :
746 :
747 : ! local
748 : integer, parameter :: nghq = 2
749 : integer :: i, k, n
750 :
751 : ! Can use sqrt here once Lahey is gone.
752 : real(r8), parameter :: tworootpi = 3.5449077_r8
753 : real(r8), parameter :: root2 = 1.4142135_r8
754 : real(r8), parameter :: beta = 2.0_r8
755 :
756 : real(r8) :: const
757 : real(r8) :: dp
758 : real(r8) :: gasdiffus, gasspeed
759 : real(r8) :: freepathx2, fuchs_sutugin
760 : real(r8) :: knudsen
761 :
762 : ! initialize to zero
763 3475887360 : uptkrate(:,:,:) = 0.0_r8
764 :
765 : ! outermost loop over all bins
766 2991360 : do n = 1, nbins
767 :
768 : ! loops k and i
769 207279360 : do k=top_lev,pver
770 3148953600 : do i=1,ncol
771 3146035200 : if (wetr(i,k,n) .gt. 0.0_r8) then
772 :
773 : ! gasdiffus = h2so4 gas diffusivity from mosaic code (m^2/s)
774 : ! (pmid must be Pa)
775 2861412147 : gasdiffus = 0.557e-4_r8 * (t(i,k)**1.75_r8) / pmid(i,k)
776 : ! gasspeed = h2so4 gas mean molecular speed from mosaic code (m/s)
777 2861412147 : gasspeed = 1.470e1_r8 * sqrt(t(i,k))
778 : ! freepathx2 = 2 * (h2so4 mean free path) (m)
779 2861412147 : freepathx2 = 6.0_r8*gasdiffus/gasspeed
780 2861412147 : dp = wetr(i,k,n) * 1.e-2_r8 ! meters
781 2861412147 : const = tworootpi * num_bin(i,k,n) * 2.0_r8 * dp
782 : ! gas_conden_rate(Dp) = const * gasdiffus * F(Kn,ac)
783 : ! knudsen number
784 2861412147 : knudsen = freepathx2/dp
785 : fuchs_sutugin = (0.4875_r8*(1._r8 + knudsen)) / &
786 2861412147 : (knudsen*(1.184_r8 + knudsen) + 0.4875_r8)
787 2861412147 : uptkrate(i,k,n) = const * gasdiffus * fuchs_sutugin
788 :
789 : else
790 80335053 : uptkrate(i,k,n) = 0.0_r8
791 : end if
792 :
793 : end do ! "do i = 1, ncol"
794 : end do ! "do k = 1, pver"
795 :
796 : end do ! "do n = 1, nbins"
797 :
798 72960 : end subroutine gas_aer_uptkrates
799 :
800 : !----------------------------------------------------------------------
801 73543680 : subroutine carma_aero_soaexch( dtfull, temp, pres, &
802 : niter, niter_max, nbins, ntot_soaspec, ntot_poaspec, &
803 : mw_poa_host, mw_soa_host, &
804 73543680 : g_soa_in, a_soa_in, a_poa_in, xferrate_in, &
805 73543680 : g_soa_tend, a_soa_tend )
806 :
807 : !-----------------------------------------------------------------------
808 : !
809 : ! Purpose:
810 : !
811 : ! calculates condensation/evaporation of "soa gas"
812 : ! to/from multiple aerosol modes in 1 grid cell
813 : !
814 : ! key assumptions
815 : ! (1) ambient equilibrium vapor pressure of soa gas
816 : ! is given by p0_soa_298 and delh_vap_soa
817 : ! (2) equilibrium vapor pressure of soa gas at aerosol
818 : ! particle surface is given by raoults law in the form
819 : ! g_star = g0_soa*[a_soa/(a_soa + a_opoa)]
820 : ! (3) (oxidized poa)/(total poa) is equal to frac_opoa (constant)
821 : !
822 : !
823 : ! Author: R. Easter and R. Zaveri
824 : ! Additions to run with multiple BC, SOA and POM's: Shrivastava et al., 2015
825 : !-----------------------------------------------------------------------
826 :
827 : use mo_constants, only: rgas ! Gas constant (J/K/mol)
828 :
829 : real(r8), intent(in) :: dtfull ! full integration time step (s)
830 : real(r8), intent(in) :: temp ! air temperature (K)
831 : real(r8), intent(in) :: pres ! air pressure (Pa)
832 : integer, intent(out) :: niter ! number of iterations performed
833 : integer, intent(in) :: niter_max ! max allowed number of iterations
834 : integer, intent(in) :: nbins ! number of bins
835 : integer, intent(in) :: ntot_poaspec ! number of poa species
836 : integer, intent(in) :: ntot_soaspec ! number of soa species
837 : real(r8), intent(in) :: mw_poa_host ! molec wght of poa used in host code
838 : real(r8), intent(in) :: mw_soa_host ! molec wght of soa used in host code
839 : real(r8), intent(in) :: g_soa_in(ntot_soaspec) ! initial soa gas mixrat (mol/mol at host mw)
840 : real(r8), intent(in) :: a_soa_in(nbins,ntot_soaspec) ! initial soa aerosol mixrat (mol/mol at host mw)
841 : real(r8), intent(in) :: a_poa_in(nbins,ntot_poaspec) ! initial poa aerosol mixrat (mol/mol at host mw)
842 : real(r8), intent(in) :: xferrate_in(nbins,ntot_soaspec) ! gas-aerosol mass transfer rate (1/s)
843 : real(r8), intent(out) :: g_soa_tend(ntot_soaspec) ! soa gas mixrat tendency (mol/mol/s at host mw)
844 : real(r8), intent(out) :: a_soa_tend(nbins,ntot_soaspec) ! soa aerosol mixrat tendency (mol/mol/s at host mw)
845 :
846 : integer :: ll
847 : integer :: m
848 :
849 147087360 : logical :: skip_soamode(nbins) ! true if this bin does not have soa
850 :
851 : real(r8), parameter :: a_min1 = 1.0e-40_r8
852 : real(r8), parameter :: g_min1 = 1.0e-40_r8
853 : real(r8), parameter :: alpha = 0.05_r8 ! parameter used in calc of time step
854 : real(r8), parameter :: dtsub_fixed = -1.0_r8 ! fixed sub-step for time integration (s)
855 :
856 147087360 : real(r8) :: a_ooa_sum_tmp(nbins) ! total ooa (=soa+opoa) in a bin
857 147087360 : real(r8) :: a_opoa(nbins) ! oxidized-poa aerosol mixrat (mol/mol at actual mw)
858 147087360 : real(r8) :: a_soa(nbins,ntot_soaspec) ! soa aerosol mixrat (mol/mol at actual mw)
859 147087360 : real(r8) :: a_soa_tmp(nbins,ntot_soaspec) ! temporary soa aerosol mixrat (mol/mol)
860 147087360 : real(r8) :: beta(nbins,ntot_soaspec) ! dtcur*xferrate
861 147087360 : real(r8) :: delh_vap_soa(ntot_soaspec) ! delh_vap_soa = heat of vaporization for gas soa (J/mol)
862 147087360 : real(r8) :: del_g_soa_tmp(ntot_soaspec)
863 : real(r8) :: dtcur ! current time step (s)
864 : real(r8) :: dtmax ! = (dtfull-tcur)
865 147087360 : real(r8) :: g0_soa(ntot_soaspec) ! ambient soa gas equilib mixrat (mol/mol at actual mw)
866 147087360 : real(r8) :: g_soa(ntot_soaspec) ! soa gas mixrat (mol/mol at actual mw)
867 147087360 : real(r8) :: g_star(nbins,ntot_soaspec) ! soa gas mixrat that is in equilib
868 : ! with each aerosol mode (mol/mol)
869 : real(r8) :: mw_poa ! actual molec wght of poa
870 : real(r8) :: mw_soa ! actual molec wght of soa
871 147087360 : real(r8) :: opoa_frac(ntot_poaspec) ! fraction of poa that is opoa
872 147087360 : real(r8) :: phi(nbins,ntot_soaspec) ! "relative driving force"
873 147087360 : real(r8) :: p0_soa(ntot_soaspec) ! soa gas equilib vapor presssure (atm)
874 147087360 : real(r8) :: p0_soa_298(ntot_soaspec) ! p0_soa_298 = soa gas equilib vapor presssure (atm) at 298 k
875 147087360 : real(r8) :: sat(nbins,ntot_soaspec) ! sat(m,ll) = g0_soa(ll)/a_ooa_sum_tmp(m) = g_star(m,ll)/a_soa(m,ll)
876 : ! used by the numerical integration scheme -- it is not a saturation rato!
877 : real(r8) :: tcur ! current integration time (from 0 s)
878 : real(r8) :: tmpa, tmpb, tmpf
879 147087360 : real(r8) :: tot_soa(ntot_soaspec) ! g_soa + sum( a_soa(:) )
880 73543680 : real(r8) :: xferrate(nbins,ntot_soaspec) ! gas-aerosol mass transfer rate (1/s)
881 :
882 : ! Changed by Manish Shrivastava
883 147087360 : opoa_frac(:) = 0.0_r8 !POA does not form solution with SOA for all runs; set opoa_frac=0.0_r8 by Manish Shrivastava
884 73543680 : mw_poa = 250.0_r8
885 73543680 : mw_soa = 250.0_r8
886 :
887 : ! New SOA properties added by Manish Shrivastava on 09/27/2012
888 73543680 : if (ntot_soaspec ==1) then
889 147087360 : p0_soa_298(:) = 1.0e-12_r8
890 147087360 : delh_vap_soa(:) = 156.0e3_r8
891 147087360 : opoa_frac(:) = 0.0_r8
892 0 : elseif (ntot_soaspec ==2) then
893 : ! same for anthropogenic and biomass burning species
894 0 : p0_soa_298 (1) = 1.0e-10_r8
895 0 : p0_soa_298 (2) = 1.0e-10_r8
896 0 : delh_vap_soa(:) = 156.0e3_r8
897 0 : elseif(ntot_soaspec ==5) then
898 : ! 5 volatility bins for each of the a combined SOA classes ( including biomass burning, fossil fuel, biogenic)
899 0 : p0_soa_298 (1) = 9.7831E-13_r8 !soaff0 C*=0.01ug/m3
900 0 : p0_soa_298 (2) = 9.7831E-12_r8 !soaff1 C*=0.10ug/m3
901 0 : p0_soa_298 (3) = 9.7831E-11_r8 !soaff2 C*=1.0ug/m3
902 0 : p0_soa_298 (4) = 9.7831E-10_r8 !soaff3 C*=10.0ug/m3
903 0 : p0_soa_298 (5) = 9.7831E-9_r8 !soaff4 C*=100.0ug/m3
904 :
905 0 : delh_vap_soa(1) = 153.0e3_r8
906 0 : delh_vap_soa(2) = 142.0e3_r8
907 0 : delh_vap_soa(3) = 131.0e3_r8
908 0 : delh_vap_soa(4) = 120.0e3_r8
909 0 : delh_vap_soa(5) = 109.0e3_r8
910 0 : elseif(ntot_soaspec ==15) then
911 : !
912 : ! 5 volatility bins for each of the 3 SOA classes ( biomass burning, fossil fuel, biogenic)
913 : ! SOA species 1-5 are for anthropogenic while 6-10 are for biomass burning SOA
914 : ! SOA species 11-15 are for biogenic SOA, based on Cappa et al., Reference needs to be updated
915 : ! For MW=250.0
916 0 : p0_soa_298 (1) = 9.7831E-13_r8 !soaff0 C*=0.01ug/m3
917 0 : p0_soa_298 (2) = 9.7831E-12_r8 !soaff1 C*=0.10ug/m3
918 0 : p0_soa_298 (3) = 9.7831E-11_r8 !soaff2 C*=1.0ug/m3
919 0 : p0_soa_298 (4) = 9.7831E-10_r8 !soaff3 C*=10.0ug/m3
920 0 : p0_soa_298 (5) = 9.7831E-9_r8 !soaff4 C*=100.0ug/m3
921 0 : p0_soa_298 (6) = 9.7831E-13_r8 !soabb0 C*=0.01ug/m3
922 0 : p0_soa_298 (7) = 9.7831E-12_r8 !soabb1 C*=0.10ug/m3
923 0 : p0_soa_298 (8) = 9.7831E-11_r8 !soabb2 C*=1.0ug/m3
924 0 : p0_soa_298 (9) = 9.7831E-10_r8 !soabb3 C*=10.0ug/m3
925 0 : p0_soa_298 (10) = 9.7831E-9_r8 !soabb4 C*=100.0ug/m3
926 0 : p0_soa_298 (11) = 9.7831E-13_r8 !soabg0 C*=0.01ug/m3
927 0 : p0_soa_298 (12) = 9.7831E-12_r8 !soabg1 C*=0.1ug/m3
928 0 : p0_soa_298 (13) = 9.7831E-11_r8 !soabg2 C*=1.0ug/m3
929 0 : p0_soa_298 (14) = 9.7831E-10_r8 !soabg3 C*=10.0ug/m3
930 0 : p0_soa_298 (15) = 9.7831E-9_r8 !soabg4 C*=100.0ug/m3
931 :
932 : !
933 : ! have to be adjusted to 15 species, following the numbers by Epstein et al., 2012
934 : !
935 0 : delh_vap_soa(1) = 153.0e3_r8
936 0 : delh_vap_soa(2) = 142.0e3_r8
937 0 : delh_vap_soa(3) = 131.0e3_r8
938 0 : delh_vap_soa(4) = 120.0e3_r8
939 0 : delh_vap_soa(5) = 109.0e3_r8
940 0 : delh_vap_soa(6) = 153.0e3_r8
941 0 : delh_vap_soa(7) = 142.0e3_r8
942 0 : delh_vap_soa(8) = 131.0e3_r8
943 0 : delh_vap_soa(9) = 120.0e3_r8
944 0 : delh_vap_soa(10) = 109.0e3_r8
945 0 : delh_vap_soa(11) = 153.0e3_r8
946 0 : delh_vap_soa(12) = 142.0e3_r8
947 0 : delh_vap_soa(13) = 131.0e3_r8
948 0 : delh_vap_soa(14) = 120.0e3_r8
949 0 : delh_vap_soa(15) = 109.0e3_r8
950 : endif
951 :
952 : !BSINGH - Initialized g_soa_tend and a_soa_tend to circumvent the undefined behavior (04/16/12)
953 147087360 : g_soa_tend(:) = 0.0_r8
954 3088834560 : a_soa_tend(:,:) = 0.0_r8
955 3088834560 : xferrate(:,:) = 0.0_r8
956 :
957 : ! determine which modes have non-zero transfer rates
958 : ! and are involved in the soa gas-aerosol transfer
959 : ! for diameter = 1 nm and number = 1 #/cm3, xferrate ~= 1e-9 s-1
960 3015290880 : do m = 1, nbins
961 3015290880 : if (do_soag_any(m)) then
962 1470873600 : skip_soamode(m) = .false.
963 2941747200 : do ll = 1, ntot_soaspec
964 2941747200 : xferrate(m,ll) = xferrate_in(m,ll)
965 : end do
966 : else
967 1470873600 : skip_soamode(m) = .true.
968 : end if
969 : end do
970 :
971 : ! convert incoming mixing ratios from mol/mol at the "host-code" molec. weight (12.0 in cam5)
972 : ! to mol/mol at the "actual" molec. weight (currently assumed to be 250.0)
973 : ! also
974 : ! force things to be non-negative
975 : ! calc tot_soa(ll)
976 : ! calc a_opoa (always slightly >0)
977 147087360 : do ll = 1, ntot_soaspec
978 73543680 : tmpf = mw_soa_host/mw_soa
979 73543680 : g_soa(ll) = max( g_soa_in(ll), 0.0_r8 ) * tmpf
980 73543680 : tot_soa(ll) = g_soa(ll)
981 3088834560 : do m = 1, nbins
982 2941747200 : if ( skip_soamode(m) ) cycle
983 1470873600 : a_soa(m,ll) = max( a_soa_in(m,ll), 0.0_r8 ) * tmpf
984 3015290880 : tot_soa(ll) = tot_soa(ll) + a_soa(m,ll)
985 : end do
986 : end do
987 :
988 :
989 3015290880 : tmpf = mw_poa_host/mw_poa
990 3015290880 : do m = 1, nbins
991 2941747200 : if ( skip_soamode(m) ) cycle
992 1470873600 : a_opoa(m) = 0.0_r8
993 : !check since it seems like in the modal approach there is a bug, not summing up the values for each specie
994 3015290880 : do ll = 1, ntot_poaspec
995 1470873600 : tmpf = mw_poa_host/mw_poa
996 1470873600 : a_opoa(m) = a_opoa(m) + opoa_frac(ll)*a_poa_in(m,ll)
997 4412620800 : a_opoa(m) = max( a_opoa(m), 1.0e-40_r8 ) ! force to small non-zero value
998 : end do
999 : end do
1000 :
1001 : ! calc ambient equilibrium soa gas
1002 147087360 : do ll = 1, ntot_soaspec
1003 73543680 : p0_soa(ll) = p0_soa_298(ll) * &
1004 73543680 : exp( -(delh_vap_soa(ll)/rgas)*((1.0_r8/temp)-(1.0_r8/298.0_r8)) )
1005 147087360 : g0_soa(ll) = 1.01325e5_r8*p0_soa(ll)/pres
1006 : end do
1007 :
1008 : ! IF mw of soa EQ 12 (as in the MAM3 default case), this has to be in
1009 : ! should actully talk the mw from the chemistry mechanism and substitute with 12.0
1010 :
1011 73543680 : niter = 0
1012 73543680 : tcur = 0.0_r8
1013 73543680 : dtcur = 0.0_r8
1014 3088834560 : phi(:,:) = 0.0_r8
1015 3088834560 : g_star(:,:) = 0.0_r8
1016 :
1017 : ! integration loop -- does multiple substeps to reach dtfull
1018 746672353 : time_loop: do while (tcur < dtfull-1.0e-3_r8 )
1019 :
1020 673129387 : niter = niter + 1
1021 673129387 : if (niter > niter_max) exit
1022 :
1023 27598275593 : tmpa = 0.0_r8 ! time integration parameter for all soa species
1024 27598275593 : do m = 1, nbins
1025 26925146920 : if ( skip_soamode(m) ) cycle
1026 41060849053 : a_ooa_sum_tmp(m) = sum( a_soa(m,1:ntot_soaspec) )
1027 : end do
1028 1346257346 : do ll = 1, ntot_soaspec
1029 : tmpb = 0.0_r8 ! time integration parameter for a single soa species
1030 27598275593 : do m = 1, nbins
1031 26925146920 : if ( skip_soamode(m) ) cycle
1032 13462573460 : sat(m,ll) = g0_soa(ll)/max( a_ooa_sum_tmp(m), a_min1 )
1033 13462573460 : g_star(m,ll) = sat(m,ll)*a_soa(m,ll)
1034 13462573460 : phi(m,ll) = (g_soa(ll) - g_star(m,ll))/max( g_soa(ll), g_star(m,ll), g_min1 )
1035 27598275593 : tmpb = tmpb + xferrate(m,ll)*abs(phi(m,ll))
1036 : end do
1037 1346257346 : tmpa = max( tmpa, tmpb )
1038 : end do
1039 :
1040 : if (dtsub_fixed > 0.0_r8) then
1041 : dtcur = dtsub_fixed
1042 : tcur = tcur + dtcur
1043 : else
1044 673128673 : dtmax = dtfull-tcur
1045 673128673 : if (dtmax*tmpa <= alpha) then
1046 : ! here alpha/tmpa >= dtmax, so this is final substep
1047 : dtcur = dtmax
1048 : tcur = dtfull
1049 : else
1050 599585915 : dtcur = alpha/tmpa
1051 599585915 : tcur = tcur + dtcur
1052 : end if
1053 : end if
1054 :
1055 : ! step 1 - for modes where soa is condensing, estimate "new" a_soa(m,ll)
1056 : ! using an explicit calculation with "old" g_soa
1057 : ! and g_star(m,ll) calculated using "old" a_soa(m,ll)
1058 : ! do this to get better estimate of "new" a_soa(m,ll) and sat(m,ll)
1059 27598275593 : do m = 1, nbins
1060 26925146920 : if ( skip_soamode(m) ) cycle
1061 26925146920 : do ll = 1, ntot_soaspec
1062 : ! first ll loop calcs a_soa_tmp(m,ll) & a_ooa_sum_tmp
1063 13462573460 : a_soa_tmp(m,ll) = a_soa(m,ll)
1064 13462573460 : beta(m,ll) = dtcur*xferrate(m,ll)
1065 13462573460 : del_g_soa_tmp(ll) = g_soa(ll) - g_star(m,ll)
1066 26925146920 : if (del_g_soa_tmp(ll) > 0.0_r8) then
1067 12983394789 : a_soa_tmp(m,ll) = a_soa(m,ll) + beta(m,ll)*del_g_soa_tmp(ll)
1068 : end if
1069 : end do
1070 26925146920 : a_ooa_sum_tmp(m) = sum( a_soa_tmp(m,1:ntot_soaspec) )
1071 27598275593 : do ll = 1, ntot_soaspec
1072 : ! second ll loop calcs sat & g_star
1073 40387720380 : if (del_g_soa_tmp(ll) > 0.0_r8) then
1074 12983394789 : sat(m,ll) = g0_soa(ll)/max( a_ooa_sum_tmp(m), a_min1 )
1075 12983394789 : g_star(m,ll) = sat(m,ll)*a_soa_tmp(m,ll) ! this just needed for diagnostics
1076 : end if
1077 : end do
1078 : end do
1079 :
1080 : ! step 2 - implicit in g_soa and semi-implicit in a_soa,
1081 : ! with g_star(m,ll) calculated semi-implicitly
1082 1419801026 : do ll = 1, ntot_soaspec
1083 : tmpa = 0.0_r8
1084 : tmpb = 0.0_r8
1085 27598275593 : do m = 1, nbins
1086 26925146920 : if ( skip_soamode(m) ) cycle
1087 13462573460 : tmpa = tmpa + a_soa(m,ll)/(1.0_r8 + beta(m,ll)*sat(m,ll))
1088 27598275593 : tmpb = tmpb + beta(m,ll)/(1.0_r8 + beta(m,ll)*sat(m,ll))
1089 : end do
1090 :
1091 673128673 : g_soa(ll) = (tot_soa(ll) - tmpa)/(1.0_r8 + tmpb)
1092 673128673 : g_soa(ll) = max( 0.0_r8, g_soa(ll) )
1093 28271404266 : do m = 1, nbins
1094 26925146920 : if ( skip_soamode(m) ) cycle
1095 : a_soa(m,ll) = (a_soa(m,ll) + beta(m,ll)*g_soa(ll))/ &
1096 27598275593 : (1.0_r8 + beta(m,ll)*sat(m,ll))
1097 : end do
1098 : end do
1099 :
1100 : end do time_loop
1101 :
1102 : ! calculate outgoing tendencies (at the host-code molec. weight)
1103 : ! (a_soa & g_soa are at actual mw, but a_soa_in & g_soa_in are at host-code mw)
1104 147087360 : do ll = 1, ntot_soaspec
1105 73543680 : tmpf = mw_soa/mw_soa_host
1106 73543680 : g_soa_tend(ll) = (g_soa(ll)*tmpf - g_soa_in(ll))/dtfull
1107 3088834560 : do m = 1, nbins
1108 2941747200 : if ( skip_soamode(m) ) cycle
1109 3015290880 : a_soa_tend(m,ll) = (a_soa(m,ll)*tmpf - a_soa_in(m,ll))/dtfull
1110 : end do
1111 : end do
1112 :
1113 73543680 : end subroutine carma_aero_soaexch
1114 :
1115 : !----------------------------------------------------------------------
1116 :
1117 : end module carma_aero_gasaerexch
|
