Line data Source code
1 : module nucleate_ice
2 :
3 : !-------------------------------------------------------------------------------
4 : ! Purpose:
5 : ! A parameterization of ice nucleation.
6 : !
7 : ! *** This module is intended to be a "portable" code layer. Ideally it should
8 : ! *** not contain any use association of modules that belong to the model framework.
9 : !
10 : !
11 : ! Method:
12 : ! The current method is based on Liu & Penner (2005) & Liu et al. (2007)
13 : ! It related the ice nucleation with the aerosol number, temperature and the
14 : ! updraft velocity. It includes homogeneous freezing of sulfate & immersion
15 : ! freezing on mineral dust (soot disabled) in cirrus clouds, and
16 : ! Meyers et al. (1992) deposition nucleation in mixed-phase clouds
17 : !
18 : ! The effect of preexisting ice crystals on ice nucleation in cirrus clouds is included,
19 : ! and also consider the sub-grid variability of temperature in cirrus clouds,
20 : ! following X. Shi et al. ACP (2014).
21 : !
22 : ! Ice nucleation in mixed-phase clouds now uses classical nucleation theory (CNT),
23 : ! follows Y. Wang et al. ACP (2014), Hoose et al. (2010).
24 : !
25 : ! Authors:
26 : ! Xiaohong Liu, 01/2005, modifications by A. Gettelman 2009-2010
27 : ! Xiangjun Shi & Xiaohong Liu, 01/2014.
28 : !
29 : ! With help from C. C. Chen and B. Eaton (2014)
30 : !-------------------------------------------------------------------------------
31 :
32 : use wv_saturation, only: svp_water, svp_ice
33 :
34 : implicit none
35 : private
36 : save
37 :
38 : integer, parameter :: r8 = selected_real_kind(12)
39 :
40 : public :: nucleati_init, nucleati
41 :
42 : logical :: use_preexisting_ice
43 : logical :: use_hetfrz_classnuc
44 : logical :: use_incloud_nuc
45 : integer :: iulog
46 : real(r8) :: pi
47 : real(r8) :: mincld
48 :
49 : real(r8), parameter :: Shet = 1.3_r8 ! het freezing threshold
50 : real(r8), parameter :: rhoice = 0.5e3_r8 ! kg/m3, Wpice is not sensitive to rhoice
51 : real(r8), parameter :: minweff= 0.001_r8 ! m/s
52 : real(r8), parameter :: gamma1=1.0_r8
53 : real(r8), parameter :: gamma2=1.0_r8
54 : real(r8), parameter :: gamma3=2.0_r8
55 : real(r8), parameter :: gamma4=6.0_r8
56 :
57 : real(r8) :: ci
58 :
59 : !===============================================================================
60 : contains
61 : !===============================================================================
62 :
63 1536 : subroutine nucleati_init( &
64 : use_preexisting_ice_in, use_hetfrz_classnuc_in, use_incloud_nuc_in, iulog_in, pi_in, &
65 : mincld_in)
66 :
67 : logical, intent(in) :: use_preexisting_ice_in
68 : logical, intent(in) :: use_hetfrz_classnuc_in
69 : logical, intent(in) :: use_incloud_nuc_in
70 : integer, intent(in) :: iulog_in
71 : real(r8), intent(in) :: pi_in
72 : real(r8), intent(in) :: mincld_in
73 :
74 1536 : use_preexisting_ice = use_preexisting_ice_in
75 1536 : use_hetfrz_classnuc = use_hetfrz_classnuc_in
76 1536 : use_incloud_nuc = use_incloud_nuc_in
77 1536 : iulog = iulog_in
78 1536 : pi = pi_in
79 1536 : mincld = mincld_in
80 :
81 1536 : ci = rhoice*pi/6._r8
82 :
83 1536 : end subroutine nucleati_init
84 :
85 : !===============================================================================
86 :
87 87518271 : subroutine nucleati( &
88 : wbar, tair, pmid, relhum, cldn, &
89 : qc, qi, ni_in, rhoair, &
90 : so4_num, dst_num, soot_num, subgrid, &
91 : nuci, onihf, oniimm, onidep, onimey, &
92 : wpice, weff, fhom, regm, &
93 : oso4_num, odst_num, osoot_num, &
94 : call_frm_zm_in, add_preexisting_ice_in)
95 :
96 : ! Input Arguments
97 : real(r8), intent(in) :: wbar ! grid cell mean vertical velocity (m/s)
98 : real(r8), intent(in) :: tair ! temperature (K)
99 : real(r8), intent(in) :: pmid ! pressure at layer midpoints (pa)
100 : real(r8), intent(in) :: relhum ! relative humidity with respective to liquid
101 : real(r8), intent(in) :: cldn ! new value of cloud fraction (fraction)
102 : real(r8), intent(in) :: qc ! liquid water mixing ratio (kg/kg)
103 : real(r8), intent(in) :: qi ! grid-mean preexisting cloud ice mass mixing ratio (kg/kg)
104 : real(r8), intent(in) :: ni_in ! grid-mean preexisting cloud ice number conc (#/kg)
105 : real(r8), intent(in) :: rhoair ! air density (kg/m3)
106 : real(r8), intent(in) :: so4_num ! so4 aerosol number (#/cm^3)
107 : real(r8), intent(in) :: dst_num ! total dust aerosol number (#/cm^3)
108 : real(r8), intent(in) :: soot_num ! soot (hydrophilic) aerosol number (#/cm^3)
109 : real(r8), intent(in) :: subgrid ! subgrid saturation scaling factor
110 :
111 : ! Output Arguments
112 : real(r8), intent(out) :: nuci ! ice number nucleated (#/kg)
113 : real(r8), intent(out) :: onihf ! nucleated number from homogeneous freezing of so4
114 : real(r8), intent(out) :: oniimm ! nucleated number from immersion freezing
115 : real(r8), intent(out) :: onidep ! nucleated number from deposition nucleation
116 : real(r8), intent(out) :: onimey ! nucleated number from deposition nucleation (meyers: mixed phase)
117 : real(r8), intent(out) :: wpice ! diagnosed Vertical velocity Reduction caused by preexisting ice (m/s), at Shom
118 : real(r8), intent(out) :: weff ! effective Vertical velocity for ice nucleation (m/s); weff=wbar-wpice
119 : real(r8), intent(out) :: fhom ! how much fraction of cloud can reach Shom
120 : real(r8), intent(out) :: regm ! nucleation regime indiator
121 : real(r8), intent(out) :: oso4_num ! so4 aerosol number (#/cm^3)
122 : real(r8), intent(out) :: odst_num ! total dust aerosol number (#/cm^3)
123 : real(r8), intent(out) :: osoot_num ! soot (hydrophilic) aerosol number (#/cm^3)
124 :
125 : ! Optional Arguments
126 : logical, intent(in), optional :: call_frm_zm_in ! true if called from ZM convection scheme
127 : logical, intent(in), optional :: add_preexisting_ice_in ! only false if called with pumas_v1.21+
128 :
129 : ! Local workspace
130 : real(r8) :: nihf ! nucleated number from homogeneous freezing of so4
131 : real(r8) :: niimm ! nucleated number from immersion freezing
132 : real(r8) :: nidep ! nucleated number from deposition nucleation
133 : real(r8) :: nimey ! nucleated number from deposition nucleation (meyers)
134 : real(r8) :: n1, ni ! nucleated number
135 : real(r8) :: tc, A, B ! work variable
136 : real(r8) :: esl, esi, deles ! work variable
137 : real(r8) :: wbar1, wbar2
138 :
139 : ! used in SUBROUTINE Vpreice
140 : real(r8) :: Ni_preice ! cloud ice number conc (1/m3)
141 : real(r8) :: lami,Ri_preice ! mean cloud ice radius (m)
142 : real(r8) :: Shom ! initial ice saturation ratio; if <1, use hom threshold Si
143 : real(r8) :: detaT,RHimean ! temperature standard deviation, mean cloudy RHi
144 : real(r8) :: wpicehet ! diagnosed Vertical velocity Reduction caused by preexisting ice (m/s), at shet
145 :
146 : real(r8) :: weffhet ! effective Vertical velocity for ice nucleation (m/s) weff=wbar-wpicehet
147 :
148 : logical :: call_frm_zm, add_preexisting_ice
149 : !-------------------------------------------------------------------------------
150 :
151 87518271 : nuci = 0._r8
152 :
153 87518271 : RHimean = relhum*svp_water(tair)/svp_ice(tair)*subgrid
154 :
155 : ! temp variables that depend on use_preexisting_ice
156 87518271 : wbar1 = wbar
157 87518271 : wbar2 = wbar
158 :
159 : ! If not using prexisting ice, the homogeneous freezing happens in the
160 : ! entire gridbox.
161 87518271 : fhom = 1._r8
162 :
163 87518271 : if (present(call_frm_zm_in)) then
164 0 : call_frm_zm = call_frm_zm_in
165 : else
166 : call_frm_zm = .false.
167 : end if
168 :
169 87518271 : if (present(add_preexisting_ice_in)) then
170 0 : add_preexisting_ice = add_preexisting_ice_in
171 : else
172 : add_preexisting_ice = .true.
173 : end if
174 :
175 87518271 : if (use_preexisting_ice .and. (.not. call_frm_zm)) then
176 :
177 87518271 : Ni_preice = ni_in*rhoair ! (convert from #/kg -> #/m3)
178 87518271 : Ni_preice = Ni_preice / max(mincld,cldn) ! in-cloud ice number density
179 :
180 87518271 : if (Ni_preice > 10.0_r8 .and. qi > 1.e-10_r8) then ! > 0.01/L = 10/m3
181 20285357 : Shom = -1.5_r8 ! if Shom<1 , Shom will be recalculated in SUBROUTINE Vpreice, according to Ren & McKenzie, 2005
182 20285357 : lami = (gamma4*ci*ni_in/qi)**(1._r8/3._r8)
183 20285357 : Ri_preice = 0.5_r8/lami ! radius
184 20285357 : Ri_preice = max(Ri_preice, 1e-8_r8) ! >0.01micron
185 20285357 : call Vpreice(pmid, tair, Ri_preice, Ni_preice, Shom, wpice)
186 20285357 : call Vpreice(pmid, tair, Ri_preice, Ni_preice, Shet, wpicehet)
187 : else
188 67232914 : wpice = 0.0_r8
189 67232914 : wpicehet = 0.0_r8
190 : endif
191 :
192 87518271 : weff = max(wbar-wpice, minweff)
193 87518271 : wpice = min(wpice, wbar)
194 87518271 : weffhet = max(wbar-wpicehet,minweff)
195 87518271 : wpicehet = min(wpicehet, wbar)
196 :
197 87518271 : wbar1 = weff
198 87518271 : wbar2 = weffhet
199 :
200 87518271 : detaT = wbar/0.23_r8
201 87518271 : if (use_incloud_nuc) then
202 0 : call frachom(tair, 1._r8, detaT, fhom)
203 : else
204 87518271 : call frachom(tair, RHimean, detaT, fhom)
205 : end if
206 : end if
207 :
208 87518271 : ni = 0._r8
209 87518271 : tc = tair - 273.15_r8
210 :
211 : ! initialize
212 87518271 : niimm = 0._r8
213 87518271 : nidep = 0._r8
214 87518271 : nihf = 0._r8
215 87518271 : deles = 0._r8
216 87518271 : esi = 0._r8
217 87518271 : regm = 0._r8
218 :
219 87518271 : oso4_num = 0._r8
220 87518271 : odst_num = 0._r8
221 87518271 : osoot_num = 0._r8
222 :
223 87518271 : if ((so4_num >= 1.0e-10_r8 .or. (soot_num+dst_num) >= 1.0e-10_r8) .and. cldn > 0._r8) then
224 :
225 87021369 : if (RHimean.ge.1.2_r8) then
226 :
227 11807471 : if ( ((tc.le.0.0_r8).and.(tc.ge.-37.0_r8).and.(qc.lt.1.e-12_r8)).or.(tc.le.-37.0_r8)) then
228 :
229 7340226 : if ( (soot_num+dst_num) > 0._r8) then
230 7338104 : A = -1.4938_r8 * log(soot_num+dst_num) + 12.884_r8
231 7338104 : B = -10.41_r8 * log(soot_num+dst_num) - 67.69_r8
232 7338104 : regm = A * log(wbar1) + B
233 : end if
234 :
235 : ! heterogeneous nucleation only
236 7340226 : if (tc .gt. regm .or. so4_num < 1.0e-10_r8) then
237 :
238 5084045 : if(tc.lt.-40._r8 .and. wbar1.gt.1._r8 .and. so4_num >= 1.0e-10_r8) then ! exclude T<-40 & W>1m/s from hetero. nucleation
239 :
240 0 : call hf(tc,wbar1,relhum*subgrid,so4_num,nihf)
241 0 : niimm=0._r8
242 0 : nidep=0._r8
243 :
244 : ! If some homogeneous nucleation happened, assume all of the that heterogeneous
245 : ! and coarse mode sulfate particles nucleated.
246 0 : if (nihf.gt.1e-3_r8) then ! hom occur, add preexisting ice
247 0 : niimm = dst_num + soot_num ! assuming dst_num freeze firstly
248 0 : odst_num = dst_num
249 0 : osoot_num = soot_num
250 :
251 0 : oso4_num = nihf
252 : endif
253 :
254 0 : nihf = nihf * fhom
255 0 : oso4_num = oso4_num * fhom
256 :
257 0 : n1 = nihf + niimm
258 : else
259 :
260 5084045 : call hetero(tc,wbar2,soot_num+dst_num,niimm,nidep)
261 :
262 5084045 : nihf = 0._r8
263 5084045 : n1 = niimm + nidep
264 :
265 5084045 : if ( (soot_num+dst_num) > 0._r8) then
266 5084045 : osoot_num = soot_num * (niimm + nidep) / (soot_num + dst_num)
267 5084045 : odst_num = dst_num * (niimm + nidep) / (soot_num + dst_num)
268 : end if
269 :
270 : endif
271 :
272 : ! homogeneous nucleation only
273 2256181 : else if (tc.lt.regm-5._r8 .or. (soot_num+dst_num) < 1.0e-10_r8) then
274 :
275 2071144 : call hf(tc,wbar1,relhum*subgrid,so4_num,nihf)
276 2071144 : niimm=0._r8
277 2071144 : nidep=0._r8
278 :
279 : ! If some homogeneous nucleation happened, assume all of the that
280 : ! heterogeneous and coarse mode sulfate particles nucleated.
281 2071144 : if (nihf.gt.1e-3_r8) then ! hom occur, add preexisting ice
282 9459 : niimm = dst_num + soot_num ! assuming dst_num freeze firstly
283 9459 : odst_num = dst_num
284 9459 : osoot_num = soot_num
285 :
286 9459 : oso4_num = nihf
287 : endif
288 :
289 2071144 : nihf = nihf * fhom
290 2071144 : oso4_num = oso4_num * fhom
291 :
292 2071144 : n1 = nihf + niimm
293 :
294 : ! transition between homogeneous and heterogeneous: interpolate in-between
295 : else
296 :
297 185037 : if (tc.lt.-40._r8 .and. wbar1.gt.1._r8) then ! exclude T<-40 & W>1m/s from hetero. nucleation
298 :
299 0 : call hf(tc, wbar1, relhum*subgrid, so4_num, nihf)
300 0 : niimm = 0._r8
301 0 : nidep = 0._r8
302 :
303 : ! If some homogeneous nucleation happened, assume all of the
304 : ! that heterogeneous and coarse mode sulfate particles nucleated.
305 0 : if (nihf.gt.1e-3_r8) then ! hom occur, add preexisting ice
306 0 : niimm = dst_num + soot_num ! assuming dst_num freeze firstly
307 0 : odst_num = dst_num
308 0 : osoot_num = soot_num
309 :
310 0 : oso4_num = nihf
311 : endif
312 :
313 0 : nihf = nihf * fhom
314 0 : oso4_num = oso4_num * fhom
315 :
316 0 : n1 = nihf + niimm
317 :
318 : else
319 :
320 185037 : call hf(regm-5._r8,wbar1,relhum*subgrid,so4_num,nihf)
321 185037 : call hetero(regm,wbar2,soot_num+dst_num,niimm,nidep)
322 :
323 : ! If some homogeneous nucleation happened, assume all of the
324 : ! heterogeneous particles nucleated and add in a fraction of
325 : ! the homogeneous freezing.
326 185037 : if (nihf.gt.1e-3_r8) then ! hom occur, add preexisting ice
327 4866 : oso4_num = nihf
328 : endif
329 :
330 185037 : if ( (soot_num+dst_num) > 0._r8) then
331 185037 : osoot_num = soot_num * (niimm + nidep) / (soot_num + dst_num)
332 185037 : odst_num = dst_num * (niimm + nidep) / (soot_num + dst_num)
333 : end if
334 :
335 185037 : nihf = nihf * fhom * ((regm - tc) / 5._r8)**2
336 185037 : oso4_num = oso4_num * fhom * ((regm - tc) / 5._r8)**2
337 :
338 185037 : n1 = niimm + nidep + nihf
339 :
340 : end if
341 : end if
342 :
343 : ! Scale the rates for in-cloud number, since this is what
344 : ! MG is expecting to find.
345 7340226 : ni = n1
346 :
347 : ! If using prexsiting ice, and allowed to add, then add it to the total.
348 7340226 : if (use_preexisting_ice) then
349 7340226 : if (add_preexisting_ice .and. (.not. call_frm_zm)) then
350 7340226 : ni = ni + Ni_preice * 1e-6_r8
351 : end if
352 : end if
353 : end if
354 : end if
355 : end if
356 :
357 : ! deposition/condensation nucleation in mixed clouds (-37<T<0C) (Meyers, 1992)
358 87518271 : if(tc.lt.0._r8 .and. tc.gt.-37._r8 .and. qc.gt.1.e-12_r8) then
359 5562797 : esl = svp_water(tair) ! over water in mixed clouds
360 5562797 : esi = svp_ice(tair) ! over ice
361 5562797 : deles = (esl - esi)
362 5562797 : nimey=1.e-3_r8*exp(12.96_r8*deles/esi - 0.639_r8)
363 : else
364 : nimey=0._r8
365 : endif
366 :
367 87518271 : if (use_hetfrz_classnuc) nimey = 0._r8
368 :
369 87518271 : nuci=ni + nimey
370 :
371 87518271 : if(nuci.gt.9999._r8.or.nuci.lt.0._r8) then
372 0 : write(iulog, *) 'Warning: incorrect ice nucleation number (nuci reset =0)'
373 0 : write(iulog, *) ni, tair, relhum, wbar, nihf, niimm, nidep,deles,esi,dst_num,so4_num
374 0 : nuci=0._r8
375 : endif
376 :
377 87518271 : nuci = nuci*1.e+6_r8/rhoair ! change unit from #/cm3 to #/kg
378 87518271 : onimey = nimey*1.e+6_r8/rhoair
379 87518271 : onidep = nidep*1.e+6_r8/rhoair
380 87518271 : oniimm = niimm*1.e+6_r8/rhoair
381 87518271 : onihf = nihf*1.e+6_r8/rhoair
382 87518271 : end subroutine nucleati
383 :
384 : !===============================================================================
385 :
386 5269082 : subroutine hetero(T,ww,Ns,Nis,Nid)
387 :
388 : real(r8), intent(in) :: T, ww, Ns
389 : real(r8), intent(out) :: Nis, Nid
390 :
391 : real(r8) A11,A12,A21,A22,B11,B12,B21,B22
392 : real(r8) B,C
393 :
394 : !---------------------------------------------------------------------
395 : ! parameters
396 :
397 5269082 : A11 = 0.0263_r8
398 5269082 : A12 = -0.0185_r8
399 5269082 : A21 = 2.758_r8
400 5269082 : A22 = 1.3221_r8
401 5269082 : B11 = -0.008_r8
402 5269082 : B12 = -0.0468_r8
403 5269082 : B21 = -0.2667_r8
404 5269082 : B22 = -1.4588_r8
405 :
406 : ! ice from immersion nucleation (cm^-3)
407 :
408 5269082 : B = (A11+B11*log(Ns)) * log(ww) + (A12+B12*log(Ns))
409 5269082 : C = A21+B21*log(Ns)
410 :
411 5269082 : Nis = exp(A22) * Ns**B22 * exp(B*T) * ww**C
412 5269082 : Nis = min(Nis,Ns)
413 :
414 5269082 : Nid = 0.0_r8 ! don't include deposition nucleation for cirrus clouds when T<-37C
415 :
416 5269082 : end subroutine hetero
417 :
418 : !===============================================================================
419 :
420 2256181 : subroutine hf(T,ww,RH,Na,Ni)
421 :
422 : real(r8), intent(in) :: T, ww, RH, Na
423 : real(r8), intent(out) :: Ni
424 :
425 : real(r8) A1_fast,A21_fast,A22_fast,B1_fast,B21_fast,B22_fast
426 : real(r8) A2_fast,B2_fast
427 : real(r8) C1_fast,C2_fast,k1_fast,k2_fast
428 : real(r8) A1_slow,A2_slow,B1_slow,B2_slow,B3_slow
429 : real(r8) C1_slow,C2_slow,k1_slow,k2_slow
430 : real(r8) regm
431 : real(r8) A,B,C
432 : real(r8) RHw
433 :
434 : !---------------------------------------------------------------------
435 : ! parameters
436 :
437 2256181 : A1_fast =0.0231_r8
438 2256181 : A21_fast =-1.6387_r8 !(T>-64 deg)
439 2256181 : A22_fast =-6.045_r8 !(T<=-64 deg)
440 2256181 : B1_fast =-0.008_r8
441 2256181 : B21_fast =-0.042_r8 !(T>-64 deg)
442 2256181 : B22_fast =-0.112_r8 !(T<=-64 deg)
443 2256181 : C1_fast =0.0739_r8
444 2256181 : C2_fast =1.2372_r8
445 :
446 2256181 : A1_slow =-0.3949_r8
447 2256181 : A2_slow =1.282_r8
448 2256181 : B1_slow =-0.0156_r8
449 2256181 : B2_slow =0.0111_r8
450 2256181 : B3_slow =0.0217_r8
451 2256181 : C1_slow =0.120_r8
452 2256181 : C2_slow =2.312_r8
453 :
454 2256181 : Ni = 0.0_r8
455 :
456 : !----------------------------
457 : !RHw parameters
458 2256181 : A = 6.0e-4_r8*log(ww)+6.6e-3_r8
459 2256181 : B = 6.0e-2_r8*log(ww)+1.052_r8
460 2256181 : C = 1.68_r8 *log(ww)+129.35_r8
461 2256181 : RHw=(A*T*T+B*T+C)*0.01_r8
462 :
463 2256181 : if((T.le.-37.0_r8) .and. ((RH).ge.RHw)) then
464 :
465 15076 : regm = 6.07_r8*log(ww)-55.0_r8
466 :
467 15076 : if(T.ge.regm) then ! fast-growth regime
468 :
469 12734 : if(T.gt.-64.0_r8) then
470 : A2_fast=A21_fast
471 : B2_fast=B21_fast
472 : else
473 675 : A2_fast=A22_fast
474 675 : B2_fast=B22_fast
475 : endif
476 :
477 12734 : k1_fast = exp(A2_fast + B2_fast*T + C2_fast*log(ww))
478 12734 : k2_fast = A1_fast+B1_fast*T+C1_fast*log(ww)
479 :
480 12734 : Ni = k1_fast*Na**(k2_fast)
481 12734 : Ni = min(Ni,Na)
482 :
483 : else ! slow-growth regime
484 :
485 2342 : k1_slow = exp(A2_slow + (B2_slow+B3_slow*log(ww))*T + C2_slow*log(ww))
486 2342 : k2_slow = A1_slow+B1_slow*T+C1_slow*log(ww)
487 :
488 2342 : Ni = k1_slow*Na**(k2_slow)
489 2342 : Ni = min(Ni,Na)
490 :
491 : endif
492 :
493 : end if
494 :
495 2256181 : end subroutine hf
496 :
497 : !===============================================================================
498 :
499 40570714 : SUBROUTINE Vpreice(P_in, T_in, R_in, C_in, S_in, V_out)
500 :
501 : ! based on Karcher et al. (2006)
502 : ! VERTICAL VELOCITY CALCULATED FROM DEPOSITIONAL LOSS TERM
503 :
504 : ! SUBROUTINE arguments
505 : REAL(r8), INTENT(in) :: P_in ! [Pa],INITIAL AIR pressure
506 : REAL(r8), INTENT(in) :: T_in ! [K] ,INITIAL AIR temperature
507 : REAL(r8), INTENT(in) :: R_in ! [m],INITIAL MEAN ICE CRYSTAL NUMBER RADIUS
508 : REAL(r8), INTENT(in) :: C_in ! [m-3],INITIAL TOTAL ICE CRYSTAL NUMBER DENSITY, [1/cm3]
509 : REAL(r8), INTENT(in) :: S_in ! [-],INITIAL ICE SATURATION RATIO;; if <1, use hom threshold Si
510 : REAL(r8), INTENT(out) :: V_out ! [m/s], VERTICAL VELOCITY REDUCTION (caused by preexisting ice)
511 :
512 : ! SUBROUTINE parameters
513 : REAL(r8), PARAMETER :: ALPHAc = 0.5_r8 ! density of ice (g/cm3), !!!V is not related to ALPHAc
514 : REAL(r8), PARAMETER :: FA1c = 0.601272523_r8
515 : REAL(r8), PARAMETER :: FA2c = 0.000342181855_r8
516 : REAL(r8), PARAMETER :: FA3c = 1.49236645E-12_r8
517 : REAL(r8), PARAMETER :: WVP1c = 3.6E+10_r8
518 : REAL(r8), PARAMETER :: WVP2c = 6145.0_r8
519 : REAL(r8), PARAMETER :: FVTHc = 11713803.0_r8
520 : REAL(r8), PARAMETER :: THOUBKc = 7.24637701E+18_r8
521 : REAL(r8), PARAMETER :: SVOLc = 3.23E-23_r8 ! SVOL=XMW/RHOICE
522 : REAL(r8), PARAMETER :: FDc = 249.239822_r8
523 : REAL(r8), PARAMETER :: FPIVOLc = 3.89051704E+23_r8
524 : REAL(r8) :: T,P,S,R,C
525 : REAL(r8) :: A1,A2,A3,B1,B2
526 : REAL(r8) :: T_1,PICE,FLUX,ALP4,CISAT,DLOSS,VICE
527 :
528 40570714 : T = T_in ! K , K
529 40570714 : P = P_in*1e-2_r8 ! Pa , hpa
530 :
531 40570714 : IF (S_in.LT.1.0_r8) THEN
532 20285357 : S = 2.349_r8 - (T/259.0_r8) ! homogeneous freezing threshold, according to Ren & McKenzie, 2005
533 : ELSE
534 : S = S_in ! INPUT ICE SATURATION RATIO, -, >1
535 : ENDIF
536 :
537 40570714 : R = R_in*1e2_r8 ! m => cm
538 40570714 : C = C_in*1e-6_r8 ! m-3 => cm-3
539 40570714 : T_1 = 1.0_r8/ T
540 40570714 : PICE = WVP1c * EXP(-(WVP2c*T_1))
541 40570714 : ALP4 = 0.25_r8 * ALPHAc
542 40570714 : FLUX = ALP4 * SQRT(FVTHc*T)
543 40570714 : CISAT = THOUBKc * PICE * T_1
544 40570714 : A1 = ( FA1c * T_1 - FA2c ) * T_1
545 40570714 : A2 = 1.0_r8/ CISAT
546 40570714 : A3 = FA3c * T_1 / P
547 40570714 : B1 = FLUX * SVOLc * CISAT * ( S-1.0_r8 )
548 40570714 : B2 = FLUX * FDc * P * T_1**1.94_r8
549 40570714 : DLOSS = FPIVOLc * C * B1 * R**2 / ( 1.0_r8+ B2 * R )
550 40570714 : VICE = ( A2 + A3 * S ) * DLOSS / ( A1 * S ) ! 2006,(19)
551 40570714 : V_out = VICE*1e-2_r8 ! cm/s => m/s
552 :
553 40570714 : END SUBROUTINE Vpreice
554 :
555 87518271 : subroutine frachom(Tmean,RHimean,detaT,fhom)
556 : ! How much fraction of cirrus might reach Shom
557 : ! base on "A cirrus cloud scheme for general circulation models",
558 : ! B. Karcher and U. Burkhardt 2008
559 :
560 : real(r8), intent(in) :: Tmean, RHimean, detaT
561 : real(r8), intent(out) :: fhom
562 :
563 : real(r8), parameter :: seta = 6132.9_r8 ! K
564 : integer, parameter :: Nbin=200 ! (Tmean - 3*detaT, Tmean + 3*detaT)
565 :
566 : real(r8) :: PDF_T(Nbin) ! temperature PDF; ! PDF_T=0 outside (Tmean-3*detaT, Tmean+3*detaT)
567 : real(r8) :: Sbin(Nbin) ! the fluctuations of Si that are driven by the T variations
568 : real(r8) :: Sihom, deta
569 : integer :: i
570 :
571 87518271 : Sihom = 2.349_r8-Tmean/259.0_r8 ! homogeneous freezing threshold, according to Ren & McKenzie, 2005
572 87518271 : fhom = 0.0_r8
573 :
574 907947389 : do i = Nbin, 1, -1
575 :
576 906495936 : deta = (i - 0.5_r8 - Nbin/2)*6.0_r8/Nbin ! PDF_T=0 outside (Tmean-3*detaT, Tmean+3*detaT)
577 906495936 : Sbin(i) = RHimean*exp(deta*detaT*seta/Tmean**2.0_r8)
578 906495936 : PDF_T(i) = exp(-deta**2.0_r8/2.0_r8)*6.0_r8/(sqrt(2.0_r8*Pi)*Nbin)
579 :
580 :
581 907947389 : if (Sbin(i).ge.Sihom) then
582 820429118 : fhom = fhom + PDF_T(i)
583 : else
584 : exit
585 : end if
586 : end do
587 :
588 87518271 : fhom = min(1.0_r8, fhom/0.997_r8) ! accounting for the finite limits (-3 , 3)
589 87518271 : end subroutine frachom
590 :
591 : end module nucleate_ice
|
