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