Line data Source code
1 : module modal_aero_wateruptake
2 :
3 : ! RCE 07.04.13: Adapted from MIRAGE2 code
4 :
5 : use shr_kind_mod, only: r8 => shr_kind_r8
6 : use physconst, only: pi, rhoh2o, rair
7 : use ppgrid, only: pcols, pver
8 : use physics_types, only: physics_state
9 : use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_old_tim_idx, pbuf_get_field
10 :
11 : use wv_saturation, only: qsat_water
12 : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_mmr, rad_cnst_get_aer_props, &
13 : rad_cnst_get_mode_props
14 : use cam_history, only: addfld, add_default, outfld, horiz_only
15 : use cam_logfile, only: iulog
16 : use ref_pres, only: top_lev => clim_modal_aero_top_lev
17 : use phys_control, only: phys_getopts
18 : use cam_abortutils, only: endrun
19 :
20 : implicit none
21 : private
22 : save
23 :
24 : public :: &
25 : modal_aero_wateruptake_init, &
26 : modal_aero_wateruptake_dr, &
27 : modal_aero_wateruptake_sub, &
28 : modal_aero_kohler
29 :
30 : public :: modal_aero_wateruptake_reg
31 :
32 : real(r8), parameter :: third = 1._r8/3._r8
33 : real(r8), parameter :: pi43 = pi*4.0_r8/3.0_r8
34 :
35 :
36 : ! Physics buffer indices
37 : integer :: cld_idx = 0
38 : integer :: dgnum_idx = 0
39 : integer :: dgnumwet_idx = 0
40 : integer :: sulfeq_idx = 0
41 : integer :: wetdens_ap_idx = 0
42 : integer :: qaerwat_idx = 0
43 : integer :: hygro_idx = 0
44 : integer :: dryvol_idx = 0
45 : integer :: dryrad_idx = 0
46 : integer :: drymass_idx = 0
47 : integer :: so4dryvol_idx = 0
48 : integer :: naer_idx = 0
49 :
50 :
51 : logical, public :: modal_strat_sulfate = .false. ! If .true. then MAM sulfate surface area density used in stratospheric heterogeneous chemistry
52 :
53 : !===============================================================================
54 : contains
55 : !===============================================================================
56 :
57 1536 : subroutine modal_aero_wateruptake_reg()
58 :
59 : use physics_buffer, only: pbuf_add_field, dtype_r8
60 : use rad_constituents, only: rad_cnst_get_info
61 :
62 : integer :: nmodes
63 :
64 1536 : call rad_cnst_get_info(0, nmodes=nmodes)
65 6144 : call pbuf_add_field('DGNUMWET', 'global', dtype_r8, (/pcols, pver, nmodes/), dgnumwet_idx)
66 6144 : call pbuf_add_field('WETDENS_AP', 'physpkg', dtype_r8, (/pcols, pver, nmodes/), wetdens_ap_idx)
67 :
68 : ! 1st order rate for direct conversion of strat. cloud water to precip (1/s)
69 6144 : call pbuf_add_field('QAERWAT', 'physpkg', dtype_r8, (/pcols, pver, nmodes/), qaerwat_idx)
70 :
71 1536 : if (modal_strat_sulfate) then
72 0 : call pbuf_add_field('MAMH2SO4EQ', 'global', dtype_r8, (/pcols, pver, nmodes/), sulfeq_idx)
73 : end if
74 :
75 :
76 1536 : end subroutine modal_aero_wateruptake_reg
77 :
78 : !===============================================================================
79 : !===============================================================================
80 :
81 1536 : subroutine modal_aero_wateruptake_init(pbuf2d)
82 1536 : use time_manager, only: is_first_step
83 : use physics_buffer,only: pbuf_set_field
84 : use infnan, only : nan, assignment(=)
85 :
86 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
87 : real(r8) :: real_nan
88 :
89 : integer :: m, nmodes
90 : logical :: history_aerosol ! Output the MAM aerosol variables and tendencies
91 :
92 : character(len=3) :: trnum ! used to hold mode number (as characters)
93 : !----------------------------------------------------------------------------
94 :
95 1536 : real_nan = nan
96 :
97 1536 : cld_idx = pbuf_get_index('CLD')
98 1536 : dgnum_idx = pbuf_get_index('DGNUM')
99 :
100 1536 : hygro_idx = pbuf_get_index('HYGRO')
101 1536 : dryvol_idx = pbuf_get_index('DRYVOL')
102 1536 : dryrad_idx = pbuf_get_index('DRYRAD')
103 1536 : drymass_idx = pbuf_get_index('DRYMASS')
104 1536 : so4dryvol_idx = pbuf_get_index('SO4DRYVOL')
105 1536 : naer_idx = pbuf_get_index('NAER')
106 :
107 : ! assume for now that will compute wateruptake for climate list modes only
108 :
109 1536 : call rad_cnst_get_info(0, nmodes=nmodes)
110 :
111 7680 : do m = 1, nmodes
112 6144 : write(trnum, '(i3.3)') m
113 : call addfld('dgnd_a'//trnum(2:3), (/ 'lev' /), 'A', 'm', &
114 12288 : 'dry dgnum, interstitial, mode '//trnum(2:3))
115 : call addfld('dgnw_a'//trnum(2:3), (/ 'lev' /), 'A', 'm', &
116 12288 : 'wet dgnum, interstitial, mode '//trnum(2:3))
117 : call addfld('wat_a'//trnum(3:3), (/ 'lev' /), 'A', 'm', &
118 12288 : 'aerosol water, interstitial, mode '//trnum(2:3))
119 :
120 : ! determine default variables
121 6144 : call phys_getopts(history_aerosol_out = history_aerosol)
122 :
123 7680 : if (history_aerosol) then
124 0 : call add_default('dgnd_a'//trnum(2:3), 1, ' ')
125 0 : call add_default('dgnw_a'//trnum(2:3), 1, ' ')
126 0 : call add_default('wat_a'//trnum(3:3), 1, ' ')
127 : endif
128 :
129 : end do
130 :
131 3072 : call addfld('PM25', (/ 'lev' /), 'A', 'kg/m3', 'PM2.5 concentration')
132 1536 : call addfld('PM25_SRF', horiz_only, 'A', 'kg/m3', 'surface PM2.5 concentration')
133 :
134 1536 : if (is_first_step()) then
135 : ! initialize fields in physics buffer
136 768 : call pbuf_set_field(pbuf2d, dgnumwet_idx, 0.0_r8)
137 768 : if (modal_strat_sulfate) then
138 : ! initialize fields in physics buffer to NaN (not a number)
139 : ! so model will crash if used before initialization
140 0 : call pbuf_set_field(pbuf2d, sulfeq_idx, real_nan)
141 : endif
142 : endif
143 :
144 1536 : end subroutine modal_aero_wateruptake_init
145 :
146 : !===============================================================================
147 :
148 :
149 61920 : subroutine modal_aero_wateruptake_dr(state, pbuf, list_idx_in, dgnumdry_m, dgnumwet_m, &
150 : qaerwat_m, wetdens_m, hygro_m, dryvol_m, dryrad_m, drymass_m,&
151 : so4dryvol_m, naer_m)
152 : !-----------------------------------------------------------------------
153 : !
154 : ! CAM specific driver for modal aerosol water uptake code.
155 : !
156 : ! *** N.B. *** The calculation has been enabled for diagnostic mode lists
157 : ! via optional arguments. If the list_idx arg is present then
158 : ! all the optional args must be present.
159 : !
160 : !-----------------------------------------------------------------------
161 :
162 1536 : use time_manager, only: is_first_step
163 : use cam_history, only: outfld, fieldname_len
164 : use tropopause, only: tropopause_find, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE
165 : ! Arguments
166 : type(physics_state), target, intent(in) :: state ! Physics state variables
167 : type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer
168 :
169 : integer, optional, intent(in) :: list_idx_in
170 : real(r8), optional, pointer :: dgnumdry_m(:,:,:)
171 : real(r8), optional, pointer :: dgnumwet_m(:,:,:)
172 : real(r8), optional, pointer :: qaerwat_m(:,:,:)
173 : real(r8), optional, pointer :: wetdens_m(:,:,:)
174 : real(r8), optional, pointer :: hygro_m(:,:,:)
175 : real(r8), optional, pointer :: dryvol_m(:,:,:)
176 : real(r8), optional, pointer :: dryrad_m(:,:,:)
177 : real(r8), optional, pointer :: drymass_m(:,:,:)
178 : real(r8), optional, pointer :: so4dryvol_m(:,:,:)
179 : real(r8), optional, pointer :: naer_m(:,:,:)
180 :
181 : ! local variables
182 :
183 : integer :: lchnk ! chunk index
184 : integer :: ncol ! number of columns
185 : integer :: list_idx ! radiative constituents list index
186 : integer :: stat
187 :
188 : integer :: i, k, l, m
189 : integer :: itim_old
190 : integer :: nmodes
191 : integer :: nspec
192 : integer :: tropLev(pcols)
193 :
194 : character(len=fieldname_len+3) :: fieldname
195 :
196 61920 : real(r8), pointer :: h2ommr(:,:) ! specific humidity
197 61920 : real(r8), pointer :: t(:,:) ! temperatures (K)
198 61920 : real(r8), pointer :: pmid(:,:) ! layer pressure (Pa)
199 :
200 61920 : real(r8), pointer :: cldn(:,:) ! layer cloud fraction (0-1)
201 61920 : real(r8), pointer :: dgncur_a(:,:,:)
202 61920 : real(r8), pointer :: dgncur_awet(:,:,:)
203 61920 : real(r8), pointer :: wetdens(:,:,:)
204 61920 : real(r8), pointer :: qaerwat(:,:,:)
205 :
206 61920 : real(r8), pointer :: raer(:,:) ! aerosol species MRs (kg/kg)
207 61920 : real(r8), pointer :: maer(:,:,:) ! accumulated aerosol mode MRs
208 61920 : real(r8), pointer :: hygro(:,:,:) ! volume-weighted mean hygroscopicity (--)
209 61920 : real(r8), pointer :: naer(:,:,:) ! aerosol number MR (bounded!) (#/kg-air)
210 61920 : real(r8), pointer :: dryvol(:,:,:) ! single-particle-mean dry volume (m3)
211 61920 : real(r8), pointer :: so4dryvol(:,:,:) ! single-particle-mean so4 dry volume (m3)
212 61920 : real(r8), pointer :: drymass(:,:,:) ! single-particle-mean dry mass (kg)
213 61920 : real(r8), pointer :: dryrad(:,:,:) ! dry volume mean radius of aerosol (m)
214 :
215 61920 : real(r8), allocatable :: wetrad(:,:,:) ! wet radius of aerosol (m)
216 61920 : real(r8), allocatable :: wetvol(:,:,:) ! single-particle-mean wet volume (m3)
217 61920 : real(r8), allocatable :: wtrvol(:,:,:) ! single-particle-mean water volume in wet aerosol (m3)
218 :
219 61920 : real(r8), allocatable :: rhcrystal(:)
220 61920 : real(r8), allocatable :: rhdeliques(:)
221 61920 : real(r8), allocatable :: specdens_1(:)
222 :
223 61920 : real(r8), pointer :: sulfeq(:,:,:) ! H2SO4 equilibrium mixing ratios over particles (mol/mol)
224 61920 : real(r8), allocatable :: wtpct(:,:,:) ! sulfate aerosol composition, weight % H2SO4
225 61920 : real(r8), allocatable :: sulden(:,:,:) ! sulfate aerosol mass density (g/cm3)
226 :
227 : real(r8) :: specdens, so4specdens
228 : real(r8) :: sigmag
229 61920 : real(r8), allocatable :: alnsg(:)
230 : real(r8) :: rh(pcols,pver) ! relative humidity (0-1)
231 : real(r8) :: dmean, qh2so4_equilib, wtpct_mode, sulden_mode
232 :
233 : real(r8) :: es(pcols) ! saturation vapor pressure
234 : real(r8) :: qs(pcols) ! saturation specific humidity
235 :
236 : real(r8) :: pm25(pcols,pver) ! PM2.5 diagnostics
237 : real(r8) :: rhoair(pcols,pver)
238 :
239 : character(len=3) :: trnum ! used to hold mode number (as characters)
240 : character(len=32) :: spectype
241 : !-----------------------------------------------------------------------
242 :
243 61920 : lchnk = state%lchnk
244 61920 : ncol = state%ncol
245 :
246 61920 : list_idx = 0
247 61920 : if (present(list_idx_in)) then
248 0 : list_idx = list_idx_in
249 :
250 : ! check that all optional args are present
251 : if (.not. present(dgnumdry_m) .or. .not. present(dgnumwet_m) .or. &
252 0 : .not. present(qaerwat_m) .or. .not. present(wetdens_m)) then
253 : call endrun('modal_aero_wateruptake_dr called for'// &
254 0 : 'diagnostic list but required args not present')
255 : end if
256 :
257 : ! arrays for diagnostic calculations must be associated
258 : if (.not. associated(dgnumdry_m) .or. .not. associated(dgnumwet_m) .or. &
259 0 : .not. associated(qaerwat_m) .or. .not. associated(wetdens_m)) then
260 : call endrun('modal_aero_wateruptake_dr called for'// &
261 0 : 'diagnostic list but required args not associated')
262 : end if
263 :
264 0 : if (modal_strat_sulfate) then
265 : call endrun('modal_aero_wateruptake_dr cannot be called with optional arguments and'// &
266 0 : ' having modal_strat_sulfate set to true')
267 : end if
268 : end if
269 :
270 : ! loop over all aerosol modes
271 61920 : call rad_cnst_get_info(list_idx, nmodes=nmodes)
272 :
273 61920 : if (modal_strat_sulfate) then
274 0 : call pbuf_get_field(pbuf, sulfeq_idx, sulfeq )
275 : endif
276 :
277 : allocate( &
278 0 : wetrad(pcols,pver,nmodes), &
279 0 : wetvol(pcols,pver,nmodes), &
280 0 : wtrvol(pcols,pver,nmodes), &
281 0 : wtpct(pcols,pver,nmodes), &
282 0 : sulden(pcols,pver,nmodes), &
283 0 : rhcrystal(nmodes), &
284 0 : rhdeliques(nmodes), &
285 : specdens_1(nmodes), &
286 743040 : alnsg(nmodes) )
287 :
288 391891680 : wtpct(:,:,:) = 75._r8
289 391891680 : sulden(:,:,:) = 1.923_r8
290 :
291 61920 : if (list_idx == 0) then
292 61920 : call pbuf_get_field(pbuf, dgnum_idx, dgncur_a )
293 61920 : call pbuf_get_field(pbuf, dgnumwet_idx, dgncur_awet )
294 61920 : call pbuf_get_field(pbuf, wetdens_ap_idx, wetdens)
295 61920 : call pbuf_get_field(pbuf, qaerwat_idx, qaerwat)
296 61920 : call pbuf_get_field(pbuf, hygro_idx, hygro)
297 61920 : call pbuf_get_field(pbuf, dryvol_idx, dryvol)
298 61920 : call pbuf_get_field(pbuf, dryrad_idx, dryrad)
299 61920 : call pbuf_get_field(pbuf, drymass_idx, drymass)
300 61920 : call pbuf_get_field(pbuf, so4dryvol_idx, so4dryvol)
301 61920 : call pbuf_get_field(pbuf, naer_idx, naer)
302 :
303 61920 : if (is_first_step()) then
304 78372144 : dgncur_awet(:,:,:) = dgncur_a(:,:,:)
305 : end if
306 : else
307 0 : dgncur_a => dgnumdry_m
308 0 : dgncur_awet => dgnumwet_m
309 0 : qaerwat => qaerwat_m
310 0 : wetdens => wetdens_m
311 0 : hygro => hygro_m
312 0 : dryvol => dryvol_m
313 0 : dryrad => dryrad_m
314 0 : drymass => drymass_m
315 0 : so4dryvol => so4dryvol_m
316 0 : naer => naer_m
317 : end if
318 :
319 61920 : if (modal_strat_sulfate) then
320 : ! get tropopause level
321 0 : call tropopause_find(state, tropLev, primary=TROP_ALG_HYBSTOB, backup=TROP_ALG_CLIMATE)
322 : endif
323 :
324 61920 : h2ommr => state%q(:,:,1)
325 61920 : t => state%t
326 61920 : pmid => state%pmid
327 :
328 185760 : allocate( maer(pcols,pver,nmodes))
329 391891680 : maer(:,:,:) = 0._r8
330 :
331 309600 : do m = 1, nmodes
332 :
333 : call rad_cnst_get_mode_props(list_idx, m, sigmag=sigmag, &
334 247680 : rhcrystal=rhcrystal(m), rhdeliques=rhdeliques(m))
335 :
336 : ! get mode info
337 247680 : call rad_cnst_get_info(list_idx, m, nspec=nspec)
338 :
339 1176480 : do l = 1, nspec
340 :
341 : ! accumulate the aerosol masses of each mode
342 928800 : call rad_cnst_get_aer_mmr(0, m, l, 'a', state, pbuf, raer)
343 2885565600 : maer(:ncol,:,m)= maer(:ncol,:,m) + raer(:ncol,:)
344 :
345 : ! get species interstitial mixing ratio ('a')
346 : call rad_cnst_get_aer_props(list_idx, m, l, density_aer=specdens, &
347 928800 : spectype=spectype)
348 :
349 928800 : if (modal_strat_sulfate .and. (trim(spectype).eq.'sulfate')) then
350 0 : so4specdens=specdens
351 : end if
352 :
353 2105280 : if (l == 1) then
354 : ! save off these values to be used as defaults
355 247680 : specdens_1(m) = specdens
356 : end if
357 :
358 : end do
359 :
360 247680 : alnsg(m) = log(sigmag)
361 :
362 557280 : if (modal_strat_sulfate) then
363 0 : do k = top_lev, pver
364 0 : do i = 1, ncol
365 0 : dmean = dgncur_awet(i,k,m)*exp(1.5_r8*alnsg(m)**2)
366 0 : call calc_h2so4_equilib_mixrat( t(i,k), pmid(i,k), h2ommr(i,k), dmean, &
367 0 : qh2so4_equilib, wtpct_mode, sulden_mode )
368 0 : sulfeq(i,k,m) = qh2so4_equilib
369 0 : wtpct(i,k,m) = wtpct_mode
370 0 : sulden(i,k,m) = sulden_mode
371 : end do ! i = 1, ncol
372 : end do ! k = top_lev, pver
373 :
374 0 : fieldname = ' '
375 0 : write(fieldname,fmt='(a,i1)') 'wtpct_a',m
376 0 : call outfld(fieldname,wtpct(1:ncol,1:pver,m), ncol, lchnk )
377 :
378 0 : fieldname = ' '
379 0 : write(fieldname,fmt='(a,i1)') 'sulfeq_a',m
380 0 : call outfld(fieldname,sulfeq(1:ncol,1:pver,m), ncol, lchnk )
381 :
382 0 : fieldname = ' '
383 0 : write(fieldname,fmt='(a,i1)') 'sulden_a',m
384 0 : call outfld(fieldname,sulden(1:ncol,1:pver,m), ncol, lchnk )
385 :
386 : end if
387 :
388 : end do ! m = 1, nmodes
389 :
390 : ! relative humidity calc
391 :
392 61920 : itim_old = pbuf_old_tim_idx()
393 247680 : call pbuf_get_field(pbuf, cld_idx, cldn, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) )
394 :
395 5820480 : do k = top_lev, pver
396 5758560 : call qsat_water(t(1:ncol,k), pmid(1:ncol,k), es(1:ncol), qs(1:ncol), ncol)
397 96216480 : do i = 1, ncol
398 90396000 : if (qs(i) > h2ommr(i,k)) then
399 88785055 : rh(i,k) = h2ommr(i,k)/qs(i)
400 : else
401 1610945 : rh(i,k) = 0.98_r8
402 : endif
403 90396000 : rh(i,k) = max(rh(i,k), 0.0_r8)
404 90396000 : rh(i,k) = min(rh(i,k), 0.98_r8)
405 90396000 : if (cldn(i,k) .lt. 1.0_r8) then
406 89676082 : rh(i,k) = (rh(i,k) - cldn(i,k)) / (1.0_r8 - cldn(i,k)) ! clear portion
407 : end if
408 90396000 : rh(i,k) = max(rh(i,k), 0.0_r8)
409 96154560 : rhoair(i,k) = pmid(i,k)/(rair*t(i,k))
410 : end do
411 : end do
412 :
413 : call modal_aero_wateruptake_sub( &
414 : ncol, nmodes, rhcrystal, rhdeliques, dryrad, &
415 : hygro, rh, dryvol, so4dryvol, so4specdens, tropLev, &
416 61920 : wetrad, wetvol, wtrvol, sulden, wtpct)
417 :
418 391891680 : qaerwat = 0.0_r8
419 :
420 309600 : do m = 1, nmodes
421 :
422 23343840 : do k = top_lev, pver
423 384865920 : do i = 1, ncol
424 :
425 361584000 : dgncur_awet(i,k,m) = dgncur_a(i,k,m) * (wetrad(i,k,m)/dryrad(i,k,m))
426 361584000 : qaerwat(i,k,m) = rhoh2o*naer(i,k,m)*wtrvol(i,k,m)
427 :
428 : ! compute aerosol wet density (kg/m3)
429 384618240 : if (wetvol(i,k,m) > 1.0e-30_r8) then
430 361584000 : wetdens(i,k,m) = (drymass(i,k,m) + rhoh2o*wtrvol(i,k,m))/wetvol(i,k,m)
431 : else
432 0 : wetdens(i,k,m) = specdens_1(m)
433 : end if
434 : end do
435 : end do
436 :
437 : end do ! modes
438 :
439 61920 : if (list_idx == 0) then
440 :
441 61920 : pm25(:,:)=0._r8
442 :
443 309600 : do m = 1, nmodes
444 : ! output to history
445 247680 : write( trnum, '(i3.3)' ) m
446 247680 : call outfld( 'wat_a'//trnum(3:3), qaerwat(:,:,m), pcols, lchnk)
447 247680 : call outfld( 'dgnd_a'//trnum(2:3), dgncur_a(:,:,m), pcols, lchnk)
448 247680 : call outfld( 'dgnw_a'//trnum(2:3), dgncur_awet(:,:,m), pcols, lchnk)
449 :
450 : ! calculate PM2.5 diagnostics -- dgncur_a is zero above top_lev
451 23343840 : do k = top_lev, pver
452 384865920 : do i=1,ncol
453 1084752000 : pm25(i,k) = pm25(i,k)+maer(i,k,m)*(1._r8-(0.5_r8 - 0.5_r8*erf(log(2.5e-6_r8/dgncur_a(i,k,m))/ &
454 1469370240 : (2._r8**0.5_r8*alnsg(m)))))*rhoair(i,k)
455 : end do
456 : end do
457 : end do
458 :
459 61920 : call outfld('PM25', pm25(:,:), pcols, lchnk)
460 61920 : call outfld('PM25_SRF', pm25(:,pver), pcols, lchnk)
461 :
462 : end if
463 :
464 61920 : deallocate(maer, alnsg)
465 0 : deallocate( &
466 61920 : wetrad, wetvol, wtrvol, wtpct, sulden, rhcrystal, rhdeliques, specdens_1 )
467 123840 : end subroutine modal_aero_wateruptake_dr
468 :
469 : !===============================================================================
470 :
471 61920 : subroutine modal_aero_wateruptake_sub( &
472 61920 : ncol, nmodes, rhcrystal, rhdeliques, dryrad, &
473 185760 : hygro, rh, dryvol, so4dryvol, so4specdens, troplev, &
474 123840 : wetrad, wetvol, wtrvol, sulden, wtpct)
475 :
476 : !-----------------------------------------------------------------------
477 : !
478 : ! Purpose: Compute aerosol wet radius
479 : !
480 : ! Method: Kohler theory
481 : !
482 : ! Author: S. Ghan
483 : !
484 : !-----------------------------------------------------------------------
485 :
486 :
487 : ! Arguments
488 : integer, intent(in) :: ncol ! number of columns
489 : integer, intent(in) :: nmodes
490 : integer, intent(in) :: troplev(:)
491 :
492 : real(r8), intent(in) :: rhcrystal(:)
493 : real(r8), intent(in) :: rhdeliques(:)
494 : real(r8), intent(in) :: dryrad(:,:,:) ! dry volume mean radius of aerosol (m)
495 : real(r8), intent(in) :: hygro(:,:,:) ! volume-weighted mean hygroscopicity (--)
496 : real(r8), intent(in) :: rh(:,:) ! relative humidity (0-1)
497 : real(r8), intent(in) :: dryvol(:,:,:) ! dry volume of single aerosol (m3)
498 : real(r8), intent(in) :: so4dryvol(:,:,:) ! dry volume of sulfate in single aerosol (m3)
499 : real(r8), intent(in) :: so4specdens ! mass density sulfate in single aerosol (kg/m3)
500 : real(r8), intent(in) :: wtpct(:,:,:) ! sulfate aerosol composition, weight % H2SO4
501 : real(r8), intent(in) :: sulden(:,:,:) ! sulfate aerosol mass density (g/cm3)
502 :
503 : real(r8), intent(out) :: wetrad(:,:,:) ! wet radius of aerosol (m)
504 : real(r8), intent(out) :: wetvol(:,:,:) ! single-particle-mean wet volume (m3)
505 : real(r8), intent(out) :: wtrvol(:,:,:) ! single-particle-mean water volume in wet aerosol (m3)
506 :
507 : ! local variables
508 :
509 : integer :: i, k, m
510 :
511 : real(r8) :: hystfac ! working variable for hysteresis
512 : !-----------------------------------------------------------------------
513 :
514 :
515 : ! loop over all aerosol modes
516 309600 : do m = 1, nmodes
517 :
518 247680 : hystfac = 1.0_r8 / max(1.0e-5_r8, (rhdeliques(m) - rhcrystal(m)))
519 :
520 23343840 : do k = top_lev, pver
521 384865920 : do i = 1, ncol
522 :
523 384618240 : if ( modal_strat_sulfate .and. (k<troplev(i))) then
524 0 : wetvol(i,k,m) = dryvol(i,k,m)-so4dryvol(i,k,m)
525 0 : wetvol(i,k,m) = wetvol(i,k,m)+so4dryvol(i,k,m)*so4specdens/sulden(i,k,m)/wtpct(i,k,m)/10._r8
526 0 : wetvol(i,k,m) = max(wetvol(i,k,m), dryvol(i,k,m))
527 0 : wetrad(i,k,m) = (wetvol(i,k,m)/pi43)**third
528 0 : wetrad(i,k,m) = max(wetrad(i,k,m), dryrad(i,k,m))
529 0 : wtrvol(i,k,m) = wetvol(i,k,m) - dryvol(i,k,m)
530 0 : wtrvol(i,k,m) = max(wtrvol(i,k,m), 0.0_r8)
531 : else
532 : ! compute wet radius for each mode
533 361584000 : call modal_aero_kohler(dryrad(i:i,k,m), hygro(i:i,k,m), rh(i:i,k), wetrad(i:i,k,m), 1)
534 :
535 361584000 : wetrad(i,k,m) = max(wetrad(i,k,m), dryrad(i,k,m))
536 361584000 : wetvol(i,k,m) = pi43*wetrad(i,k,m)**3
537 361584000 : wetvol(i,k,m) = max(wetvol(i,k,m), dryvol(i,k,m))
538 361584000 : wtrvol(i,k,m) = wetvol(i,k,m) - dryvol(i,k,m)
539 361584000 : wtrvol(i,k,m) = max(wtrvol(i,k,m), 0.0_r8)
540 :
541 : ! apply simple treatment of deliquesence/crystallization hysteresis
542 : ! for rhcrystal < rh < rhdeliques, aerosol water is a fraction of
543 : ! the "upper curve" value, and the fraction is a linear function of rh
544 361584000 : if (rh(i,k) < rhcrystal(m)) then
545 278012608 : wetrad(i,k,m) = dryrad(i,k,m)
546 278012608 : wetvol(i,k,m) = dryvol(i,k,m)
547 278012608 : wtrvol(i,k,m) = 0.0_r8
548 83571392 : else if (rh(i,k) < rhdeliques(m)) then
549 50497732 : wtrvol(i,k,m) = wtrvol(i,k,m)*hystfac*(rh(i,k) - rhcrystal(m))
550 50497732 : wtrvol(i,k,m) = max(wtrvol(i,k,m), 0.0_r8)
551 50497732 : wetvol(i,k,m) = dryvol(i,k,m) + wtrvol(i,k,m)
552 50497732 : wetrad(i,k,m) = (wetvol(i,k,m)/pi43)**third
553 : end if
554 : end if
555 :
556 : end do ! columns
557 : end do ! levels
558 :
559 : end do ! modes
560 :
561 61920 : end subroutine modal_aero_wateruptake_sub
562 :
563 : !-----------------------------------------------------------------------
564 361584000 : subroutine modal_aero_kohler( &
565 361584000 : rdry_in, hygro, s, rwet_out, im )
566 :
567 : ! calculates equlibrium radius r of haze droplets as function of
568 : ! dry particle mass and relative humidity s using kohler solution
569 : ! given in pruppacher and klett (eqn 6-35)
570 :
571 : ! for multiple aerosol types, assumes an internal mixture of aerosols
572 :
573 : implicit none
574 :
575 : ! arguments
576 : integer :: im ! number of grid points to be processed
577 : real(r8) :: rdry_in(:) ! aerosol dry radius (m)
578 : real(r8) :: hygro(:) ! aerosol volume-mean hygroscopicity (--)
579 : real(r8) :: s(:) ! relative humidity (1 = saturated)
580 : real(r8) :: rwet_out(:) ! aerosol wet radius (m)
581 :
582 : ! local variables
583 : integer, parameter :: imax=200
584 : integer :: i, n, nsol
585 :
586 : real(r8) :: a, b
587 : real(r8) :: p40(imax),p41(imax),p42(imax),p43(imax) ! coefficients of polynomial
588 : real(r8) :: p30(imax),p31(imax),p32(imax) ! coefficients of polynomial
589 : real(r8) :: p
590 : real(r8) :: r3, r4
591 723168000 : real(r8) :: r(im) ! wet radius (microns)
592 : real(r8) :: rdry(imax) ! radius of dry particle (microns)
593 : real(r8) :: ss ! relative humidity (1 = saturated)
594 : real(r8) :: slog(imax) ! log relative humidity
595 : real(r8) :: vol(imax) ! total volume of particle (microns**3)
596 : real(r8) :: xi, xr
597 :
598 : complex(r8) :: cx4(4,imax),cx3(3,imax)
599 :
600 : real(r8), parameter :: eps = 1.e-4_r8
601 : real(r8), parameter :: mw = 18._r8
602 : real(r8), parameter :: pi = 3.14159_r8
603 : real(r8), parameter :: rhow = 1._r8
604 : real(r8), parameter :: surften = 76._r8
605 : real(r8), parameter :: tair = 273._r8
606 : real(r8), parameter :: third = 1._r8/3._r8
607 : real(r8), parameter :: ugascon = 8.3e7_r8
608 :
609 :
610 : ! effect of organics on surface tension is neglected
611 361584000 : a=2.e4_r8*mw*surften/(ugascon*tair*rhow)
612 :
613 723168000 : do i=1,im
614 361584000 : rdry(i) = rdry_in(i)*1.0e6_r8 ! convert (m) to (microns)
615 361584000 : vol(i) = rdry(i)**3 ! vol is r**3, not volume
616 361584000 : b = vol(i)*hygro(i)
617 :
618 : ! quartic
619 361584000 : ss=min(s(i),1._r8-eps)
620 361584000 : ss=max(ss,1.e-10_r8)
621 361584000 : slog(i)=log(ss)
622 361584000 : p43(i)=-a/slog(i)
623 361584000 : p42(i)=0._r8
624 361584000 : p41(i)=b/slog(i)-vol(i)
625 361584000 : p40(i)=a*vol(i)/slog(i)
626 : ! cubic for rh=1
627 361584000 : p32(i)=0._r8
628 361584000 : p31(i)=-b/a
629 723168000 : p30(i)=-vol(i)
630 : end do
631 :
632 :
633 723168000 : do 100 i=1,im
634 :
635 : ! if(vol(i).le.1.e-20)then
636 361584000 : if(vol(i).le.1.e-12_r8)then
637 0 : r(i)=rdry(i)
638 0 : go to 100
639 : endif
640 :
641 361584000 : p=abs(p31(i))/(rdry(i)*rdry(i))
642 361584000 : if(p.lt.eps)then
643 : ! approximate solution for small particles
644 90443636 : r(i)=rdry(i)*(1._r8+p*third/(1._r8-slog(i)*rdry(i)/a))
645 : else
646 271140364 : call makoh_quartic(cx4(1,i),p43(i),p42(i),p41(i),p40(i),1)
647 : ! find smallest real(r8) solution
648 271140364 : r(i)=1000._r8*rdry(i)
649 271140364 : nsol=0
650 1355701820 : do n=1,4
651 1084561456 : xr=real(cx4(n,i))
652 1084561456 : xi=aimag(cx4(n,i))
653 1084561456 : if(abs(xi).gt.abs(xr)*eps) cycle
654 542488700 : if(xr.gt.r(i)) cycle
655 542488700 : if(xr.lt.rdry(i)*(1._r8-eps)) cycle
656 271140364 : if(xr.ne.xr) cycle
657 271140364 : r(i)=xr
658 1355701820 : nsol=n
659 : end do
660 271140364 : if(nsol.eq.0)then
661 : write(iulog,*) &
662 0 : 'ccm kohlerc - no real(r8) solution found (quartic)'
663 0 : write(iulog,*)'roots =', (cx4(n,i),n=1,4)
664 0 : write(iulog,*)'p0-p3 =', p40(i), p41(i), p42(i), p43(i)
665 0 : write(iulog,*)'rh=',s(i)
666 0 : write(iulog,*)'setting radius to dry radius=',rdry(i)
667 0 : r(i)=rdry(i)
668 : ! stop
669 : endif
670 : endif
671 :
672 361584000 : if(s(i).gt.1._r8-eps)then
673 : ! save quartic solution at s=1-eps
674 0 : r4=r(i)
675 : ! cubic for rh=1
676 0 : p=abs(p31(i))/(rdry(i)*rdry(i))
677 0 : if(p.lt.eps)then
678 0 : r(i)=rdry(i)*(1._r8+p*third)
679 : else
680 0 : call makoh_cubic(cx3,p32,p31,p30,im)
681 : ! find smallest real(r8) solution
682 0 : r(i)=1000._r8*rdry(i)
683 0 : nsol=0
684 0 : do n=1,3
685 0 : xr=real(cx3(n,i))
686 0 : xi=aimag(cx3(n,i))
687 0 : if(abs(xi).gt.abs(xr)*eps) cycle
688 0 : if(xr.gt.r(i)) cycle
689 0 : if(xr.lt.rdry(i)*(1._r8-eps)) cycle
690 0 : if(xr.ne.xr) cycle
691 0 : r(i)=xr
692 0 : nsol=n
693 : end do
694 0 : if(nsol.eq.0)then
695 : write(iulog,*) &
696 0 : 'ccm kohlerc - no real(r8) solution found (cubic)'
697 0 : write(iulog,*)'roots =', (cx3(n,i),n=1,3)
698 0 : write(iulog,*)'p0-p2 =', p30(i), p31(i), p32(i)
699 0 : write(iulog,*)'rh=',s(i)
700 0 : write(iulog,*)'setting radius to dry radius=',rdry(i)
701 0 : r(i)=rdry(i)
702 : ! stop
703 : endif
704 : endif
705 0 : r3=r(i)
706 : ! now interpolate between quartic, cubic solutions
707 0 : r(i)=(r4*(1._r8-s(i))+r3*(s(i)-1._r8+eps))/eps
708 : endif
709 :
710 361584000 : 100 continue
711 :
712 : ! bound and convert from microns to m
713 723168000 : do i=1,im
714 361584000 : r(i) = min(r(i),30._r8) ! upper bound based on 1 day lifetime
715 723168000 : rwet_out(i) = r(i)*1.e-6_r8
716 : end do
717 :
718 361584000 : return
719 : end subroutine modal_aero_kohler
720 :
721 :
722 : !-----------------------------------------------------------------------
723 0 : subroutine makoh_cubic( cx, p2, p1, p0, im )
724 : !
725 : ! solves x**3 + p2 x**2 + p1 x + p0 = 0
726 : ! where p0, p1, p2 are real
727 : !
728 : integer, parameter :: imx=200
729 : integer :: im
730 : real(r8) :: p0(imx), p1(imx), p2(imx)
731 : complex(r8) :: cx(3,imx)
732 :
733 : integer :: i
734 : real(r8) :: eps, q(imx), r(imx), sqrt3, third
735 : complex(r8) :: ci, cq, crad(imx), cw, cwsq, cy(imx), cz(imx)
736 :
737 : save eps
738 : data eps/1.e-20_r8/
739 :
740 0 : third=1._r8/3._r8
741 0 : ci=cmplx(0._r8,1._r8,r8)
742 0 : sqrt3=sqrt(3._r8)
743 0 : cw=0.5_r8*(-1+ci*sqrt3)
744 0 : cwsq=0.5_r8*(-1-ci*sqrt3)
745 :
746 0 : do i=1,im
747 0 : if(p1(i).eq.0._r8)then
748 : ! completely insoluble particle
749 0 : cx(1,i)=(-p0(i))**third
750 0 : cx(2,i)=cx(1,i)
751 0 : cx(3,i)=cx(1,i)
752 : else
753 0 : q(i)=p1(i)/3._r8
754 0 : r(i)=p0(i)/2._r8
755 0 : crad(i)=r(i)*r(i)+q(i)*q(i)*q(i)
756 0 : crad(i)=sqrt(crad(i))
757 :
758 0 : cy(i)=r(i)-crad(i)
759 0 : if (abs(cy(i)).gt.eps) cy(i)=cy(i)**third
760 0 : cq=q(i)
761 0 : cz(i)=-cq/cy(i)
762 :
763 0 : cx(1,i)=-cy(i)-cz(i)
764 0 : cx(2,i)=-cw*cy(i)-cwsq*cz(i)
765 0 : cx(3,i)=-cwsq*cy(i)-cw*cz(i)
766 : endif
767 : enddo
768 :
769 0 : return
770 : end subroutine makoh_cubic
771 :
772 :
773 : !-----------------------------------------------------------------------
774 271140364 : subroutine makoh_quartic( cx, p3, p2, p1, p0, im )
775 :
776 : ! solves x**4 + p3 x**3 + p2 x**2 + p1 x + p0 = 0
777 : ! where p0, p1, p2, p3 are real
778 : !
779 : integer, parameter :: imx=200
780 : integer :: im
781 : real(r8) :: p0(imx), p1(imx), p2(imx), p3(imx)
782 : complex(r8) :: cx(4,imx)
783 :
784 : integer :: i
785 : real(r8) :: third, q(imx), r(imx)
786 : complex(r8) :: cb(imx), cb0(imx), cb1(imx), &
787 : crad(imx), cy(imx), czero
788 :
789 :
790 271140364 : czero=cmplx(0.0_r8,0.0_r8,r8)
791 271140364 : third=1._r8/3._r8
792 :
793 542280728 : do 10 i=1,im
794 :
795 271140364 : q(i)=-p2(i)*p2(i)/36._r8+(p3(i)*p1(i)-4*p0(i))/12._r8
796 : r(i)=-(p2(i)/6)**3+p2(i)*(p3(i)*p1(i)-4*p0(i))/48._r8 &
797 271140364 : +(4*p0(i)*p2(i)-p0(i)*p3(i)*p3(i)-p1(i)*p1(i))/16
798 :
799 271140364 : crad(i)=r(i)*r(i)+q(i)*q(i)*q(i)
800 271140364 : crad(i)=sqrt(crad(i))
801 :
802 271140364 : cb(i)=r(i)-crad(i)
803 271140364 : if(cb(i).eq.czero)then
804 : ! insoluble particle
805 0 : cx(1,i)=(-p1(i))**third
806 0 : cx(2,i)=cx(1,i)
807 0 : cx(3,i)=cx(1,i)
808 0 : cx(4,i)=cx(1,i)
809 : else
810 271140364 : cb(i)=cb(i)**third
811 :
812 271140364 : cy(i)=-cb(i)+q(i)/cb(i)+p2(i)/6
813 :
814 271140364 : cb0(i)=sqrt(cy(i)*cy(i)-p0(i))
815 271140364 : cb1(i)=(p3(i)*cy(i)-p1(i))/(2*cb0(i))
816 :
817 271140364 : cb(i)=p3(i)/2+cb1(i)
818 271140364 : crad(i)=cb(i)*cb(i)-4*(cy(i)+cb0(i))
819 271140364 : crad(i)=sqrt(crad(i))
820 271140364 : cx(1,i)=(-cb(i)+crad(i))/2._r8
821 271140364 : cx(2,i)=(-cb(i)-crad(i))/2._r8
822 :
823 271140364 : cb(i)=p3(i)/2-cb1(i)
824 271140364 : crad(i)=cb(i)*cb(i)-4*(cy(i)-cb0(i))
825 271140364 : crad(i)=sqrt(crad(i))
826 271140364 : cx(3,i)=(-cb(i)+crad(i))/2._r8
827 271140364 : cx(4,i)=(-cb(i)-crad(i))/2._r8
828 : endif
829 271140364 : 10 continue
830 :
831 271140364 : return
832 : end subroutine makoh_quartic
833 :
834 : !----------------------------------------------------------------------
835 0 : subroutine calc_h2so4_equilib_mixrat( temp, pres, qh2o, dmean, &
836 : qh2so4_equilib, wtpct, sulden )
837 :
838 : implicit none
839 :
840 : real(r8), intent(in) :: temp ! temperature (K)
841 : real(r8), intent(in) :: pres ! pressure (Pa)
842 : real(r8), intent(in) :: qh2o ! water vapor specific humidity (kg/kg)
843 : real(r8), intent(in) :: dmean ! mean diameter of particles in each mode
844 : real(r8), intent(out) :: qh2so4_equilib ! h2so4 saturation mixing ratios over the particles (mol/mol)
845 : real(r8), intent(out) :: wtpct ! sulfate composition, weight % H2SO4
846 : real(r8), intent(out) :: sulden ! sulfate aerosol mass density (g/cm3)
847 :
848 : ! Local declarations
849 : real(r8) :: qh2o_kelvin ! water vapor specific humidity adjusted for Kelvin effect (kg/kg)
850 : real(r8) :: wtpct_flat ! sulfate composition over a flat surface, weight % H2SO4
851 : real(r8) :: fk1, fk4, fk4_1, fk4_2
852 : real(r8) :: factor_kulm ! Kulmala correction terms
853 : real(r8) :: en, t, sig1, sig2, frac, surf_tens, surf_tens_mode, dsigma_dwt
854 : real(r8) :: den1, den2, sulfate_density, drho_dwt
855 : real(r8) :: akelvin, expon, akas, rkelvinH2O, rkelvinH2O_a, rkelvinH2O_b
856 : real(r8) :: sulfequil, r
857 : real(r8), parameter :: t0_kulm = 340._r8 ! T0 set in the low end of the Ayers measurement range (338-445K)
858 : real(r8), parameter :: t_crit_kulm = 905._r8 ! Critical temperature = 1.5 * Boiling point
859 : real(r8), parameter :: fk0 = -10156._r8 / t0_kulm + 16.259_r8 ! Log(Kulmala correction factor)
860 : real(r8), parameter :: fk2 = 1._r8 / t0_kulm
861 : real(r8), parameter :: fk3 = 0.38_r8 / (t_crit_kulm - t0_kulm)
862 : real(r8), parameter :: RGAS = 8.31430e+07_r8 ! ideal gas constant (erg/mol/K)
863 : real(r8), parameter :: wtmol_h2so4 = 98.078479_r8 ! molecular weight of sulphuric acid
864 : real(r8), parameter :: wtmol_h2o = 18.015280_r8 ! molecular weight of water vapor
865 : real(r8), parameter :: wtmol_air = 28.9644_r8 ! molecular weight of dry air
866 : real(r8) :: gv_wt_abcd_en(6,95), gvh2ovp(95)
867 : real(r8) :: dnwtp(46), dnc0(46), dnc1(46)
868 : real(r8) :: stwtp(15), stc0(15), stc1(15)
869 : integer :: i, k
870 :
871 :
872 : data stwtp/0._r8, 23.8141_r8, 38.0279_r8, 40.6856_r8, 45.335_r8, 52.9305_r8, &
873 : 56.2735_r8, 59.8557_r8, 66.2364_r8, 73.103_r8, 79.432_r8, 85.9195_r8, &
874 : 91.7444_r8, 97.6687_r8, 100._r8/
875 :
876 : data stc0/117.564_r8, 103.303_r8, 101.796_r8, 100.42_r8, 98.4993_r8, 91.8866_r8, &
877 : 88.3033_r8, 86.5546_r8, 84.471_r8, 81.2939_r8, 79.3556_r8, 75.608_r8, &
878 : 70.0777_r8, 63.7412_r8, 61.4591_r8 /
879 :
880 : data stc1/-0.153641_r8, -0.0982007_r8, -0.0872379_r8, -0.0818509_r8, &
881 : -0.0746702_r8, -0.0522399_r8, -0.0407773_r8, -0.0357946_r8, -0.0317062_r8, &
882 : -0.025825_r8, -0.0267212_r8, -0.0269204_r8, -0.0276187_r8, -0.0302094_r8, &
883 : -0.0303081_r8 /
884 :
885 :
886 : data dnwtp / 0._r8, 1._r8, 5._r8, 10._r8, 20._r8, 25._r8, 30._r8, 35._r8, 40._r8, &
887 : 41._r8, 45._r8, 50._r8, 53._r8, 55._r8, 56._r8, 60._r8, 65._r8, 66._r8, 70._r8, &
888 : 72._r8, 73._r8, 74._r8, 75._r8, 76._r8, 78._r8, 79._r8, 80._r8, 81._r8, 82._r8, &
889 : 83._r8, 84._r8, 85._r8, 86._r8, 87._r8, 88._r8, 89._r8, 90._r8, 91._r8, 92._r8, &
890 : 93._r8, 94._r8, 95._r8, 96._r8, 97._r8, 98._r8, 100._r8 /
891 :
892 : data dnc0 / 1._r8, 1.13185_r8, 1.17171_r8, 1.22164_r8, 1.3219_r8, 1.37209_r8, &
893 : 1.42185_r8, 1.4705_r8, 1.51767_r8, 1.52731_r8, 1.56584_r8, 1.61834_r8, 1.65191_r8, &
894 : 1.6752_r8, 1.68708_r8, 1.7356_r8, 1.7997_r8, 1.81271_r8, 1.86696_r8, 1.89491_r8, &
895 : 1.9092_r8, 1.92395_r8, 1.93904_r8, 1.95438_r8, 1.98574_r8, 2.00151_r8, 2.01703_r8, &
896 : 2.03234_r8, 2.04716_r8, 2.06082_r8, 2.07363_r8, 2.08461_r8, 2.09386_r8, 2.10143_r8,&
897 : 2.10764_r8, 2.11283_r8, 2.11671_r8, 2.11938_r8, 2.12125_r8, 2.1219_r8, 2.12723_r8, &
898 : 2.12654_r8, 2.12621_r8, 2.12561_r8, 2.12494_r8, 2.12093_r8 /
899 :
900 : data dnc1 / 0._r8, -0.000435022_r8, -0.000479481_r8, -0.000531558_r8, -0.000622448_r8, &
901 : -0.000660866_r8, -0.000693492_r8, -0.000718251_r8, -0.000732869_r8, -0.000735755_r8, &
902 : -0.000744294_r8, -0.000761493_r8, -0.000774238_r8, -0.00078392_r8, -0.000788939_r8, &
903 : -0.00080946_r8, -0.000839848_r8, -0.000845825_r8, -0.000874337_r8, -0.000890074_r8, &
904 : -0.00089873_r8, -0.000908778_r8, -0.000920012_r8, -0.000932184_r8, -0.000959514_r8, &
905 : -0.000974043_r8, -0.000988264_r8, -0.00100258_r8, -0.00101634_r8, -0.00102762_r8, &
906 : -0.00103757_r8, -0.00104337_r8, -0.00104563_r8, -0.00104458_r8, -0.00104144_r8, &
907 : -0.00103719_r8, -0.00103089_r8, -0.00102262_r8, -0.00101355_r8, -0.00100249_r8, &
908 : -0.00100934_r8, -0.000998299_r8, -0.000990961_r8, -0.000985845_r8, -0.000984529_r8, &
909 : -0.000989315_r8 /
910 :
911 : ! Saturation vapor pressure of sulfuric acid
912 : !
913 : ! Limit extrapolation at extreme temperatures
914 0 : t=min(max(temp,140._r8),450._r8)
915 :
916 : !! Calculate the weight % H2SO4 composition of sulfate
917 0 : call calc_h2so4_wtpct(t, pres, qh2o, wtpct_flat)
918 :
919 : !! Calculate surface tension (erg/cm2) of sulfate of
920 : !! different compositions as a linear function of temperature.
921 0 : i = 2 ! minimum wt% is 29.517
922 0 : do while (wtpct_flat.gt.stwtp(i))
923 0 : i = i + 1
924 : end do
925 0 : sig1 = stc0(i-1) + stc1(i-1) * t
926 0 : sig2 = stc0(i) + stc1(i) * t
927 : ! calculate derivative needed later for kelvin factor for h2o
928 0 : dsigma_dwt = (sig2-sig1) / (stwtp(i)-stwtp(i-1))
929 0 : surf_tens = sig1 + dsigma_dwt*(wtpct_flat-stwtp(i))
930 :
931 : !! Calculate density (g/cm3) of sulfate of
932 : !! different compositions as a linear function of temperature.
933 0 : i = 6 ! minimum wt% is 29.517
934 0 : do while (wtpct_flat .gt. dnwtp(i))
935 0 : i = i + 1
936 : end do
937 0 : den1 = dnc0(i-1) + dnc1(i-1) * t
938 0 : den2 = dnc0(i) + dnc1(i) * t
939 : ! Calculate derivative needed later for Kelvin factor for H2O
940 0 : drho_dwt = (den2-den1) / (dnwtp(i)-dnwtp(i-1))
941 0 : sulfate_density = den1 + drho_dwt * (wtpct_flat-dnwtp(i-1))
942 :
943 0 : r=dmean*100._r8/2._r8 ! calcuate mode radius (cm) from diameter (m)
944 :
945 : ! Adjust for Kelvin effect for water
946 : rkelvinH2O_b = 1 + wtpct_flat * drho_dwt / &
947 0 : sulfate_density - 3._r8 * wtpct_flat * dsigma_dwt / (2._r8*surf_tens)
948 :
949 : rkelvinH2O_a = 2._r8 * wtmol_h2so4 * surf_tens / &
950 0 : (sulfate_density * RGAS * t * r)
951 :
952 0 : rkelvinH2O = exp (rkelvinH2O_a*rkelvinH2O_b)
953 :
954 0 : qh2o_kelvin = qh2o/rkelvinH2O
955 0 : call calc_h2so4_wtpct(t, pres, qh2o_kelvin, wtpct)
956 :
957 :
958 0 : wtpct=max(wtpct,wtpct_flat)
959 :
960 : ! Parameterized fit to Giauque's (1959) enthalpies v. wt %:
961 0 : en = 4.184_r8 * (23624.8_r8 - 1.14208e8_r8 / ((wtpct - 105.318_r8)**2 + 4798.69_r8))
962 0 : en = max(en, 0.0_r8)
963 :
964 : !! Calculate surface tension (erg/cm2) of sulfate of
965 : !! different compositions as a linear function of temperature.
966 0 : i = 2 ! minimum wt% is 29.517
967 0 : do while (wtpct.gt.stwtp(i))
968 0 : i=i+1
969 : end do
970 0 : sig1=stc0(i-1)+stc1(i-1)*t
971 0 : sig2=stc0(i)+stc1(i)*t
972 0 : frac=(stwtp(i)-wtpct)/(stwtp(i)-stwtp(i-1))
973 0 : surf_tens_mode=sig1*frac+sig2*(1.0_r8-frac)
974 :
975 : !! Calculate density (g/cm3) of sulfate of
976 : !! different compositions as a linear function of temperature.
977 0 : i = 6 ! minimum wt% is 29.517
978 0 : do while (wtpct .gt. dnwtp(i))
979 0 : i=i+1
980 : end do
981 0 : den1=dnc0(i-1)+dnc1(i-1)*t
982 0 : den2=dnc0(i)+dnc1(i)*t
983 0 : frac=(dnwtp(i)-wtpct)/(dnwtp(i)-dnwtp(i-1))
984 0 : sulden=den1*frac+den2*(1.0_r8-frac)
985 :
986 : ! Ayers' (1980) fit to sulfuric acid equilibrium vapor pressure:
987 : ! (Remember this is the log)
988 : ! SULFEQ=-10156/Temp+16.259-En/(8.3143*Temp)
989 : !
990 : ! Kulmala correction (J. CHEM. PHYS. V.93, No.1, 1 July 1990)
991 0 : fk1 = -1._r8 / t
992 0 : fk4_1 = log(t0_kulm / t)
993 0 : fk4_2 = t0_kulm / t
994 0 : fk4 = 1.0_r8 + fk4_1 - fk4_2
995 0 : factor_kulm = fk1 + fk2 + fk3 * fk4
996 :
997 : ! This is for pure H2SO4
998 0 : sulfequil = fk0 + 10156._r8 * factor_kulm
999 :
1000 : ! Adjust for WTPCT composition:
1001 0 : sulfequil = sulfequil - en / (8.3143_r8 * t)
1002 :
1003 : ! Take the exponential:
1004 0 : sulfequil = exp(sulfequil)
1005 :
1006 : ! Convert atmospheres ==> Pa
1007 0 : sulfequil = sulfequil * 1.01325e5_r8
1008 :
1009 : ! Convert Pa ==> mol/mol
1010 0 : sulfequil = sulfequil / pres
1011 :
1012 : ! Calculate Kelvin curvature factor for H2SO4 interactively with temperature:
1013 : ! (g/mol)*(erg/cm2)/(K * g/cm3 * erg/mol/K) = cm
1014 0 : akelvin = 2._r8 * wtmol_h2so4 * surf_tens_mode / (t * sulden * RGAS)
1015 :
1016 0 : expon = akelvin / r ! divide by mode radius (cm)
1017 0 : expon = max(-100._r8, expon)
1018 0 : expon = min(100._r8, expon)
1019 0 : akas = exp( expon )
1020 0 : qh2so4_equilib = sulfequil * akas ! reduce H2SO4 equilibrium mixing ratio by Kelvin curvature factor
1021 :
1022 0 : return
1023 : end subroutine calc_h2so4_equilib_mixrat
1024 :
1025 :
1026 : !----------------------------------------------------------------------
1027 0 : subroutine calc_h2so4_wtpct( temp, pres, qh2o, wtpct )
1028 :
1029 : !! This function calculates the weight % H2SO4 composition of
1030 : !! sulfate aerosol, using Tabazadeh et. al. (GRL, 1931, 1997).
1031 : !! Rated for T=185-260K, activity=0.01-1.0
1032 : !!
1033 : !! Argument list input:
1034 : !! temp = temperature (K)
1035 : !! pres = atmospheric pressure (Pa)
1036 : !! qh2o = water specific humidity (kg/kg)
1037 : !!
1038 : !! Output:
1039 : !! wtpct = weight % H2SO4 in H2O/H2SO4 particle (0-100)
1040 : !!
1041 : !! @author Mike Mills
1042 : !! @ version October 2013
1043 :
1044 : use wv_saturation, only: qsat_water
1045 :
1046 : implicit none
1047 :
1048 : real(r8), intent(in) :: temp ! temperature (K)
1049 : real(r8), intent(in) :: pres ! pressure (Pa)
1050 : real(r8), intent(in) :: qh2o ! water vapor specific humidity (kg/kg)
1051 : real(r8), intent(out) :: wtpct ! sulfate weight % H2SO4 composition
1052 :
1053 : ! Local declarations
1054 : real(r8) :: atab1,btab1,ctab1,dtab1,atab2,btab2,ctab2,dtab2
1055 : real(r8) :: contl, conth, contt, conwtp
1056 : real(r8) :: activ
1057 : real(r8) :: es ! saturation vapor pressure over water (Pa) (dummy)
1058 : real(r8) :: qs ! saturation specific humidity over water (kg/kg)
1059 :
1060 : ! calculate saturation specific humidity over pure water, qs (kg/kg)
1061 0 : call qsat_water(temp, pres, es, qs)
1062 :
1063 : ! Activity = water specific humidity (kg/kg) / equilibrium water (kg/kg)
1064 0 : activ = qh2o/qs
1065 :
1066 0 : if (activ.lt.0.05_r8) then
1067 0 : activ = max(activ,1.e-6_r8) ! restrict minimum activity
1068 0 : atab1 = 12.37208932_r8
1069 0 : btab1 = -0.16125516114_r8
1070 0 : ctab1 = -30.490657554_r8
1071 0 : dtab1 = -2.1133114241_r8
1072 0 : atab2 = 13.455394705_r8
1073 0 : btab2 = -0.1921312255_r8
1074 0 : ctab2 = -34.285174607_r8
1075 0 : dtab2 = -1.7620073078_r8
1076 0 : elseif (activ.ge.0.05_r8.and.activ.le.0.85_r8) then
1077 : atab1 = 11.820654354_r8
1078 : btab1 = -0.20786404244_r8
1079 : ctab1 = -4.807306373_r8
1080 : dtab1 = -5.1727540348_r8
1081 : atab2 = 12.891938068_r8
1082 : btab2 = -0.23233847708_r8
1083 : ctab2 = -6.4261237757_r8
1084 : dtab2 = -4.9005471319_r8
1085 0 : elseif (activ.gt.0.85_r8) then
1086 0 : activ = min(activ,1._r8) ! restrict maximum activity
1087 0 : atab1 = -180.06541028_r8
1088 0 : btab1 = -0.38601102592_r8
1089 0 : ctab1 = -93.317846778_r8
1090 0 : dtab1 = 273.88132245_r8
1091 0 : atab2 = -176.95814097_r8
1092 0 : btab2 = -0.36257048154_r8
1093 0 : ctab2 = -90.469744201_r8
1094 0 : dtab2 = 267.45509988_r8
1095 : else
1096 0 : write(iulog,*) 'calc_h2so4_wtpct: invalid activity: activ,qh2o,qs,temp,pres=',activ,qh2o,qs,temp,pres
1097 0 : call endrun( 'calc_h2so4_wtpct error' )
1098 0 : return
1099 : endif
1100 :
1101 0 : contl = atab1*(activ**btab1)+ctab1*activ+dtab1
1102 0 : conth = atab2*(activ**btab2)+ctab2*activ+dtab2
1103 :
1104 0 : contt = contl + (conth-contl) * ((temp -190._r8)/70._r8)
1105 0 : conwtp = (contt*98._r8) + 1000._r8
1106 :
1107 0 : wtpct = (100._r8*contt*98._r8)/conwtp
1108 0 : wtpct = min(max(wtpct,25._r8),100._r8) ! restrict between 1 and 100 %
1109 :
1110 0 : return
1111 : end subroutine calc_h2so4_wtpct
1112 :
1113 :
1114 : !----------------------------------------------------------------------
1115 :
1116 : end module modal_aero_wateruptake
1117 :
1118 :
|
