Line data Source code
1 : module hetfrz_classnuc
2 :
3 : !-----------------------------------------------------------------------
4 : !
5 : ! Purpose: Calculate heterogeneous freezing rates from classical nucleation theory
6 : !
7 : ! Public interfaces:
8 : !
9 : ! hetfrz_classnuc_init
10 : ! hetfrz_classnuc_calc
11 : !
12 : ! Author:
13 : ! Corinna Hoose, UiO, May 2009
14 : ! Yong Wang and Xiaohong Liu, UWyo, 12/2012,
15 : ! implement in CAM5 and constrain uncertain parameters using natural dust and
16 : ! BC(soot) datasets.
17 : ! Yong Wang and Xiaohong Liu, UWyo, 05/2013, implement the PDF-contact angle approach:
18 : ! Y. Wang et al., Atmos. Chem. Phys., 2014. https://doi.org/10.5194/acp-14-10411-2014
19 : ! Jack Chen, NCAR, 09/2015, modify calculation of dust activation fraction.
20 : !
21 : !-----------------------------------------------------------------------
22 :
23 : use shr_kind_mod, only: r8 => shr_kind_r8
24 : use wv_saturation, only: svp_water
25 : use shr_spfn_mod, only: erf => shr_spfn_erf
26 :
27 : use physconst, only: pi, planck, boltz, mwso4, amu, pstd
28 :
29 : implicit none
30 : private
31 :
32 : public :: hetfrz_classnuc_init, hetfrz_classnuc_calc
33 :
34 : real(r8) :: rair
35 : real(r8) :: cpair
36 : real(r8) :: rh2o
37 : real(r8) :: rhoh2o
38 : real(r8) :: mwh2o
39 : real(r8) :: tmelt
40 :
41 : !*****************************************************************************
42 : ! PDF theta model
43 : !*****************************************************************************
44 : ! some variables for PDF theta model
45 : ! immersion freezing
46 : !
47 : ! With the original value of pdf_n_theta set to 101 the dust activation
48 : ! fraction between -15 and 0 C could be overestimated. This problem was
49 : ! eliminated by increasing pdf_n_theta to 301. To reduce the expense of
50 : ! computing the dust activation fraction the integral is only evaluated
51 : ! where dim_theta is non-zero. This was determined to be between
52 : ! dim_theta index values of 53 through 113. These loop bounds are
53 : ! hardcoded in the variables i1 and i2.
54 :
55 : logical, parameter :: pdf_imm_in = .true.
56 : integer, parameter :: pdf_n_theta = 301
57 : integer, parameter :: i1 = 53
58 : integer, parameter :: i2 = 113
59 :
60 : real(r8) :: dim_theta(pdf_n_theta) = -huge(1._r8)
61 : real(r8) :: pdf_imm_theta(pdf_n_theta) = 0.0_r8
62 : real(r8) :: pdf_d_theta = -huge(1._r8)
63 : real(r8) :: dim_f_imm(pdf_n_theta) = 0.0_r8
64 :
65 : integer :: iulog
66 :
67 : real(r8), parameter :: n1 = 1.e19_r8 ! number of water molecules in contact with unit area of substrate [m-2]
68 : real(r8), parameter :: rhplanck = 1._r8/planck
69 : real(r8), parameter :: nus = 1.e13_r8 ! frequ. of vibration [s-1] higher freq. (as in P&K, consistent with Anupam's data)
70 : real(r8), parameter :: rhwincloud = 0.98_r8 ! 98% RH in mixed-phase clouds (Korolev & Isaac, JAS 2006)
71 :
72 : logical, parameter :: tot_in = .false.
73 :
74 : real(r8) :: bc_limfac = -huge(1._r8) ! soot ice nucleating fraction
75 : real(r8) :: dust_limfac = -huge(1._r8) ! dust ice nucleating fraction
76 :
77 : !===================================================================================================
78 : contains
79 : !===================================================================================================
80 :
81 1536 : subroutine hetfrz_classnuc_init( &
82 : rair_in, cpair_in, rh2o_in, rhoh2o_in, mwh2o_in, &
83 : tmelt_in, iulog_in, bc_limfac_in, dust_limfac_in)
84 :
85 : real(r8), intent(in) :: rair_in
86 : real(r8), intent(in) :: cpair_in
87 : real(r8), intent(in) :: rh2o_in
88 : real(r8), intent(in) :: rhoh2o_in
89 : real(r8), intent(in) :: mwh2o_in
90 : real(r8), intent(in) :: tmelt_in
91 : integer, intent(in) :: iulog_in
92 : real(r8), intent(in) :: bc_limfac_in
93 : real(r8), intent(in) :: dust_limfac_in
94 :
95 1536 : rair = rair_in
96 1536 : cpair = cpair_in
97 1536 : rh2o = rh2o_in
98 1536 : rhoh2o = rhoh2o_in
99 1536 : mwh2o = mwh2o_in
100 1536 : tmelt = tmelt_in
101 1536 : iulog = iulog_in
102 :
103 1536 : bc_limfac = bc_limfac_in
104 1536 : dust_limfac = dust_limfac_in
105 :
106 : ! Initialize all the PDF theta variables:
107 1536 : if (pdf_imm_in) then
108 : call hetfrz_classnuc_init_pdftheta()
109 : end if
110 :
111 1536 : end subroutine hetfrz_classnuc_init
112 :
113 : !===================================================================================================
114 :
115 1270492981 : subroutine hetfrz_classnuc_calc(ntypes, types,&
116 : deltat, T, p, supersatice, &
117 1270492981 : fn, &
118 : r3lx, icnlx, &
119 : frzbcimm, frzduimm, &
120 : frzbccnt, frzducnt, &
121 : frzbcdep, frzdudep, &
122 1270492981 : hetraer, wact_factor, dstcoat, &
123 1270492981 : total_aer_num, uncoated_aer_num, &
124 1270492981 : total_interstitial_aer_num, total_cloudborne_aer_num, errstring)
125 :
126 : integer, intent(in) :: ntypes
127 : character(len=*), intent(in) :: types(ntypes)
128 : real(r8), intent(in) :: deltat ! timestep [s]
129 : real(r8), intent(in) :: T ! temperature [K]
130 : real(r8), intent(in) :: p ! pressure [Pa]
131 : real(r8), intent(in) :: supersatice ! supersaturation ratio wrt ice at 100%rh over water [ ]
132 : real(r8), intent(in) :: r3lx ! volume mean drop radius [m]
133 : real(r8), intent(in) :: icnlx ! in-cloud droplet concentration [cm-3]
134 :
135 : real(r8), intent(in) :: fn(ntypes) ! fraction activated [ ] for cloud borne aerosol number
136 : ! index values are 1:bc, 2:dust_a1, 3:dust_a3
137 : real(r8), intent(in) :: hetraer(ntypes) ! bc and dust mass mean radius [m]
138 : real(r8), intent(in) :: wact_factor(ntypes) ! water activity factor -- density*(1.-(OC+BC)/(OC+BC+SO4)) [mug m-3]
139 : real(r8), intent(in) :: dstcoat(ntypes) ! coated fraction
140 : real(r8), intent(in) :: total_aer_num(ntypes) ! total bc and dust number concentration(interstitial+cloudborne) [#/cm^3]
141 : real(r8), intent(in) :: uncoated_aer_num(ntypes) ! uncoated bc and dust number concentration(interstitial)
142 : real(r8), intent(in) :: total_interstitial_aer_num(ntypes) ! total bc and dust concentration(interstitial)
143 : real(r8), intent(in) :: total_cloudborne_aer_num(ntypes) ! total bc and dust concentration(cloudborne)
144 :
145 : real(r8), target, intent(out) :: frzbcimm ! het. frz by BC immersion nucleation [cm-3 s-1]
146 : real(r8), target, intent(out) :: frzduimm ! het. frz by dust immersion nucleation [cm-3 s-1]
147 : real(r8), target, intent(out) :: frzbccnt ! het. frz by BC contact nucleation [cm-3 s-1]
148 : real(r8), target, intent(out) :: frzducnt ! het. frz by dust contact nucleation [cm-3 s-1]
149 : real(r8), target, intent(out) :: frzbcdep ! het. frz by BC deposition nucleation [cm-3 s-1]
150 : real(r8), target, intent(out) :: frzdudep ! het. frz by dust deposition nucleation [cm-3 s-1]
151 :
152 : character(len=*), intent(out) :: errstring
153 :
154 : ! local variables
155 :
156 : real(r8) :: tc
157 : real(r8) :: vwice
158 : real(r8) :: rhoice
159 : real(r8) :: sigma_iw ! [J/m2]
160 : real(r8) :: sigma_iv ! [J/m2]
161 : real(r8) :: eswtr ! [Pa]
162 :
163 : real(r8) :: rgimm ! critical germ size
164 : real(r8) :: rgdep
165 : real(r8) :: dg0dep ! homogeneous energy of germ formation
166 : real(r8) :: dg0cnt
167 : real(r8) :: Adep ! prefactors
168 : real(r8) :: Acnt
169 :
170 : !********************************************************
171 : ! Hoose et al., 2010 fitting parameters
172 : !********************************************************
173 : !freezing parameters for immersion freezing
174 : !real(r8),parameter :: theta_imm_bc = 40.17 ! contact angle [deg], converted to rad later
175 : !real(r8),parameter :: dga_imm_bc = 14.4E-20 ! activation energy [J]
176 : !real(r8),parameter :: theta_imm_dust = 30.98 ! contact angle [deg], converted to rad later
177 : !real(r8),parameter :: dga_imm_dust = 15.7E-20 ! activation energy [J]
178 : !freezing parameters for deposition nucleation
179 : !real(r8),parameter :: theta_dep_dust = 12.7 ! contact angle [deg], converted to rad later !Zimmermann et al (2008), illite
180 : !real(r8),parameter :: dga_dep_dust = -6.21E-21 ! activation energy [J]
181 : !real(r8),parameter :: theta_dep_bc = 28. ! contact angle [deg], converted to rad later !Moehler et al (2005), soot
182 : !real(r8),parameter :: dga_dep_bc = -2.E-19 ! activation energy [J]
183 : !********************************************************
184 : ! Wang et al., 2014 fitting parameters
185 : !********************************************************
186 : ! freezing parameters for immersion freezing
187 : real(r8),parameter :: theta_imm_bc = 48.0_r8 ! contact angle [deg], converted to rad later !DeMott et al (1990)
188 : real(r8),parameter :: dga_imm_bc = 14.15E-20_r8 ! activation energy [J]
189 : real(r8),parameter :: theta_imm_dust = 46.0_r8 ! contact angle [deg], converted to rad later !DeMott et al (2011) SD
190 : real(r8),parameter :: dga_imm_dust = 14.75E-20_r8 ! activation energy [J]
191 : ! freezing parameters for deposition nucleation
192 : real(r8),parameter :: theta_dep_dust = 20.0_r8 ! contact angle [deg], converted to rad later !Koehler et al (2010) SD
193 : real(r8),parameter :: dga_dep_dust = -8.1E-21_r8 ! activation energy [J]
194 : real(r8),parameter :: theta_dep_bc = 28._r8 ! contact angle [deg], converted to rad later !Moehler et al (2005), soot
195 : real(r8),parameter :: dga_dep_bc = -2.E-19_r8 ! activation energy [J]
196 :
197 : ! form factor
198 : ! only consider flat surfaces due to uncertainty of curved surfaces
199 : real(r8),parameter :: m_depcnt_bc = COS(theta_dep_bc*pi/180._r8)
200 : real(r8),parameter :: f_depcnt_bc = (2+m_depcnt_bc)*(1-m_depcnt_bc)**2/4._r8
201 :
202 : real(r8),parameter :: m_depcnt_dst = COS(theta_dep_dust*pi/180._r8)
203 : real(r8),parameter :: f_depcnt_dust = (2+m_depcnt_dst)*(1-m_depcnt_dst)**2/4._r8
204 :
205 : real(r8),parameter :: m_imm_bc = COS(theta_imm_bc*pi/180._r8)
206 : real(r8),parameter :: f_imm_bc = (2+m_imm_bc)*(1-m_imm_bc)**2/4._r8
207 :
208 : real(r8),parameter :: m_imm_dust = COS(theta_imm_dust*pi/180._r8)
209 : real(r8),parameter :: f_imm_dust = (2+m_imm_dust)*(1-m_imm_dust)**2/4._r8
210 :
211 : real(r8) :: f_dep, f_cnt, f_imm
212 : real(r8) :: dga_dep, dga_imm
213 : real(r8) :: limfac
214 : real(r8) :: frzimm, frzcnt, frzdep
215 : real(r8), pointer :: frzimm_ptr, frzcnt_ptr, frzdep_ptr
216 :
217 : logical :: pdf_imm
218 :
219 : integer :: ispc
220 :
221 1270492981 : real(r8) :: ktherm(ntypes), kcoll(ntypes)
222 :
223 : real(r8), parameter :: Ktherm_bc = 4.2_r8 ! black carbon thermal conductivity [J/(m s K)]
224 : real(r8), parameter :: Ktherm_dst = 0.72_r8 ! clay thermal conductivity [J/(m s K)]
225 :
226 : !------------------------------------------------------------------------------------------------
227 :
228 1270492981 : errstring = ' '
229 :
230 1270492981 : nullify(frzimm_ptr)
231 1270492981 : nullify(frzcnt_ptr)
232 1270492981 : nullify(frzdep_ptr)
233 :
234 1270492981 : frzbcimm = 0._r8
235 1270492981 : frzbccnt= 0._r8
236 1270492981 : frzbcdep = 0._r8
237 1270492981 : frzduimm = 0._r8
238 1270492981 : frzducnt= 0._r8
239 1270492981 : frzdudep = 0._r8
240 :
241 : ! get saturation vapor pressure
242 1270492981 : eswtr = svp_water(t) ! 0 for liquid
243 :
244 1270492981 : tc = T - tmelt
245 1270492981 : rhoice = 916.7_r8-0.175_r8*tc-5.e-4_r8*tc**2
246 1270492981 : vwice = mwh2o*amu/rhoice
247 1270492981 : sigma_iw = (28.5_r8+0.25_r8*tc)*1E-3_r8
248 1270492981 : sigma_iv = (76.1_r8-0.155_r8*tc + 28.5_r8+0.25_r8*tc)*1E-3_r8
249 :
250 : ! critical germ size
251 1270492981 : rgimm = 2*vwice*sigma_iw/(boltz*T*LOG(supersatice))
252 :
253 : ! critical germ size
254 : ! assume 98% RH in mixed-phase clouds (Korolev & Isaac, JAS 2006)
255 1270492981 : rgdep=2*vwice*sigma_iv/(boltz*T*LOG(rhwincloud*supersatice))
256 :
257 : ! homogeneous energy of germ formation
258 1270492981 : dg0dep = 4*pi/3._r8*sigma_iv*rgdep**2
259 :
260 : ! prefactor
261 : ! attention: division of small numbers
262 1270492981 : Adep = (rhwincloud*eswtr)**2*(vwice/(mwh2o*amu))/(boltz*T*nus)*SQRT(sigma_iv/(boltz*T))
263 :
264 : ! homogeneous energy of germ formation
265 1270492981 : dg0cnt = 4*pi/3._r8*sigma_iv*rgimm**2
266 :
267 : ! prefactor
268 : ! attention: division of small numbers
269 1270492981 : Acnt = rhwincloud*eswtr*4*pi/(nus*SQRT(2*pi*mwh2o*amu*boltz*T))
270 :
271 6352464905 : do ispc = 1, ntypes
272 :
273 11434436829 : select case (trim(types(ispc)))
274 : case ('black-c')
275 2540985962 : ktherm(ispc) = ktherm_bc
276 : case ('dust')
277 2540985962 : ktherm(ispc) = ktherm_dst
278 : case default
279 0 : errstring = 'hetfrz_classnuc_calc ERROR: unrecognized aerosol type: '//trim(types(ispc))
280 10163943848 : return
281 : end select
282 : end do
283 :
284 1270492981 : call collkernel(T, p, eswtr, rhwincloud, r3lx, hetraer, Ktherm, Kcoll)
285 :
286 6352464905 : do ispc = 1, ntypes
287 :
288 10163943848 : select case (trim(types(ispc)))
289 : case ('black-c')
290 2540985962 : f_dep = f_depcnt_bc
291 2540985962 : f_cnt = f_depcnt_bc
292 2540985962 : f_imm = f_imm_bc
293 2540985962 : dga_dep = dga_dep_bc
294 2540985962 : dga_imm = dga_imm_bc
295 2540985962 : pdf_imm = .false.
296 2540985962 : limfac = bc_limfac
297 2540985962 : frzimm_ptr => frzbcimm
298 2540985962 : frzcnt_ptr => frzbccnt
299 2540985962 : frzdep_ptr => frzbcdep
300 : case ('dust')
301 2540985962 : f_dep = f_depcnt_dust
302 2540985962 : f_cnt = f_depcnt_dust
303 2540985962 : f_imm = f_imm_dust
304 2540985962 : dga_dep = dga_dep_dust
305 2540985962 : dga_imm = dga_imm_dust
306 2540985962 : pdf_imm = .true.
307 2540985962 : limfac = dust_limfac
308 2540985962 : frzimm_ptr => frzduimm
309 2540985962 : frzcnt_ptr => frzducnt
310 2540985962 : frzdep_ptr => frzdudep
311 : case default
312 0 : errstring = 'hetfrz_classnuc_calc ERROR: unrecognized aerosol type: '//trim(types(ispc))
313 10163943848 : return
314 : end select
315 :
316 : call hetfrz_classnuc_calc_rates( f_dep, f_cnt, f_imm, dga_dep, dga_imm, pdf_imm, limfac, &
317 : kcoll(ispc), hetraer(ispc), icnlx, r3lx, T, supersatice, sigma_iw, &
318 : rgimm, rgdep, dg0dep, Adep, dg0cnt, Acnt, vwice, deltat, &
319 : fn(ispc), wact_factor(ispc), dstcoat(ispc), &
320 : total_aer_num(ispc), total_interstitial_aer_num(ispc), total_cloudborne_aer_num(ispc), uncoated_aer_num(ispc), &
321 5081971924 : frzimm, frzcnt, frzdep, errstring )
322 :
323 : ! accumulate dust and bc frz rates
324 5081971924 : frzimm_ptr = frzimm_ptr + frzimm
325 5081971924 : frzcnt_ptr = frzcnt_ptr + frzcnt
326 6352464905 : frzdep_ptr = frzdep_ptr + frzdep
327 :
328 : end do
329 :
330 1270492981 : end subroutine hetfrz_classnuc_calc
331 :
332 5081971924 : subroutine hetfrz_classnuc_calc_rates( f_dep, f_cnt, f_imm, dga_dep, dga_imm, pdf_imm, limfac, &
333 : kcoll, mradius, icnlx, r3lx, T, supersatice, sigma_iw, &
334 : rgimm, rgdep, dg0dep, Adep, dg0cnt, Acnt, vwice, deltat, &
335 : fn, wact_factor, dstcoat, &
336 : total_aer_num, total_interstitial_aer_num, total_cloudborne_aer_num, uncoated_aer_num, &
337 : frzimm, frzcnt, frzdep, errstring )
338 :
339 : ! input
340 : real(r8), intent(in) :: f_dep ! deposition form factor
341 : real(r8), intent(in) :: f_cnt ! contact form factor
342 : real(r8), intent(in) :: f_imm ! immersion form factor
343 : real(r8), intent(in) :: dga_dep ! deposition activation energy [J]
344 : real(r8), intent(in) :: dga_imm ! immersion activation energy [J]
345 : logical, intent(in) :: pdf_imm ! PDF theta model switch (TRUE for dust)
346 : real(r8), intent(in) :: limfac ! Limit to 1% of available potential IN (for BC), no limit for dust
347 : real(r8), intent(in) :: kcoll ! collision kernel [cm3 s-1]
348 : real(r8), intent(in) :: mradius ! mass mean radius [m]
349 : real(r8), intent(in) :: icnlx ! in-cloud droplet concentration [cm-3]
350 : real(r8), intent(in) :: r3lx ! volume mean drop radius [m]
351 : real(r8), intent(in) :: T ! temperature [K]
352 : real(r8), intent(in) :: supersatice ! supersaturation ratio wrt ice at 100%rh over water [ ]
353 : real(r8), intent(in) :: sigma_iw ! [J/m2]
354 : real(r8), intent(in) :: rgimm ! critical germ size
355 : real(r8), intent(in) :: rgdep ! critical germ size
356 : real(r8), intent(in) :: dg0dep ! homogeneous energy of germ formation
357 : real(r8), intent(in) :: Adep ! deposition nucleation prefactor
358 : real(r8), intent(in) :: dg0cnt ! homogeneous energy of germ formation
359 : real(r8), intent(in) :: Acnt ! contact nucleation prefactor
360 :
361 : real(r8), intent(in) :: vwice
362 : real(r8), intent(in) :: deltat ! timestep [s]
363 : real(r8), intent(in) :: fn ! fraction activated [ ] for cloud borne aerosol number
364 : real(r8), intent(in) :: wact_factor ! water activity factor -- density*(1.-(OC+BC)/(OC+BC+SO4)) [mug m-3]
365 : real(r8), intent(in) :: dstcoat ! coated fraction
366 : real(r8), intent(in) :: total_aer_num ! total bc and dust number concentration(interstitial+cloudborne) [#/cm^3]
367 : real(r8), intent(in) :: total_interstitial_aer_num ! total bc and dust concentration(interstitial)
368 : real(r8), intent(in) :: total_cloudborne_aer_num ! total bc and dust concentration(cloudborne)
369 : real(r8), intent(in) :: uncoated_aer_num ! uncoated bc and dust number concentration(interstitial)
370 : ! output
371 : real(r8), intent(out) :: frzimm ! het. frz by immersion nucleation [cm-3 s-1]
372 : real(r8), intent(out) :: frzcnt ! het. frz by contact nucleation [cm-3 s-1]
373 : real(r8), intent(out) :: frzdep ! het. frz by deposition nucleation [cm-3 s-1]
374 :
375 : character(len=*), intent(out) :: errstring
376 :
377 : ! local vars
378 : real(r8) :: aw ! water activity [ ]
379 : real(r8) :: molal ! molality [moles/kg]
380 :
381 : real(r8) :: Aimm
382 : real(r8) :: Jdep
383 : real(r8) :: Jimm
384 : real(r8) :: Jcnt
385 : real(r8) :: dg0imm
386 : real(r8) :: rgimm_aer
387 : real(r8) :: sum_imm
388 : real(r8) :: dim_Jimm(pdf_n_theta)
389 :
390 : logical :: do_frz
391 : logical :: do_imm
392 :
393 : integer :: i
394 :
395 : !*****************************************************************************
396 : ! take water activity into account
397 : !*****************************************************************************
398 : ! solute effect
399 5081971924 : aw = 1._r8
400 5081971924 : molal = 0._r8
401 :
402 : ! The heterogeneous ice freezing temperatures of all IN generally decrease with
403 : ! increasing total solute mole fraction. Therefore, the large solution concentration
404 : ! will cause the freezing point depression and the ice freezing temperatures of all
405 : ! IN will get close to the homogeneous ice freezing temperatures. Since we take into
406 : ! account water activity for three heterogeneous freezing modes(immersion, deposition,
407 : ! and contact), we utilize interstitial aerosols(not cloudborne aerosols) to calculate
408 : ! water activity.
409 : ! If the index of IN is 0, it means three freezing modes of this aerosol are depressed.
410 :
411 : !calculate molality
412 5081971924 : if ( total_interstitial_aer_num > 0._r8 ) then
413 : molal = (1.e-6_r8*wact_factor/(mwso4*total_interstitial_aer_num*1.e6_r8))/ &
414 5081971924 : (4*pi/3*rhoh2o*(MAX(r3lx,4.e-6_r8))**3)
415 5081971924 : aw = 1._r8/(1._r8+2.9244948e-2_r8*molal+2.3141243e-3_r8*molal**2+7.8184854e-7_r8*molal**3)
416 : end if
417 :
418 : !*****************************************************************************
419 : ! immersion freezing begin
420 : !*****************************************************************************
421 :
422 5081971924 : frzimm = 0._r8
423 5081971924 : frzcnt = 0._r8
424 5081971924 : frzdep = 0._r8
425 :
426 : ! take solute effect into account
427 5081971924 : rgimm_aer = rgimm
428 :
429 : ! if aw*Si<=1, the freezing point depression is strong enough to prevent freezing
430 :
431 5081971924 : do_frz = aw*supersatice > 1._r8
432 5081971924 : if (do_frz) then
433 4800189048 : rgimm_aer = 2*vwice*sigma_iw/(boltz*T*LOG(aw*supersatice))
434 : else
435 281782876 : return
436 : endif
437 :
438 4800189048 : do_imm = T <= 263.15_r8 ! temperature threshold for immersion freezing (-10 C)
439 :
440 4800189048 : if (do_imm) then
441 : ! homogeneous energy of germ formation
442 3927906824 : dg0imm = 4*pi/3._r8*sigma_iw*rgimm_aer**2
443 :
444 : ! prefactor
445 3927906824 : Aimm = n1*((vwice*rhplanck)/(rgimm_aer**3)*SQRT(3._r8/pi*boltz*T*dg0imm))
446 :
447 : ! nucleation rate per particle
448 :
449 3927906824 : if (pdf_imm) then
450 1834215884 : dim_Jimm(:) = 0._r8
451 >11372*10^7 : do i = i1,i2
452 : ! 1/sqrt(f)
453 >11188*10^7 : dim_Jimm(i) = Aimm*mradius**2/SQRT(dim_f_imm(i))*EXP((-dga_imm-dim_f_imm(i)*dg0imm)/(boltz*T))
454 >11372*10^7 : dim_Jimm(i) = max(dim_Jimm(i), 0._r8)
455 : end do
456 :
457 : sum_imm = 0._r8
458 >11188*10^7 : do i = i1,i2-1
459 >11005*10^7 : sum_imm = sum_imm + 0.5_r8*((pdf_imm_theta(i )*exp(-dim_Jimm(i )*deltat)+ &
460 >22194*10^7 : pdf_imm_theta(i+1)*exp(-dim_Jimm(i+1)*deltat)))*pdf_d_theta
461 : end do
462 1834215884 : if (sum_imm > 0.99_r8) then
463 724871142 : sum_imm = 1.0_r8
464 : end if
465 : else
466 2093690940 : Jimm = Aimm*mradius**2/SQRT(f_imm)*EXP(( -dga_imm - f_imm*dg0imm )/(boltz*T))
467 2093690940 : sum_imm = exp(-Jimm*deltat)
468 : end if
469 : end if
470 :
471 : !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
472 : ! Deposition nucleation
473 : !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
474 :
475 : ! nucleation rate per particle
476 4800189048 : if (rgdep > 0) then
477 4800189048 : Jdep = Adep*mradius**2/SQRT(f_dep)*EXP((-dga_dep-f_dep*dg0dep)/(boltz*T))
478 : else
479 : Jdep = 0._r8
480 : end if
481 :
482 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
483 : ! contact nucleation
484 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
485 :
486 : ! nucleation rate per particle
487 4800189048 : Jcnt = Acnt*mradius**2*EXP((-dga_dep-f_cnt*dg0cnt)/(boltz*T))*Kcoll*icnlx
488 :
489 : ! Limit to 1% of available potential IN (for BC), no limit for dust
490 : if (tot_in) then
491 : if (do_imm) then
492 : frzimm = MIN(limfac*fn*total_aer_num/deltat, fn*total_aer_num/deltat*(1._r8-sum_imm))
493 : end if
494 : frzdep = MIN(limfac*(1._r8-fn)*(1._r8-dstcoat)*total_aer_num/deltat, &
495 : (1._r8-fn)*(1._r8-dstcoat)*total_aer_num/deltat*(1._r8-exp(-Jdep*deltat)))
496 : frzcnt = MIN(limfac*(1._r8-fn)*(1._r8-dstcoat)*total_aer_num/deltat, &
497 : (1._r8-fn)*(1._r8-dstcoat)*total_aer_num/deltat*(1._r8-exp(-Jcnt*deltat)))
498 : else
499 4800189048 : if (do_imm) then
500 3927906824 : frzimm = MIN(limfac*total_cloudborne_aer_num /deltat, total_cloudborne_aer_num/deltat*(1._r8-sum_imm))
501 : end if
502 4800189048 : frzdep = MIN(limfac*uncoated_aer_num/deltat, uncoated_aer_num/deltat*(1._r8-exp(-Jdep*deltat)))
503 4800189048 : frzcnt = MIN(limfac*uncoated_aer_num/deltat, uncoated_aer_num/deltat*(1._r8-exp(-Jcnt*deltat)))
504 : end if
505 :
506 4800189048 : if (frzcnt <= -1._r8) then
507 0 : write(iulog,*) 'hetfrz_classnuc_calc: frzcnt, Jcnt, Kcoll: ', frzcnt, Jcnt, Kcoll
508 0 : errstring = 'ERROR in hetfrz_classnuc_calc::frzducnt'
509 0 : return
510 : end if
511 :
512 5081971924 : end subroutine hetfrz_classnuc_calc_rates
513 :
514 : !===================================================================================================
515 :
516 : !-----------------------------------------------------------------------
517 : !
518 : ! Purpose: calculate collision kernels as a function of environmental parameters and aerosol/droplet sizes
519 : !
520 : ! Author: Corinna Hoose, UiO, October 2009
521 : !
522 : ! Modifications: Yong Wang and Xiaohong Liu, UWyo, 12/2012
523 : !
524 : ! "Seinfeld & Pandis" referenced in several places in this routine is:
525 : ! Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 3rd Edition
526 : ! John H. Seinfeld, Spyros N. Pandis ISBN: 978-1-118-94740-1
527 : !-----------------------------------------------------------------------
528 :
529 1270492981 : subroutine collkernel( temp, pres, eswtr, rhwincloud, r3lx, rad, Ktherm, Kcoll )
530 :
531 : real(r8), intent(in) :: temp ! temperature [K]
532 : real(r8), intent(in) :: pres ! pressure [Pa]
533 : real(r8), intent(in) :: eswtr ! saturation vapor pressure of water [Pa]
534 : real(r8), intent(in) :: r3lx ! volume mean drop radius [m]
535 : real(r8), intent(in) :: rhwincloud ! in-cloud relative humidity over water [ ]
536 : real(r8), intent(in) :: rad(:) ! aerosol radius [m]
537 : real(r8), intent(in) :: Ktherm(:) ! thermal conductivity of aerosol [J/(m s K)]
538 : real(r8), intent(out) :: Kcoll(:) ! collision kernel [cm3 s-1]
539 :
540 : ! local variables
541 : real(r8) :: a, b, c, a_f, b_f, c_f, f
542 : real(r8) :: tc ! temperature [deg C]
543 : real(r8) :: rho_air ! air density [kg m-3]
544 : real(r8) :: viscos_air ! dynamic viscosity of air [kg m-1 s-1]
545 : real(r8) :: Ktherm_air ! thermal conductivity of air [J/(m s K)]
546 : real(r8) :: lambda ! mean free path [m]
547 : real(r8) :: Kn ! Knudsen number [ ]
548 : real(r8) :: Re ! Reynolds number [ ]
549 : real(r8) :: Pr ! Prandtl number [ ]
550 : real(r8) :: Sc ! Schmidt number [ ]
551 : real(r8) :: vterm ! terminal velocity [m s-1]
552 : real(r8) :: Dvap ! water vapor diffusivity [m2 s-1]
553 : real(r8) :: Daer ! aerosol diffusivity [m2 s-1]
554 : real(r8) :: latvap ! latent heat of vaporization [J kg-1]
555 : real(r8) :: G ! thermodynamic function in Cotton et al. [kg m-1 s-1]
556 : real(r8) :: f_t ! factor by Waldmann & Schmidt [ ]
557 : real(r8) :: Q_heat ! heat flux [J m-2 s-1]
558 : real(r8) :: Tdiff_cotton ! temperature difference between droplet and environment [K]
559 : real(r8) :: K_brownian,K_thermo_cotton,K_diffusio_cotton ! collision kernels [m3 s-1]
560 :
561 : integer :: ntot, idx
562 :
563 : !------------------------------------------------------------------------------------------------
564 :
565 1270492981 : ntot = size(ktherm)
566 :
567 6352464905 : Kcoll(:) = 0._r8
568 :
569 1270492981 : tc = temp - tmelt
570 :
571 : ! air viscosity for tc<0, from depvel_part.F90
572 1270492981 : viscos_air = (1.718_r8+0.0049_r8*tc-1.2e-5_r8*tc*tc)*1.e-5_r8
573 : ! air density
574 1270492981 : rho_air = pres/(rair*temp)
575 : ! mean free path: Seinfeld & Pandis 8.6 (Book: ISBN: 978-1-118-94740-1)
576 1270492981 : lambda = 2*viscos_air/(pres*SQRT(8/(pi*rair*temp)))
577 : ! latent heat of vaporization, varies with T
578 1270492981 : latvap = 1000*(-0.0000614342_r8*tc**3 + 0.00158927_r8*tc**2 - 2.36418_r8*tc + 2500.79_r8)
579 : ! droplet terminal velocity after Chen & Liu, QJRMS 2004 (https://doi-org.cuucar.idm.oclc.org/10.1256/qj.03.41)
580 1270492981 : a = 8.8462e2_r8
581 1270492981 : b = 9.7593e7_r8
582 1270492981 : c = -3.4249e-11_r8
583 1270492981 : a_f = 3.1250e-1_r8
584 1270492981 : b_f = 1.0552e-3_r8
585 1270492981 : c_f = -2.4023_r8
586 1270492981 : f = EXP(EXP(a_f + b_f*(LOG(r3lx))**3 + c_f*rho_air**1.5_r8))
587 1270492981 : vterm = (a+ (b + c*r3lx)*r3lx)*r3lx*f
588 :
589 : ! Reynolds number
590 1270492981 : Re = 2*vterm*r3lx*rho_air/viscos_air
591 : ! thermal conductivity of air: Seinfeld & Pandis eq. 15.75 (Book: ISBN: 978-1-118-94740-1)
592 1270492981 : Ktherm_air = 1.e-3_r8*(4.39_r8+0.071_r8*temp) !J/(m s K)
593 : ! Prandtl number
594 1270492981 : Pr = viscos_air*cpair/Ktherm_air
595 : ! water vapor diffusivity: Pruppacher & Klett 13-3 (https://link.springer.com/book/10.1007/978-0-306-48100-0)
596 1270492981 : Dvap = 0.211e-4_r8*(temp/tmelt)*(pstd/pres)
597 : ! G-factor = rhoh2o*Xi in Rogers & Yau, p. 104
598 : G = rhoh2o/((latvap/(rh2o*temp) - 1)*latvap*rhoh2o/(Ktherm_air*temp) &
599 1270492981 : + rhoh2o*rh2o*temp/(Dvap*eswtr))
600 :
601 6352464905 : do idx = 1,ntot
602 6352464905 : if (rad(idx)>0._r8) then
603 : ! Knudsen number (Seinfeld & Pandis 8.1) (Book: ISBN: 978-1-118-94740-1)
604 5081971924 : Kn = lambda/rad(idx)
605 : ! aerosol diffusivity
606 5081971924 : Daer = boltz*temp*(1 + Kn)/(6*pi*rad(idx)*viscos_air)
607 :
608 : ! Schmidt number
609 5081971924 : Sc = viscos_air/(Daer*rho_air)
610 :
611 : ! Young (1974) first equ. on page 771 (doi: 10.1175/1520-0469(1974)031<0768:TROCNI>2.0.CO;2)
612 5081971924 : K_brownian = 4*pi*r3lx*Daer*(1 + 0.3_r8*Re**0.5_r8*Sc**0.33_r8)
613 :
614 : ! thermal conductivities from Seinfeld & Pandis, Table 8.6 (Book: ISBN: 978-1-118-94740-1)
615 : ! form factor
616 : f_t = 0.4_r8*(1._r8 + 1.45_r8*Kn + 0.4_r8*Kn*EXP(-1._r8/Kn)) &
617 5081971924 : *(Ktherm_air + 2.5_r8*Kn*Ktherm(idx)) &
618 5081971924 : /((1._r8 + 3._r8*Kn)*(2._r8*Ktherm_air + 5._r8*Kn*Ktherm(idx)+Ktherm(idx)))
619 :
620 : ! calculate T-Tc as in Cotton et al.
621 5081971924 : Tdiff_cotton = -G*(rhwincloud - 1._r8)*latvap/Ktherm_air
622 5081971924 : Q_heat = Ktherm_air/r3lx*(1._r8 + 0.3_r8*Re**0.5_r8*Pr**0.33_r8)*Tdiff_cotton
623 5081971924 : K_thermo_cotton = 4._r8*pi*r3lx*r3lx*f_t*Q_heat/pres
624 5081971924 : K_diffusio_cotton = -(1._r8/f_t)*(rh2o*temp/latvap)*K_thermo_cotton
625 5081971924 : Kcoll(idx) = 1.e6_r8*(K_brownian + K_thermo_cotton + K_diffusio_cotton) ! convert m3/s -> cm3/s
626 :
627 : ! set K to 0 if negative
628 5081971924 : if (Kcoll(idx) < 0._r8) Kcoll(idx) = 0._r8
629 : else
630 0 : Kcoll(idx) = 0._r8
631 : endif
632 : end do
633 :
634 1270492981 : end subroutine collkernel
635 :
636 : !===================================================================================================
637 :
638 1536 : subroutine hetfrz_classnuc_init_pdftheta()
639 :
640 : ! Local variables:
641 : real(r8) :: theta_min, theta_max
642 : real(r8) :: x1_imm, x2_imm
643 : real(r8) :: norm_theta_imm
644 : real(r8) :: imm_dust_mean_theta
645 : real(r8) :: imm_dust_var_theta
646 : integer :: i
647 : real(r8) :: m
648 : real(r8) :: temp
649 : !----------------------------------------------------------------------------
650 :
651 1536 : theta_min = pi/180._r8
652 1536 : theta_max = 179._r8/180._r8*pi
653 1536 : imm_dust_mean_theta = 46.0_r8/180.0_r8*pi
654 1536 : imm_dust_var_theta = 0.01_r8
655 :
656 1536 : pdf_d_theta = (179._r8-1._r8)/180._r8*pi/(pdf_n_theta-1)
657 :
658 1536 : x1_imm = (LOG(theta_min) - LOG(imm_dust_mean_theta))/(sqrt(2.0_r8)*imm_dust_var_theta)
659 1536 : x2_imm = (LOG(theta_max) - LOG(imm_dust_mean_theta))/(sqrt(2.0_r8)*imm_dust_var_theta)
660 1536 : norm_theta_imm = (ERF(x2_imm) - ERF(x1_imm))*0.5_r8
661 1536 : dim_theta = 0.0_r8
662 1536 : pdf_imm_theta = 0.0_r8
663 95232 : do i = i1, i2
664 93696 : dim_theta(i) = 1._r8/180._r8*pi + (i-1)*pdf_d_theta
665 : pdf_imm_theta(i) = exp(-((LOG(dim_theta(i)) - LOG(imm_dust_mean_theta))**2._r8) / &
666 : (2._r8*imm_dust_var_theta**2._r8) ) / &
667 95232 : (dim_theta(i)*imm_dust_var_theta*SQRT(2*pi))/norm_theta_imm
668 : end do
669 :
670 95232 : do i = i1, i2
671 93696 : m = cos(dim_theta(i))
672 93696 : temp = (2+m)*(1-m)**2/4._r8
673 95232 : dim_f_imm(i) = temp
674 : end do
675 :
676 1536 : end subroutine hetfrz_classnuc_init_pdftheta
677 :
678 : !===================================================================================================
679 :
680 : end module hetfrz_classnuc
|
