CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
aer_cloud.F
1
6
13 MODULE aer_cloud
14
15#ifdef GEOS5
16 use mapl_constantsmod, r8 => mapl_r8
17#endif
18#ifdef NEMS_GSM
19 use physcons, only : mapl_pi => con_pi
20 use machine, only : r8 => kind_phys
21#endif
22
23! according to the models of Nenes & Seinfeld (2003), Fountoukis and Nenes (2005) and Barahona and Nenes (2008, 2009).
24! *** Code Developer: Donifan Barahona donifan.o.barahona@nasa.gov
25!
26!=======================================================================
27!
28
29 implicit none
30 private
31
32 public :: aerosol_activate
33 public :: aerconversion
34 public :: aerconversion1
35 public :: aerprops
36 public :: getinsubset
37 public :: init_aer
38 public :: aer_cloud_init
39 public :: vertical_vel_variance
40
41
42 integer, parameter :: nsmx_par=20, npgauss=10
43
44
45
46 type :: aerprops
47! sequence
48 real, dimension(nsmx_par) :: num, dpg, sig, den, kap
49 &, fdust, fsoot, forg
50 integer :: nmods
51 end type aerprops
52
53 interface assignment (=)
54 module procedure copy_aer
55 end interface
56
57
58!==================================================================
59
60
61!
62!=================================================================
63
64!
65 real*8, dimension(npgauss) :: xgs_par, wgs_par
66!
67
68!Global aux variables
69
70 type(aerprops) :: aerpr_base_clean, aerpr_base_polluted
71
72 real*8 :: base_mass_so4_polluted, base_mass_so4_clean,
73 & base_mass_ss, frac_dust(5), frac_bc, frac_org,
74 & ahet_dust(5), ahet_bc
75
76!==================================================================
77
78!
79!=================================================================
80
81 real*8 :: sh_ice, doin_ice,vmin_ice, acorr_dust, acorr_bc
82 &, denw_par, cpair_par, dhv_par
83
84 integer :: typeofspec_ice
85 logical :: purehet_ice, purehom_ice, is_gocart
86
87!==================================================================
88!
89!==========================
90
91 integer, parameter :: maxit_par=100
92
93 real, parameter :: amw_par=18d-3, ama_par=29d-3
94 &, grav_par=9.81d0, rgas_par=8.31d0
95 &, accom_par=1.0d0, eps_par=1d-6
96 &, zero_par=1.0e-20, great_par=1d20
97 &, pi_par=mapl_pi, sq2pi_par=sqrt(pi_par)
98! &, pi_par=3.1415927d0, sq2pi_par=sqrt(pi_par)
99 &, sq2_par=1.41421356237d0
100!
101 &, wmw_ice=018d0, amw_ice=0.029d0
102 &, rgas_ice=8.314d0, grav_ice=9.81d0
103 &, cpa_ice=1005.1d0, pi_ice=pi_par
104 &, depcoef_ice=0.1d0, thaccom_ice=0.7d0
105!
106 &, to_ice=272.15d0, tmin_ice=185.d0
107 &, pmin_ice=100.0d0, thom=236.0d0
108 &, rv_ice=rgas_ice/wmw_ice
109
110
111 CONTAINS
112
115 subroutine aer_cloud_init()
116
117 real*8 :: daux, sigaux
118 integer ::ix
119
121
122 acorr_dust = 2.7e7
123 acorr_bc = 8.0e7
124
125 do ix = 1, 5
126 daux = aerpr_base_polluted%dpg(ix)
127 sigaux = aerpr_base_polluted%sig(ix)
128 frac_dust(ix) = 0.5d0*(1d0
129 & - erfapp(log(0.1e-6/daux)/(sigaux*sq2_par)))
130
131 ahet_dust(ix) = daux*daux*daux*0.52*acorr_dust
132 & * exp(4.5*sigaux*sigaux)
133
134 end do
135
136
137 daux = aerpr_base_polluted%dpg(12)
138 sigaux = aerpr_base_polluted%sig(12)
139 frac_bc = 0.5d0*(1d0-erfapp(log(0.1e-6/daux)/(sigaux*sq2_par)))
140 ahet_bc = daux*daux*daux*0.52*acorr_bc* exp(4.5*sigaux*sigaux)
141
142 daux = aerpr_base_polluted%dpg(13)
143 sigaux = aerpr_base_polluted%sig(13)
144 frac_org = 0.5d0*(1d0-erfapp(log(0.1e-6/daux)/(sigaux*sq2_par)))
145
146 end subroutine aer_cloud_init
147
148
185 subroutine aerosol_activate(tparc_in, pparc_in, sigwparc_in, &
186 & wparc_ls, Aer_Props, npre_in, dpre_in, ccn_diagr8, Ndropr8, &
187 & cdncr8, smaxliqr8, incr8, smaxicer8, nheticer8, INimmr8, &
188 & dINimmr8, Ncdepr8, Ncdhfr8, sc_icer8, fdust_immr8, fdust_depr8, &
189 & fdust_dhfr8, nlimr8, use_average_v, CCN_param, IN_param, fd_dust,&
190 & fd_soot, pfrz_inc_r8, sigma_nuc, rhi_cell,nccn)
191! & fd_soot, pfrz_inc_r8, sigma_nuc, rhi_cell,nccn, lprnt)
192
193
194
195
196 type(aerprops), intent(in) :: aer_props
197
198 logical :: use_average_v
199! logical :: use_average_v, lprnt
200
201 real(r8), intent(in) :: tparc_in, pparc_in, sigwparc_in, &
202 & wparc_ls, npre_in, dpre_in, Ndropr8, fd_soot, fd_dust, &
203 & sigma_nuc, rhi_cell
204 integer, intent(in) :: ccn_param, in_param, nccn
205
206 real(r8), dimension(:), intent(inout) :: ccn_diagr8
207
208 real(r8), intent(out) :: cdncr8, smaxliqr8, incr8, smaxicer8, &
209 & nheticer8, INimmr8, dINimmr8, Ncdepr8, Ncdhfr8, sc_icer8, &
210 & fdust_immr8, fdust_depr8, fdust_dhfr8, nlimr8, pfrz_inc_r8
211
212 type(aerprops) :: aeraux
213
214
215 integer :: k, n, i, j, naux
216
217
218 real*8 :: nact, wparc, tparc,pparc, accom,sigw, smax, antot, &
219 & ccn_at_s, sigwparc
220! real*8, allocatable, dimension(:) :: smax_diag
221 real*8, dimension(nccn) :: smax_diag
222
223
224 real*8 :: nhet, nice, smaxice, nlim, air_den, frac, norg, nbc, &
225 & nhom, dorg, dbc, kappa, inimm, dinimm, aux
226
227!Anning Cheng move allocable array here for thread safety, 5/4/2016
228! real*8, allocatable, dimension(:) :: sg_par, tp_par, dpg_par,
229 real*8, dimension(nsmx_par) :: sg_par, tp_par, dpg_par
230 &, sig_par, vhf_par, ams_par
231 &, dens_par, deni_par, amfs_par
232 &, kappa_par, ndust_ice, sigdust_ice
233 &, ddust_ice
234! real*8, allocatable, dimension(:) :: ndust_ice, sigdust_ice,
235! & ddust_ice
236 real*8 :: temp_par, pres_par
237 real*8 :: akoh_par, alfa_par, bet2_par
238 real*8 aka_par, dv_par, psat_par, dair_par,surt_par,ddry_ice,
239 & np_ice,nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
240 & g1_ice, g2_ice,gdoin_ice, z_ice, norg_ice, sigorg_ice,
241 & dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
242 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
243 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
244 & si0bc_ice,nbc_ice,sc_ice,nhet_dep,fdust_imm,fdust_dep,fdust_dhf,
245 & waux_ice,fdrop_dust,fdrop_bc,d_preex, n_preex,
246 & one_over_tao_preex,p_ice, t_ice,miuv_ice,nhet_dhf,denice_ice,
247 & vpresw_ice,vpresi_ice,denair_ice
248 real*8 :: ntot
249 integer :: nmodes,act_param,nbindust_ice
250 logical :: use_av_v
251
252!=============inputs================
253 tparc = tparc_in
254 pparc = pparc_in
255 sigwparc = sigwparc_in
256 miuv_ice = wparc_ls
257 air_den = pparc*28.8d-3/rgas_par/tparc
258 n_preex = max(npre_in*air_den, zero_par)
259 d_preex = max(dpre_in, 1.0e-9)
260 use_av_v = use_average_v
261 act_param = 2
262 typeofspec_ice = 5
263
264
265
266 smaxicer8 = zero_par
267 smaxice = zero_par
268 cdncr8 = zero_par
269 smaxliqr8 = zero_par
270 incr8 = zero_par
271 smaxice = max(2.349d0-(tparc/259d0) -1.0 , 0.0)
272 nheticer8 = zero_par
273 nlimr8 = zero_par
274 sc_ice = max(2.349d0-(tparc/259d0), 1.0)
275 If (tparc > thom) sc_ice =1.0
276
277 inimmr8 = zero_par
278 dinimmr8 = zero_par
279 ncdepr8 = zero_par
280 ncdhfr8 = zero_par
281 fdust_immr8 = zero_par
282 fdust_dhfr8 = zero_par
283 fdust_depr8 = zero_par
284 fdust_imm = zero_par
285 fdust_dhf = zero_par
286 fdust_dep = zero_par
287 pfrz_inc_r8 = zero_par
288
289 nact = zero_par
290 smax = zero_par
291 sc_icer8 = sc_ice
292
293 is_gocart = .true.
294
295 if (sum(aer_props%num) <= 1.0e2) then
296 return
297 end if
298
299 nmodes = max(aer_props%nmods, 1)
300
301! allocate (dpg_par(nmodes))
302! allocate (vhf_par(nmodes))
303! allocate (ams_par(nmodes))
304! allocate (dens_par(nmodes))
305! allocate (sig_par(nmodes))
306! allocate (tp_par(nmodes))
307! allocate (amfs_par(nmodes))
308! allocate (deni_par(nmodes))
309! allocate (sg_par(nmodes))
310! allocate (smax_diag(size(ccn_diagr8)))
311
312! allocate (kappa_par(nmodes))
313
314
315 smax_diag = 0.01
316 sigw = zero_par
317 do n=1,nmodes
318 dpg_par(n) = zero_par
319 vhf_par(n) = zero_par
320 ams_par(n) = zero_par
321 dens_par(n) = zero_par
322 sig_par(n) = 1d0
323 tp_par(n) = zero_par
324 amfs_par(n) = zero_par
325 deni_par(n) = zero_par
326 kappa_par(n) = zero_par
327 enddo
328
329 call init_aer(aeraux)
330
331! if (lprnt) write(0,*)' in aero Aer_Props%num='
332! &,Aer_Props%num,' nmodes=',nmodes,' air_den=',air_den
333! if (lprnt) write(0,*)' in aero Aer_Props%kap='
334! &,Aer_Props%kap
335
336 antot = 0.0
337
338 do n=1,nmodes
339! tp_par(n) = DBLE(Aer_Props%num(n))*air_den
340! dpg_par(n) = max(DBLE(Aer_Props%dpg(n)), 1.0e-10)
341! sig_par(n) = DBLE(Aer_Props%sig(n))
342! kappa_par(n) = max(DBLE(Aer_Props%kap(n)), 0.001)
343! dens_par(n) = DBLE(Aer_Props%den(n))
344
345 tp_par(n) = aer_props%num(n) * air_den
346 dpg_par(n) = max(aer_props%dpg(n), 1.0e-10)
347 sig_par(n) = aer_props%sig(n)
348 kappa_par(n) = max(aer_props%kap(n), 0.001)
349 dens_par(n) = aer_props%den(n)
350 vhf_par(n) = 3.0
351 if (kappa_par(n) > 0.01) then
352 ams_par(n) = 18.0e-3*1.7*3.0/kappa_par(n)
353 else
354 ams_par(n) = 900.0e-3
355 tp_par(n) = 0.0
356 endif
357 amfs_par(n) = 1.0
358 deni_par(n) = dens_par(n)
359 antot = antot + tp_par(n)
360
361! if (lprnt) write(0,*)' n=',n,' tp_par=',tp_par(n),' antot=',antot
362! &,' Aer_Props%num=',Aer_Props%num(n),' kappa_par=',kappa_par(n)
363! &,' air_den=',air_den
364 enddo
365
366! kappa_par = max(kappa_par, 0.001)
367! dpg_par = max(dpg_par, 1.0e-10)
368 temp_par = max(tparc, 245.0)
369 pres_par = max(pparc, 34000.0)
370
371! antot = sum(tp_par)
372 ntot = antot
373
374 wparc = max(max(0.8d0*sigwparc, 0.01)+ wparc_ls, 0.01)
375
376
377 act_param = ccn_param
378
379!============== Calculate cloud droplet number concentration===================
380
381! if (lprnt) write(0,*)' in aero tparc=',tparc,' antot=',antot
382! if (lprnt) write(0,*)' in aero tp_par=',tp_par(1:nmodes)
383
384 if (tparc > 245.0) then
385 if (antot > 1.0) then
386
387 call ccnspec (tparc,pparc,nmodes
388 &, amfs_par,dens_par,deni_par,vhf_par,ams_par
389 &, sg_par,tp_par,akoh_par,surt_par,temp_par,pres_par
390 &, dv_par,act_param,aka_par, psat_par,dair_par
391 &, ntot,dpg_par)
392
393 if (wparc >= 0.005) then
394 if (act_param > 1) then
395
396 call arg_activ (wparc,0.d0,nact,smax,nmodes,tp_par
397 &, dpg_par,kappa_par,sig_par,temp_par, pres_par)
398
399 else
400
401 call pdfactiv (wparc,0.d0,nact,smax,nmodes
402 &, alfa_par,bet2_par,akoh_par,sg_par,tp_par
403 &, temp_par, pres_par,aka_par, dv_par, psat_par
404 &, dair_par,ntot,sig_par)
405
406 endif
407 endif
408
409 cdncr8 = max(nact/air_den, zero_par)
410 smaxliqr8 = max(smax, zero_par)
411
412! if (lprnt) write(0,*)' in aero cdncr8=',cdncr8,' nact=',nact,
413! &' air_den=',air_den,' wparc=',wparc,' act_param=',act_param
414!============ Calculate diagnostic CCN number concentration==================
415
416! smax_diag = ccn_diagr8
417
418! do k =1, size (smax_diag)
419! do k =1, nccn
420! call ccn_at_super (smax_diag(k), ccn_at_s,nmodes,
421! & sig_par,sg_par,tp_par)
422! ccn_diagr8 (k) = ccn_at_s
423! end do
424
425 end if
426 end if
427
428
429! ==========================================================================================
430! ==========================================================================================
431!========================== Ice crystal nucleation parameterization ======================
432! ==========================================================================================
433 dbc_ice = 1.0e-9
434 nbc_ice = zero_par
435 norg_ice = zero_par
436 dorg_ice = 1.0e-9
437 sigbc_ice = zero_par
438 sigorg_ice = zero_par
439 ddry_ice = 1.0e-9
440 np_ice = zero_par
441 nin_ice = 0.
442 kdust_ice = 0.
443 kbc_ice = 0.
444 shdust_ice = 0.
445 shbc_ice = 0.
446 effdust_ice = 0.
447 effbc_ice = 0.
448 del1dust_ice = 0.
449 si0dust_ice = 0.
450 del1bc_ice = 0.
451 si0bc_ice = 0.
452
453 naux = 0
454 do k = 1, nmodes
455 if (kappa_par(k) > 0.1) then
456 np_ice = np_ice + tp_par(k)
457 ddry_ice = ddry_ice + dpg_par(k)
458 naux = naux + 1
459 end if
460 end do
461 ddry_ice = ddry_ice/max(naux , 1)
462 frac = 1.0
463 np_ice = frac*np_ice
464
465!get dust from input structure
466
467 call getinsubset(1, aer_props, aeraux)
468 nbindust_ice = max(aeraux%nmods, 1)
469
470! allocate(ndust_ice(nbindust_ice))
471! allocate(sigdust_ice(nbindust_ice))
472! allocate(ddust_ice(nbindust_ice))
473
474 ddust_ice = zero_par
475 ndust_ice = zero_par
476 sigdust_ice = zero_par
477 do n=1,nbindust_ice
478 ddust_ice(n) = dble(aeraux%dpg(n))
479 ndust_ice(n) = dble(aeraux%num(n))*air_den
480 sigdust_ice(n) = dble(aeraux%sig(n))
481 enddo
482
483
484!Black carbon. Only a single mode considered. Use average size and sigma
485
486 call getinsubset(2, aer_props, aeraux)
487 naux = max(aeraux%nmods, 1)
488 dbc_ice = dble(sum(aeraux%dpg(1:naux)))/naux
489 nbc_ice = dble(sum(aeraux%num(1:naux)))*air_den
490 sigbc_ice = dble(sum(aeraux%sig(1:naux)))/naux
491
492
493
494 call getinsubset(3, aer_props, aeraux)
495 naux = max(aeraux%nmods, 1)
496 dorg_ice = dble(sum(aeraux%dpg(1:naux)))/naux
497 norg_ice = dble(sum(aeraux%num(1:naux)))*air_den
498 sigorg_ice = dble(sum(aeraux%sig(1:naux)))/naux
499
500
501 nhet = zero_par
502 nice = zero_par
503 nlim = zero_par
504 inimm = zero_par
505 dinimm = zero_par
506 nhet_dep = zero_par
507 nhet_dhf = zero_par
508 antot=sum(ndust_ice)+ norg_ice+ nbc_ice+ np_ice
509
510 sigwparc=max(0.01, sigwparc)
511 sigwparc=min(5.0, sigwparc)
512 waux_ice=max(wparc_ls + sigwparc*0.8, 0.01)
513
514
515!===========Calculate nucleated crystal number. Follows Barahona and Nenes (2008, 2009)==============
516
517 purehet_ice= .false.
518 purehom_ice= .false.
519
520
521 if (antot > 1.0e2) then
522 if (tparc < to_ice) then
523
524 CALL prop_ice(tparc, pparc, denice_ice,ddry_ice
525 &, nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice
526 &, g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice
527 &, kdust_ice, kbc_ice, shdust_ice, shbc_ice
528 &, effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice
529 &, si0bc_ice,nbc_ice,d_preex, n_preex
530 &, one_over_tao_preex,p_ice, t_ice,act_param
531 &, ndust_ice,vpresw_ice,vpresi_ice,denair_ice)
532
533
534 if (tparc > thom) then
535
536 fdrop_dust = fd_dust
537 fdrop_bc = fd_soot
538
539 if (sum(ndust_ice)+ norg_ice+ nbc_ice .gt. 1.e3) then
540
541 call inimmersion(inimm, dinimm, waux_ice,dbc_ice,sigbc_ice
542 &, nbc_ice,fdust_imm,fdrop_dust,fdrop_bc,ndust_ice
543 &, sigdust_ice,ddust_ice,nbindust_ice
544 &, vpresw_ice,vpresi_ice,t_ice)
545
546
547 ndust_ice = max(ndust_ice*(1.0-fdrop_dust), 0.0)
548 nbc_ice = max(nbc_ice*(1.0-fdrop_bc), 0.0)
549
550
551 call iceparam (sigwparc, denice_ice,ddry_ice,np_ice
552 &, nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice
553 &, g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,sc_ice
554 &, norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice
555 &, kdust_ice, kbc_ice, shdust_ice, shbc_ice
556 &, effdust_ice, effbc_ice, del1dust_ice, si0dust_ice
557 &, del1bc_ice, si0bc_ice,nbc_ice,one_over_tao_preex
558 &, nhet, nice, smaxice, nlim,nhet_dep,nhet_dhf,fdust_dep
559 &, p_ice,t_ice,ndust_ice,sigdust_ice,ddust_ice,nbindust_ice
560 &, use_av_v,miuv_ice,vpresw_ice,vpresi_ice,denair_ice)
561 end if
562
563 sc_ice = 1.0
564
565 else
566
567 call iceparam (sigwparc, denice_ice,ddry_ice,np_ice
568 &, nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice
569 &, g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,sc_ice
570 &, norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice
571 &, kdust_ice, kbc_ice, shdust_ice, shbc_ice
572 &, effdust_ice, effbc_ice, del1dust_ice, si0dust_ice
573 &, del1bc_ice, si0bc_ice,nbc_ice,one_over_tao_preex
574 &, nhet, nice, smaxice, nlim,nhet_dep,nhet_dhf,fdust_dep
575 &, p_ice, t_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice
576 &, use_av_v,miuv_ice,vpresw_ice,vpresi_ice,denair_ice)
577
578 end if
579
580 aux = (sc_ice - rhi_cell)/(sq2_par*sigma_nuc)
581
582 pfrz_inc_r8 = 1.0d0- 0.5d0*(1.0d0+erf(aux))
583 pfrz_inc_r8 = min(max(pfrz_inc_r8, 0.0), 0.999)
584
585 end if
586 end if
587
588
589!======================== use sc_ice only for cirrus
590 If (tparc > thom) sc_ice =1.0
591!==========================
592
593!All # m-3 except those passed to MG later
594 smaxicer8 = min(max(smaxice, zero_par), 2.0)
595 nheticer8 = min(max(nhet, zero_par), 1e10)
596 incr8 = min(max(nice/air_den, zero_par), 1e10)
597 nlimr8 = min(max(nlim, zero_par), 1e10)
598 sc_icer8 = min(max(sc_ice, 1.0), 2.0)
599 inimmr8 = min(max(inimm, zero_par), 1e10)
600 dinimmr8 = min(max(dinimm/air_den, zero_par), 1e10)
601 ncdepr8 = min(max(nhet_dep, zero_par), 1e10)
602 ncdhfr8 = min(max(nhet_dhf, zero_par), 1e10)
603 fdust_immr8 = min(max(fdust_imm, zero_par), 1e10)
604 fdust_depr8 = min(max(fdust_dep, zero_par), 1e10)
605 fdust_dhfr8 = min(max(fdust_dhf, zero_par), 1e10)
606
607! deallocate (ndust_ice)
608! deallocate (sigdust_ice)
609! deallocate (ddust_ice)
610! deallocate (dpg_par)
611! deallocate (vhf_par)
612! deallocate (ams_par)
613! deallocate (dens_par)
614! deallocate (sig_par)
615! deallocate (tp_par)
616! deallocate (amfs_par)
617! deallocate (deni_par)
618! deallocate (sg_par)
619! deallocate (smax_diag)
620
621! deallocate (kappa_par)
622
623
624
6252033 return
626
627 END subroutine aerosol_activate
628!
629
630
631!=======================================================================
632!
636 SUBROUTINE aerconversion_base ()
637
638 integer, parameter :: NMDM = 20
639 real, dimension(NMDM) :: TPI, DPGI, SIGI, DENSI, KAPPAS, FDUST,
640 & fsoot, forg, tpi_aux, dpgi_aux, sigi_aux
641
642 real:: aux
643 integer:: nmod, K
644
645 do k=1,nmdm
646 tpi(k) = 0.0
647 dpgi(k) = 1.0e-9
648 sigi(k) = 2.0
649 densi(k) = 2200.0
650 kappas(k) = 0.01
651 fdust(k) = 0.0
652 fsoot(k) = 0.0
653 forg(k) = 0.0
654 enddo
655 nmod = 13
656
657! Gocart aerosol size distributions for dust
658
659!!!!!!!!!!!!!!======================================
660!!!!!!!!! Dust
661!!!!!!!!!!!!!!======================================
662
663
664 do k=1,5
665 sigi(k) = log(1.8)
666 densi(k) = 1700.0
667 kappas(k) = 0.0001
668 fdust(k) = 1.0
669 enddo
670
671
672
673 dpgi(1) = 1.46e-6
674
675 dpgi(2) = 2.80e-6
676
677 dpgi(3) = 4.80e-6
678!! Dust 4: 3-6
679 dpgi(4) = 9.0e-6
680!! Dust 5: 6-10
681 dpgi(5) = 16.0e-6
682
683 DO k =1 , 5
684 tpi(k) = 6.0/(densi(k)*pi_par*exp(4.5*sigi(k)*sigi(k))*dpgi(k)*
685 & dpgi(k)*dpgi(k))
686 END DO
687
688!!!!!!!!!!!!!!======================================
689!!!!!!!!! Sea Salt (Using 3 modes based on Barahona et al. GMD. 2014.
690!!!!!!!!!!!!!!======================================
691
692
693
694 densi(6:8) = 2200.0
695 kappas(6:8) = 1.28
696
697
698
699 tpi(6) = 100e6
700 dpgi(6) = .01e-6
701 sigi(6) = log(1.6)
702
703 tpi(7) = 60.0e6
704 dpgi(7) = 0.071e-6
705 sigi(7) = log(2.0)
706
707 tpi(8) = 3.1e6
708 dpgi(8) = 0.62e-6
709 sigi(8) = log(2.7)
710
711 aux = 0.
712 DO k =6 , 8
713 aux = (tpi(k)*densi(k)*pi_par*exp(4.5*sigi(k)*sigi(k))*dpgi(k)*
714 & dpgi(k)*dpgi(k))/6.0 + aux
715 END DO
716 base_mass_ss = aux
717
718
719
720
721 kappas(9:11) = 0.65
722 densi(9:11) = 1650.0
723
724! Different size distributions for polluted and clean environments
725
726
727 tpi(9) = 1.06e11
728 dpgi(9) = .014e-6
729 sigi(9) = log(1.8d0)
730
731 tpi(10) = 3.2e10
732 dpgi(10) = 0.054e-6
733 sigi(10) = log(2.16)
734
735 tpi(11) = 5.4e6
736 dpgi(11) = 0.86e-6
737 sigi(11) = log(2.21)
738
739 aux = 0.
740 DO k =9, 11
741 aux = (tpi(k)*densi(k)*pi_par*exp(4.5*sigi(k)*sigi(k))*dpgi(k)*
742 & dpgi(k)*dpgi(k))/6.0 + aux
743 END DO
744 base_mass_so4_polluted = aux
745
746
747
748
749 tpi_aux(9) = 1.0e9
750 dpgi_aux(9) = .016e-6
751 sigi_aux(9) = log(1.6d0)
752
753 tpi_aux(10) = 8.0e8
754 dpgi_aux(10) = 0.067e-6
755 sigi_aux(10) = log(2.1)
756
757 tpi_aux(11) = 2.0e6
758 dpgi_aux(11) = 0.93e-6
759 sigi_aux(11) = log(2.2)
760
761
762 aux = 0.
763 DO k =9, 11
764 aux =(tpi_aux(k)*densi(k)*pi_par*exp(4.5*sigi_aux(k)*
765 & sigi_aux(k))*dpgi_aux(k)*dpgi_aux(k)*dpgi_aux(k))/6.0 + aux
766 END DO
767 base_mass_so4_clean = aux
768
769
770
771
772 dpgi(12) = 0.0118*2.e-6
773 sigi(12) = log(2.00)
774 densi(12) = 1600.0
775 kappas(12) = 0.0001
776 fsoot(12) = 1.0
777 tpi(12) = 6.0/(densi(12)*pi_par*exp(4.5*sigi(12)*sigi(12))*
778 & dpgi(12)*dpgi(12)*dpgi(12))
779
780
781
782 dpgi(13) = 0.0212*2.e-6
783 sigi(13) = log(2.20)
784 densi(13) = 900.0
785 kappas(13) = 0.0001
786 forg(13) = 1.0
787 tpi(13) = 6.0/(densi(13)*pi_par*exp(4.5*sigi(13)*sigi(13))*
788 & dpgi(13)*dpgi(13)*dpgi(13))
789
790 call init_aer(aerpr_base_polluted)
791 call init_aer(aerpr_base_clean)
792
793
794 do k=1,nmod
795 aerpr_base_polluted%num(k) = tpi(k)
796 aerpr_base_polluted%dpg(k) = dpgi(k)
797 aerpr_base_polluted%sig(k) = sigi(k)
798 aerpr_base_polluted%kap(k) = kappas(k)
799 aerpr_base_polluted%den(k) = densi(k)
800 aerpr_base_polluted%fdust(k) = fdust(k)
801 aerpr_base_polluted%fsoot(k) = fsoot(k)
802 aerpr_base_polluted%forg(k) = forg(k)
803 enddo
804 aerpr_base_polluted%nmods = nmod
805
806 aerpr_base_clean = aerpr_base_polluted
807 aerpr_base_clean%num(9:11) = tpi_aux(9:11)
808 aerpr_base_clean%dpg(9:11) = dpgi_aux(9:11)
809 aerpr_base_clean%sig(9:11) = sigi_aux(9:11)
810
811 RETURN
812!
813 END SUBROUTINE aerconversion_base
814
815
816!=======================================================================
817!
821 SUBROUTINE aerconversion (aer_mass, AerPr, kappa, SULFATE, ORG, &
822 & BCARBON, DUST, SEASALT)
823
824
825 type(aerprops), intent (out), dimension(:,:,:) :: aerpr
826
827
828 real, dimension(:,:,:,:), intent(in) :: aer_mass
829 real, intent (out), dimension(:,:,:) :: kappa, dust, sulfate, &
830 & BCARBON, ORG, SEASALT
831 real:: aux, densso4, densorg, k_so4, k_org, k_ss, tot_mass, dens,
832 & kappa_aux, tx1
833
834 integer :: i,j,k,l
835 integer :: im, jm, lm
836 type(aerprops) :: aeroaux
837 real, dimension(size(aer_mass,4)) :: aer_mass_tmp
838
839 im = size(aer_mass,1)
840 jm = size(aer_mass,2)
841 lm = size(aer_mass,3)
842
843 call init_aer(aeroaux)
844
845 do k = 1, lm
846 do j = 1, jm
847 do i = 1, im
848 aer_mass_tmp(:) = aer_mass(i,j,k,:)
849
850
851 tot_mass = aer_mass_tmp(11) + aer_mass_tmp(15)
852 & + aer_mass_tmp(14)
853 densso4 = 1700.0
854 densorg = 1600.0
855 k_so4 = 0.65
856 k_org = 0.2
857 kappa_aux = 0.65
858
859
860
861 if (tot_mass > 2.0e-20) then
862 tx1 = 1.0 / tot_mass
863 dens = (aer_mass_tmp(11)*densso4
864 & + sum(aer_mass_tmp(14:15))*densorg) * tx1
865 kappa_aux = (aer_mass_tmp(11)*k_so4
866 & + sum(aer_mass_tmp(14:15))*k_org) * tx1
867 else
868 dens = 1750.0
869 kappa_aux = 0.65
870 end if
871
872
873 if (tot_mass > 5.0e-8) then
874 aeroaux = aerpr_base_polluted
875 aeroaux%num(9:11) = aeroaux%num(9:11)*tot_mass
876 & / base_mass_so4_polluted
877 else
878 aeroaux = aerpr_base_clean
879 aeroaux%num(9:11) = aeroaux%num(9:11)*tot_mass
880 & / base_mass_so4_clean
881 end if
882
883 aeroaux%kap(9:11) = max(kappa_aux, 0.001)
884 aeroaux%den(9:11) = dens
885 sulfate(i,j,k) = sum(aeroaux%num(9:11))
886 kappa_aux = kappa_aux*tot_mass
887
888
889
890
891 aeroaux%num(1:5) = aeroaux%num(1:5) *aer_mass_tmp(1:5)
892 kappa_aux = kappa_aux
893 & + aeroaux%kap(1)*sum(aer_mass_tmp(1:5))
894 dust(i,j,k) = sum(aeroaux%num(1:5))
895
896
897 tot_mass = sum(aer_mass_tmp(6:10))
898 aeroaux%num(6:8) = aeroaux%num(6:8)*tot_mass/base_mass_ss
899 kappa_aux = kappa_aux + aeroaux%kap(6)*tot_mass
900 seasalt(i,j,k ) = sum(aeroaux%num(6:8))
901
902
903 aeroaux%num(12) = aeroaux%num(12) *aer_mass_tmp(13)
904 kappa_aux = kappa_aux
905 & + aeroaux%kap(12)*aer_mass_tmp(13)
906 bcarbon(i,j,k) = aeroaux%num(12)
907
908
909 aeroaux%num(13) = aeroaux%num(13) *aer_mass_tmp(15)
910 org(i, j, k) = aeroaux%num(13)
911 tot_mass = sum(aer_mass_tmp)
912
913 if (tot_mass > 0.0) then
914 kappa(i,j,k) = kappa_aux/tot_mass
915 end if
916
917 aerpr(i,j,k) = aeroaux
918
919 end do
920 end do
921 end do
922
923 RETURN
924!
925 END SUBROUTINE aerconversion
926
929 SUBROUTINE aerconversion1 (aer_mass, AerPr)
930
931
932 type(aerprops), dimension(:,:) :: aerpr
933! type(AerProps), intent (out), dimension(:,:) :: AerPr
934
935
936 real, dimension(:,:,:) :: aer_mass
937! real, dimension(:,:,:), intent(in) :: aer_mass
938 real:: aux, densso4, densorg, k_so4, k_org, k_ss, tot_mass, dens,
939 & kappa_aux, tx1, tx2
940
941 integer :: i,k,l,im,lm
942 type(aerprops) :: aeroaux
943 real, dimension(size(aer_mass,3)) :: aer_mass_tmp
944
945 im = size(aer_mass,1)
946 lm = size(aer_mass,2)
947
948 call init_aer(aeroaux)
949
950 do k = 1, lm
951 do i = 1, im
952 aer_mass_tmp(:) = aer_mass(i,k,:)
953
954
955 tot_mass = aer_mass_tmp(11) + aer_mass_tmp(15)
956 & + aer_mass_tmp(14)
957 densso4 = 1700.0
958 densorg = 1600.0
959 k_so4 = 0.65
960 k_org = 0.2
961 kappa_aux = 0.65
962
963
964
965 if (tot_mass > 2.0e-20) then
966 tx1 = 1.0 / tot_mass
967 tx2 = aer_mass_tmp(14) + aer_mass_tmp(15)
968 dens = (aer_mass_tmp(11)*densso4 + tx2*densorg) * tx1
969 kappa_aux = (aer_mass_tmp(11)*k_so4 + tx2*k_org) * tx1
970 else
971 dens = 1750.0
972 kappa_aux = 0.65
973 end if
974
975
976 if (tot_mass > 5.0e-8) then
977 aeroaux = aerpr_base_polluted
978 aeroaux%num(9:11) = aeroaux%num(9:11)*tot_mass
979 & / base_mass_so4_polluted
980 else
981 aeroaux = aerpr_base_clean
982 aeroaux%num(9:11) = aeroaux%num(9:11)*tot_mass
983 & / base_mass_so4_clean
984 end if
985
986 aeroaux%kap(9:11) = max(kappa_aux, 0.001)
987 aeroaux%den(9:11) = dens
988 kappa_aux = kappa_aux*tot_mass
989
990
991
992
993 aeroaux%num(1:5) = aeroaux%num(1:5) *aer_mass_tmp(1:5)
994 kappa_aux = kappa_aux
995 & + aeroaux%kap(1)*sum(aer_mass_tmp(1:5))
996
997
998 tot_mass = sum(aer_mass_tmp(6:10))
999 aeroaux%num(6:8) = aeroaux%num(6:8)*tot_mass/base_mass_ss
1000 kappa_aux = kappa_aux + aeroaux%kap(6)*tot_mass
1001
1002
1003 aeroaux%num(12) = aeroaux%num(12)*aer_mass_tmp(13)
1004 kappa_aux = kappa_aux+aeroaux%kap(12)*aer_mass_tmp(13)
1005
1006
1007 aeroaux%num(13) = aeroaux%num(13) *aer_mass_tmp(15)
1008 tot_mass = sum(aer_mass_tmp)
1009
1010
1011 aerpr(i,k) = aeroaux
1012
1013 end do
1014 end do
1015
1016 RETURN
1017!
1018 END SUBROUTINE aerconversion1
1019
1020!=======================================================================
1021!=======================================================================
1022!=======================================================================
1023
1026 subroutine vertical_vel_variance(omeg, lc_turb, tm_gw, pm_gw,
1027 & rad_cool, uwind_gw, tausurf_gw, nm_gw, LCCIRRUS, Nct, Wct, ksa1,
1028 & fcn, KH, FRLAND, ZPBL, Z, maxkhpbl, wparc_ls, wparc_gw,
1029 & wparc_cgw, wparc_turb, LTS)
1030
1031
1032 real(r8), intent(in) :: omeg, tm_gw, lc_turb, rad_cool, uwind_gw,
1033 & pm_gw
1034 real , intent(in) :: lccirrus, kh, zpbl, z, frland, nm_gw,
1035 & tausurf_gw, ksa1, fcn, maxkhpbl, nct, wct, lts
1036
1037 real(r8), intent(out) :: wparc_ls, wparc_gw, wparc_cgw,
1038 & wparc_turb
1039
1040 real(r8) :: rho_gw, k_gw, h_gw, c2_gw, dummyw, maxkh, wbreak
1041
1042
1043
1044!!!:========= mean V Large scale and radiative cooling
1045 rho_gw = pm_gw*28.8d-3/rgas_par/tm_gw
1046
1047
1048 wparc_ls =-omeg/rho_gw/grav_ice + cpa_ice*rad_cool/grav_ice
1049
1050!!!======== Orographic Gravity gwave (and brackground) initiated (According to Barahona et al. 2013 GMD)
1051
1052 wparc_gw = 0.0
1053 k_gw = 2d0*pi_par/lccirrus
1054
1055 h_gw= k_gw*rho_gw*uwind_gw*nm_gw
1056
1057 if (h_gw .gt. 0.0) then
1058 h_gw=sqrt(2.0*tausurf_gw/h_gw)
1059 else
1060 h_gw = 0.0
1061 end if
1062
1063 wbreak = 0.133*k_gw*uwind_gw/nm_gw
1064
1065 wparc_gw = k_gw*uwind_gw*h_gw*0.133
1066 wparc_gw = min(wparc_gw, wbreak)
1067 wparc_gw = wparc_gw*wparc_gw
1068
1069
1070!!!======== Subgrid variability from Convective Sources According to Barahona et al. 2014 in prep
1071
1072 wparc_cgw = 0.0
1073 c2_gw = (nm_gw+nct)/nct
1074 wparc_cgw = sqrt(ksa1)*fcn*c2_gw*wct*0.6334
1075 wparc_cgw = min(wparc_cgw, wbreak)
1076 wparc_cgw = wparc_cgw*wparc_cgw
1077
1078!!!:=========Subgrid scale variance from turbulence
1079
1080
1081 wparc_turb=kh/lc_turb
1082
1083
1084 if (lts > 20.0) then
1085 wparc_turb=maxkhpbl/lc_turb
1086 end if
1087
1088
1089 wparc_turb= wparc_turb*wparc_turb
1090
1091
1092
1093 end subroutine vertical_vel_variance
1094!=======================================================================
1095!=======================================================================
1096!=======================================================================
1097
1098!=======================================================================
1101 subroutine getinsubset(typ, aerin, aerout)
1102
1103! typ : type of aerosol needed: 1 dust, 2 soot, 3 organics
1104! nbins: number of modes with num>0
1105
1106 type(aerprops), intent (in) :: aerin
1107 type(aerprops), intent (inout) :: aerout
1108 integer, intent(in) :: typ
1109
1110 integer:: k, bin
1111
1112 call init_aer(aerout)
1113 bin = 0
1114
1115 do k=1, aerin%nmods
1116
1117 if (typ == 1) then
1118 if (aerin%fdust(k) > 0.9) then
1119 bin = bin + 1
1120 call copy_mode(aerout,aerin, k,bin)
1121 end if
1122 elseif (typ == 2) then
1123 if (aerin%fsoot(k) > 0.9) then
1124 bin = bin + 1
1125 call copy_mode(aerout,aerin, k,bin)
1126 end if
1127 elseif (typ == 3) then
1128 if (aerin%forg(k) > 0.9) then
1129 bin = bin + 1
1130 call copy_mode(aerout,aerin, k,bin)
1131 end if
1132 end if
1133
1134 end do
1135
1136 aerout%nmods = max(bin, 1)
1137
1138 end subroutine getinsubset
1139
1140!========================subroutines to handle aer strucuture=====================================
1141
1144 subroutine copy_aer(a,b)
1145
1146 type (AerProps), intent(out) :: a
1147 type (AerProps), intent(in) :: b
1148
1149 a%num = b%num
1150 a%sig = b%sig
1151 a%dpg = b%dpg
1152 a%kap = b%kap
1153 a%den = b%den
1154 a%fdust = b%fdust
1155 a%fsoot = b%fsoot
1156 a%forg = b%forg
1157 a%nmods = b%nmods
1158
1159 end subroutine copy_aer
1160
1163 subroutine copy_mode(a_out,a_in, mode_in, mode_out)
1164 type (AerProps), intent(out) :: a_out
1165 type (AerProps), intent(in) :: a_in
1166 integer, intent (in) :: mode_in, mode_out
1167
1168 a_out%num(mode_out)= a_in%num(mode_in)
1169 a_out%sig(mode_out) = a_in%sig(mode_in)
1170 a_out%dpg(mode_out) = a_in%dpg(mode_in)
1171 a_out%kap(mode_out) = a_in%kap(mode_in)
1172 a_out%den(mode_out) = a_in%den(mode_in)
1173 a_out%fdust(mode_out) = a_in%fdust(mode_in)
1174 a_out%fsoot(mode_out) = a_in%fsoot(mode_in)
1175 a_out%forg(mode_out) = a_in%forg(mode_in)
1176
1177 end subroutine copy_mode
1178
1181 subroutine init_aer(aerout)
1182
1183 type (aerprops), intent(inout) :: aerout
1184 integer n
1185
1186 do n=1,nsmx_par
1187 aerout%num(n) = 0.
1188 aerout%dpg(n) = 1.0e-9
1189 aerout%sig(n) = 2.0
1190 aerout%kap(n) = 0.2
1191 aerout%den(n) = 2200.0
1192 aerout%fdust(n) = 0.0
1193 aerout%fsoot(n) = 0.0
1194 aerout%forg(n) = 0.0
1195 enddo
1196 aerout%nmods = 1
1197
1198 end subroutine init_aer
1199
1200
1201
1207 subroutine arg_activ (wparc,sigw,nact,smax,nmodes,tp_par, &
1208 & dpg_par,kappa_par,sig_par,temp_par, pres_par)
1209
1210
1211 real*8, intent(in) :: wparc, sigw
1212 integer, intent(in):: nmodes
1213 real*8, intent(out) :: nact, smax
1214 real*8 :: smi(nmodes),tp_par(nmodes)
1215 real*8, dimension(nmodes)::dpg_par,kappa_par,sig_par
1216
1217 real*8 :: alfa, beta, akoh, g, t, fi, gi, nui, citai, ui, aux1,
1218 & pact, auxx, aux,temp_par, pres_par
1219 integer :: I, INDEX
1220
1221 t = min(max(temp_par, 243.0), 323.0)
1222 alfa=2.8915e-08*(t*t) - 2.1328e-05*t + 4.2523e-03
1223 beta=exp(3.49996e-04*t*t - 2.27938e-01*t + 4.20901e+01)
1224 g=exp(-2.94362e-06*t*t*t + 2.77941e-03*t*t - 8.92889e-01*t +
1225 & 1.18787e+02)
1226 akoh= 0.66e-6/t
1227
1228
1229 DO i = 1, nmodes
1230 aux =0.667*akoh/dpg_par(i)
1231 smi(i) = (aux*sqrt(aux))/sqrt(2.0*kappa_par(i))
1232 END DO
1233
1234
1235
1236 smi=max(smi, 1.0e-5)
1237 aux =alfa*wparc*g
1238 aux = aux*sqrt(aux)/(2.d0*pi_par*980.d0*beta)
1239 citai = 0.667*akoh*sqrt(alfa*wparc*g)
1240
1241 auxx=0.0
1242 DO index =1, nmodes
1243 if (tp_par(index) .gt. 0.0) then
1244 fi=0.5*exp(2.5*sig_par(index))
1245 gi=1.0+0.25*sig_par(index)
1246 nui=aux/tp_par(index)
1247 aux1=fi*((citai/nui)*sqrt(citai/nui)) + gi*(smi(index)*
1248 &smi(index) /(nui+(3.0*citai)))**0.75
1249 aux1=aux1/(smi(index)*smi(index))
1250 auxx=auxx+aux1
1251
1252 end if
1253 end do
1254
1255
1256 if (auxx .gt. 0.0) then
1257 smax = 1/sqrt(auxx)
1258 else
1259 nact = 0.0
1260 return
1261 end if
1262
1263
1264 auxx = 0.0
1265
1266
1267 DO index = 1, nmodes
1268 ui = sq2_par*log(smi(index)/smax)/3.0
1269 aux1 = 0.5*tp_par(index)*(1.0-erfapp(ui))
1270 auxx = auxx + aux1
1271 END DO
1272
1273 nact = auxx
1274
1275
1276 End Subroutine arg_activ
1277
1278!===================================================================================================
1279
1280!===================================================================================================
1281
1282
1283
1284
1285!=====================================================================================================
1287 subroutine ccn_at_super (super,ccn_at_s,nmodes, &
1288 & sig_par,sg_par,tp_par)
1289
1290 integer :: j, I,nmodes
1291 real*8 :: dlgsg, dlgsp, orism5, ndl, nds, super, ccn_at_s
1292 real*8, dimension(nmodes) :: sig_par,sg_par,tp_par
1293
1294 ndl = 0d0
1295
1296 do j=1, nmodes
1297
1298 dlgsg = sig_par(j)
1299
1300 if (sg_par(j) .gt. 0.0) then
1301 if (super .gt. 0.0) then
1302 dlgsp = dlog(sg_par(j)/super)
1303 else
1304 dlgsp = dlog(sg_par(j)/0.01)
1305 end if
1306 else
1307 dlgsp = 0.0
1308
1309 end if
1310 orism5 = 2.d0*dlgsp/(3d0*sq2_par*dlgsg)
1311 ndl = (tp_par(j)/2.0)*(1.0-erf(orism5))+ndl
1312
1313 end do
1314
1315 ccn_at_s = ndl
1316
1317 end subroutine ccn_at_super
1318
1319!=======================================================================
1320
1321! subroutine ccnspec
1322!------------------------------
1323! DESCRIPTION
1324!
1325! *** subroutine ccnspec
1326! *** this subroutine calculates the ccn spectrum of the aerosol using
1327! the appropriate form of kohler theory and assuming a lognormal aerosol size dist
1328!
1329! *** written by athanasios nenes
1330!!
1331! Code Developer
1332! Donifan Barahona
1333! donifanb@umbc.edu
1334!=======================================================================
1335!
1336 subroutine ccnspec (tparc,pparc,nmodes,
1337 & amfs_par,dens_par,deni_par,vhf_par,ams_par,
1338 & sg_par,tp_par,akoh_par,surt_par,temp_par,pres_par
1339 & ,dv_par,act_param,aka_par, psat_par,dair_par,
1340 & ntot,dpg_par)
1341!
1342
1343 integer :: nmodes,i,j,k
1344 real*8 :: tparc, pparc, amfi,denp,vlfs,par1, par2
1345 real*8, dimension(nmodes)::amfs_par,dens_par,deni_par,
1346 & vhf_par,ams_par,sg_par,tp_par,dpg_par
1347 real*8 akoh_par,surt_par,temp_par,pres_par,
1348 & aka_par, dv_par, psat_par,dair_par,ntot
1349 integer act_param
1350
1351
1352 ntot = zero_par
1353 temp_par = max(tparc, 245.0)
1354 pres_par = max(pparc, 34000.0)
1355!
1356 call props(pres_par,temp_par,surt_par,dv_par,act_param,
1357 & aka_par, psat_par,dair_par)
1358
1359
1360!
1361 akoh_par = 4d0*amw_par*surt_par/rgas_par/temp_par/denw_par
1362!
1363
1364 do k=1,nmodes
1365
1366
1367 amfi = max(1.0d0-amfs_par(k),zero_par)
1368 denp = amfs_par(k)*dens_par(k) + amfi*deni_par(k)
1369 vlfs = amfs_par(k)/dens_par(k)/(amfs_par(k)/dens_par(k)+ amfi/
1370 &deni_par(k))
1371 par1 = 4d0*denw_par*ams_par(k)/27d0/vhf_par(k)/denp/amw_par/
1372 & dpg_par(k)**3
1373
1374
1375 par1 = par1/vlfs
1376 par2 = sqrt(max(par1*akoh_par*akoh_par*akoh_par, zero_par))
1377 sg_par(k)= max(exp(par2) - 1d0, zero_par)
1378 ntot = ntot + tp_par(k)
1379 enddo
1380!
1381
1382
1383! *** end of subroutine ccnspec ****************************************
1384!
1385 return
1386 end subroutine ccnspec
1387
1388
1389!=======================================================================
1390! -------------------------------
1391! DESCRIPTION
1392!
1405 subroutine pdfactiv (wparc,sigw, nact,smax,nmodes,
1406 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,
1407 & temp_par, pres_par,aka_par, dv_par, psat_par,
1408 & dair_par,ntot,sig_par)
1409!
1410!
1411 integer :: i, isec,nmodes
1412 real*8 :: tpart, nact, nacti,wparc,smax,ntot, tx1
1413 real*8 :: pdf, dpnmx,sigw,plimt,probi,whi,wlo, scal, wpi,smaxi,
1414 & alfa_par,bet2_par,akoh_par,sg_par(nmodes),tp_par(nmodes),
1415 & tparc,pparc,temp_par, pres_par,aka_par, dv_par, psat_par,
1416 & dair_par,sig_par(nmodes)
1417 real*8 :: wpdbg(npgauss), pddbg(npgauss), nadbg(npgauss),
1418 & smdbg(npgauss)
1419!
1420! *** case where updraft is very small
1421!
1422 if (wparc <= 1d-6) then
1423 smax = 0d0
1424 nact = 0d0
1425 isec = 1
1426 dpnmx = great_par
1427 return
1428 endif
1429!
1430! *** single updraft case
1431!
1432 if (sigw < 1e-10) then
1433 call activate (wparc,nact,smax,nmodes,
1434 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,
1435 & temp_par, pres_par,aka_par, dv_par, psat_par,
1436 & dair_par,ntot,sig_par)
1437 wpdbg(1) = wparc
1438 pddbg(1) = 1.0
1439 nadbg(1) = nact
1440 smdbg(1) = smax
1441!
1442! *** pdf of updrafts
1443!
1444 else
1445 nact = zero_par
1446 smax = zero_par
1447 plimt = 1e-3
1448 probi = sqrt(-2.0*log(plimt*sigw*sq2pi_par))
1449 whi = wparc + sigw*probi
1450 wlo = 0.05
1451 scal = 0.5*(whi-wlo)
1452 do i=1,npgauss
1453 wpi = wlo + scal*(1.0-xgs_par(i))
1454 call activate (wpi,nacti,smaxi,nmodes,
1455 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,
1456 & temp_par, pres_par,aka_par, dv_par, psat_par,
1457 & dair_par,ntot,sig_par)
1458 tx1 = (wpi-wparc)/ sigw
1459 pdf = (1.0/sq2pi_par/sigw) * exp(-0.5*tx1*tx1)
1460 nact = nact + wgs_par(i)*(pdf*nacti)
1461 smax = smax + wgs_par(i)*(pdf*smaxi)
1462 wpdbg(i) = wpi
1463 pddbg(i) = pdf
1464 nadbg(i) = nacti
1465 smdbg(i) = smaxi
1466 if (pdf < plimt) exit
1467 enddo
1468 nact = nact*scal
1469 smax = smax*scal
1470 endif
1471!
1472 return
1473!
1474! *** end of subroutine pdfactiv ****************************************
1475!
1476 end subroutine pdfactiv
1477
1478
1479
1480!=======================================================================
1481
1482! DESCRIPTION
1483!
1484! *** subroutine activate
1485! *** this subroutine calculates the ccn activation fraction according
1486! to the nenes and seinfeld (2003) parameterization, with
1487! modification for non-contunuum effects as proposed by fountoukis
1488! and nenes (2005).
1489!
1490
1491!=======================================================================
1492!
1493 subroutine activate (wparc,ndroplet,smax,nmodes,
1494 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,temp_par, pres_par,
1495 & aka_par, dv_par, psat_par,dair_par,ntot,sig_par)
1496
1497!
1498 integer :: i,niter,nmodes
1499
1500 real*8 :: ndrpl, wparc, beta,cf1, cf2,x1,sinteg1,sinteg2,y1, x2,
1501 &y2,x3,y3, sign,smax, ent_par, ndroplet, nrdpl,bet1_par,
1502 & alfa_par,bet2_par,akoh_par,sg_par(nmodes),tp_par(nmodes),
1503 & temp_par, pres_par,aka_par, dv_par, psat_par,dair_par,ntot,
1504 & wparcel,sig_par(nmodes)
1505!
1506! *** setup common block variables
1507!
1508 wparcel = wparc
1509
1510 sinteg1=zero_par
1511 sinteg2=zero_par
1512 nrdpl = zero_par
1513
1514! *** setup constants
1515!
1516 alfa_par = grav_par*amw_par*dhv_par/cpair_par/rgas_par/temp_par /
1517 &temp_par -grav_par*ama_par/rgas_par/temp_par
1518
1519 bet1_par = pres_par*ama_par/psat_par/amw_par + amw_par*dhv_par*
1520 & dhv_par/cpair_par/rgas_par/temp_par/temp_par
1521 bet2_par = rgas_par*temp_par*denw_par/psat_par/dv_par/amw_par/
1522 & 4d0 +dhv_par*denw_par/4d0/aka_par/temp_par*(dhv_par* amw_par/
1523 &rgas_par/temp_par - 1d0)
1524 beta = 0.5d0*pi_par*bet1_par*denw_par/bet2_par/alfa_par/ wparc/
1525 &dair_par
1526
1527
1528 cf1 = 0.5d0*(((1.d0/bet2_par)/(alfa_par*wparc))**0.5d0)
1529 cf2 = akoh_par/3d0
1530
1531!
1532! *** INITIAL VALUES FOR BISECTION **************************************
1533
1534! *** initial values for bisection **************************************
1535!
1536 x1 = 1.0d-5
1537
1538 if (ntot .gt. zero_par) then
1539 call sintegral (x1,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
1540 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
1541 end if
1542
1543 y1 = (sinteg1*cf1+sinteg2*cf2)*beta*x1 - 1d0
1544!
1545 x2 = 1d0
1546 if (ntot .gt. zero_par) then
1547 call sintegral (x2,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
1548 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
1549 end if
1550
1551 y2 = (sinteg1*cf1+sinteg2*cf2)*beta*x2 - 1d0
1552!
1553! *** perform bisection *************************************************
1554!
155520 do 30 i=1,maxit_par
1556 x3 = 0.5*(x1+x2)
1557!
1558 if (ntot .gt. zero_par) then
1559 call sintegral (x3,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
1560 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
1561 end if
1562
1563 y3 = (sinteg1*cf1+sinteg2*cf2)*beta*x3 - 1d0
1564
1565 if (sign(1.d0,y1)*sign(1.d0,y3) .le. zero_par) then
1566! ! (y1*y3 .le. zero)
1567 y2 = y3
1568 x2 = x3
1569 else
1570 y1 = y3
1571 x1 = x3
1572 endif
1573!
1574 if (abs(x2-x1) .le. eps_par*x1) goto 40
1575 niter = i
1576
1577
157830 continue
1579!
1580! *** converged ; return ************************************************
1581!
158240 x3 = 0.5*(x1+x2)
1583!
1584 if (ntot .gt. zero_par) then
1585 call sintegral (x2,ndrpl,sinteg1,sinteg2,wparcel,nmodes,
1586 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
1587 end if
1588
1589
1590 y3 = (sinteg1*cf1+sinteg2*cf2)*beta*x3 - 1d0
1591 smax = x3
1592
1593 ndroplet=ndrpl
1594 return
1595!
1596! *** end of subroutine activate ****************************************
1597!
1598 end subroutine activate
1599
1600
1601!=======================================================================
1602!
1603! *** subroutine sintegral
1604! *** this subroutine calculates the condensation integrals, according
1605! to the population splitting algorithm of nenes and seinfeld (2003)
1606! modal formulation according to fountoukis and nenes (2004)
1607!
1608! *** written by athanasios nenes
1609!
1610!! Code Developer
1611! Donifan Barahona
1612! donifanb@umbc.edu
1613!=======================================================================
1614!
1615 subroutine sintegral (spar, summa, sum, summat,wparcel,nmodes,
1616 & alfa_par,bet2_par,akoh_par,sg_par,tp_par,sig_par)
1617
1618!
1619 integer ::j,i,k,nmodes
1620 real*8 :: sum, summat, summa, nd(nsmx_par), integ1(nsmx_par),
1621 & integ2(nsmx_par),alfa_par,bet2_par,akoh_par,sg_par(nmodes),
1622 & tp_par(nmodes),sig_par(nmodes)
1623 real*8 :: descr,spar,ratio, ssplt2, ssplt1, sqrt, ssplt, sqtwo,
1624 & dlgsg,dlgsp, ekth,wparcel
1625 real*8 :: orism1, orism2, orism3, orism4, orism5
1626!
1627! *** here is where the criterion with the descriminant is put. when it
1628! is < 0, then set crit2 = .true. otherwise, set the two values of
1629! ssplt and continue.
1630!
1631 descr = 1d0 - (16d0/9d0)*alfa_par*wparcel*bet2_par* (akoh_par/
1632 &spar**2)**2
1633!
1634 if (descr.le.0d0) then
1635 ratio = (2.0d7/3.0)*akoh_par*spar**(-0.3824)
1636 if (ratio.gt.1.0) ratio = 1.0
1637 ssplt2 = spar*ratio
1638 else
1639 ssplt1 = 0.5d0*(1d0-sqrt(descr))
1640 ssplt2 = 0.5d0*(1d0+sqrt(descr))
1641 ssplt1 = sqrt(ssplt1)*spar
1642 ssplt2 = sqrt(ssplt2)*spar
1643 endif
1644!
1645 ssplt = ssplt2
1646!
1647! *** calculate integrals
1648!
1649 sum = 0
1650 summat = 0
1651 summa = 0
1652!
1653 sqtwo = sqrt(2d0)
1654
1655!
1656 do 999 j = 1, nmodes
1657!
1658 dlgsg = sig_par(j)
1659 dlgsp = dlog(sg_par(j)/spar)
1660 orism1 = 2.d0*dlog(sg_par(j)/ssplt2)/(3.d0*sqtwo*dlgsg )
1661 orism2 = orism1 - 3.d0*dlgsg/(2.d0*sqtwo)
1662 orism3 = 2.d0*dlgsp/(3.d0*sqtwo*dlgsg)-3.d0*dlgsg/(2.d0*sqtwo)
1663 orism4 = orism1 + 3.d0*dlgsg/sqtwo
1664 orism5 = 2.d0*dlgsp/(3*sqtwo*dlgsg)
1665 ekth = exp(9d0/2d0*dlgsg*dlgsg)
1666 integ1(j) = tp_par(j)*spar*((1-erf(orism1)) - 0.5d0*((sg_par(j)/
1667 &spar)**2)*ekth*(1-erf(orism4)))
1668
1669 integ2(j) = (exp(9d0/8d0*dlgsg*dlgsg)*tp_par(j)/sg_par(j))*
1670 & (erf(orism2) - erf(orism3))
1671!
1672! *** calculate number of drops
1673!
1674 nd(j) = (tp_par(j)/2.0)*(1.0-erf(orism5))
1675!
1676 sum = sum + integ1(j)
1677 summat = summat + integ2(j)
1678 summa = summa + nd(j)
1679999 continue
1680!
1681 return
1682 end subroutine sintegral
1683
1684!=======================================================================
1685
1686! DESCRIPTION
1687!
1688! *** subroutine props
1689! *** this subroutine calculates the thermophysical properties for the CCN activ param
1690!
1691! *** written by athanasios nenes
1692! Code Developer
1693! Donifan Barahona
1694! donifanb@umbc.edu
1695
1696!=======================================================================
1697!
1698 subroutine props(pres_par,temp_par,surt_par,dv_par,act_param,
1699 & aka_par, psat_par,dair_par)
1700!
1701!
1702 real*8 :: presa,dbig,dlow,coef,aka_par, dv_par, psat_par
1703 real*8 :: pres_par,temp_par,surt_par,dair_par
1704 integer act_param
1705!
1706 denw_par = 1d3
1707 dhv_par = 2.25d6
1708 cpair_par = 1.0061d3
1709 presa = pres_par/1.013d5
1710 dair_par = pres_par*ama_par/rgas_par/temp_par
1711 aka_par = (4.39+0.071*temp_par)*1d-3
1712 surt_par = sft(temp_par)
1713
1714 if (act_param .le. 1) then
1715 dv_par = (0.211d0/presa)*(temp_par/273d0)**1.94
1716 dv_par = dv_par*1d-4
1717 dbig = 5.0d-6
1718 dlow = 0.207683*((accom_par)**(-0.33048))
1719 dlow = dlow*1d-6
1720
1721
1722
1723 coef = ((2*pi_par*amw_par/(rgas_par*temp_par))**0.5)
1724
1725 dv_par = (dv_par/(dbig-dlow))*((dbig-dlow)-(2*dv_par/accom_par) *
1726 &coef*(dlog((dbig+(2*dv_par/accom_par)*coef)/(dlow+ (2*dv_par/
1727 &accom_par)*coef))))
1728
1729 psat_par = vpres(temp_par)*(1e5/1.0d3)
1730
1731
1732
1733 end if
1734!
1735 return
1736!
1737! *** end of subroutine props *******************************************
1738!
1739 end subroutine props
1740
1741
1742
1743!PHYSICAL PROPERTIES for Nenes CDNC Activation
1744
1745
1746
1747!=======================================================================
1748!
1749! *** function vpres
1750! *** this function calculates saturated water vapour pressure as a
1751! function of temperature. valid for temperatures between -50 and
1752! 50 c.
1753!
1754! ======================== arguments / usage ===========================
1755!
1756! input:
1757! [t]
1758! real variable.
1759! ambient temperature expressed in kelvin.
1760!
1761! output:
1762! [vpres]
1763! real variable.
1764! saturated vapor pressure expressed in mbar.
1765!
1766!=======================================================================
1767!
1768 real*8 function vpres (t)
1769!
1770 integer ::i
1771 real*8 :: a(0:6), t,ttemp, vp
1772 data a/6.107799610e+0, 4.436518521e-1, 1.428945805e-2,
1773 & 2.650648471e-4, 3.031240396e-6, 2.034080948e-8, 6.136820929e-11/
1774!
1775! calculate polynomial (without exponentiation).
1776!
1777 ttemp = t-273.0d0
1778 vp = a(6)*ttemp
1779 do i=5,1,-1
1780 vp = (vp + a(i))*ttemp
1781 enddo
1782 vpres = vp + a(0)
1783!
1784! end of function vpres
1785!
1786 return
1787 end function vpres
1788
1789
1790
1791!=======================================================================
1792!
1793! *** function sft
1794! *** this function calculates water surface tension as a
1795! function of temperature. valid for temperatures between -40 and
1796! 40 c.
1797!
1798! ======================== arguments / usage ===========================
1799!
1800! input:
1801! [t]
1802! real variable.
1803! ambient temperature expressed in kelvin.
1804!
1805! output:
1806! [sft]
1807! real variable.
1808! surface tension expressed in j m-2.
1809!
1810!=======================================================================
1811!
1812 real*8 function sft (t)
1813!
1814 implicit none
1815
1816!
1817 real*8 :: t,tpars
1818!
1819 tpars = t-273.15d0
1820 sft = 0.0761-1.55e-4*tpars
1821!
1822 return
1823 end function sft
1824
1825
1826! ***********************************************************************
1827!
1828 subroutine gauleg (x,w,n)
1829!
1830! calculation of points and weights for n point gauss integration
1831! ***********************************************************************
1832!
1833 integer :: n,m,i,j
1834 real*8 :: x(n), w(n),xm,xl,z,p1,p2,p3,pp,z1
1835 real*8, parameter :: eps_par=1.e-6
1836 real*8, parameter :: x1=-1.0, x2=1.0
1837!
1838! calculation
1839!
1840 m=(n+1)/2d0
1841 xm=0.5d0*(x2+x1)
1842 xl=0.5d0*(x2-x1)
1843 do 12 i=1,m
1844 z=cos(pi_par*(i-.25d0)/(n+.5d0))
1845 1 continue
1846 p1=1.d0
1847 p2=0.d0
1848 do 11 j=1,n
1849 p3=p2
1850 p2=p1
1851 p1=((2.d0*j-1.)*z*p2-(j-1.d0)*p3)/j
1852 11 continue
1853 pp=n*(z*p1-p2)/(z*z-1.d0)
1854 z1=z
1855 z=z1-p1/pp
1856 if(abs(z-z1).gt.eps_par)go to 1
1857 x(i)=xm-xl*z
1858 x(n+1-i)=xm+xl*z
1859 w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
1860 w(n+1-i)=w(i)
1861 12 continue
1862 return
1863 end subroutine gauleg
1864
1865!C=======================================================================
1866!C
1867!C *** REAL FUNCTION erf (overwrites previous versions)
1868!C *** THIS SUBROUTINE CALCULATES THE ERROR FUNCTION USING A
1869!C *** POLYNOMIAL APPROXIMATION
1870!C
1871!C=======================================================================
1872!C
1873 REAL*8 FUNCTION erf(x)
1874 real*8 :: x
1875 real*8 :: aa(4), axx, y
1876 DATA aa /0.278393d0,0.230389d0,0.000972d0,0.078108d0/
1877
1878 y = dabs(dble(x))
1879 axx = 1.d0 + y*(aa(1)+y*(aa(2)+y*(aa(3)+y*aa(4))))
1880 axx = axx*axx
1881 axx = axx*axx
1882 axx = 1.d0 - (1.d0/axx)
1883 if(x.le.0.) then
1884 erf = -axx
1885 else
1886 erf = axx
1887 endif
1888 RETURN
1889 END FUNCTION
1890
1891
1892!=======================================================================
1893!
1894! *** real function erf
1895! *** this subroutine calculates the error function
1896!
1897! *** obtained from numerical recipies
1898!
1899!=======================================================================
1900!
1901
1902!
1903! real*8 :: x
1904
1905! if(x.lt.0.)then
1906! erf=-gammp(.5d0,x**2)
1907! else
1908! erf=gammp(.5d0,x**2)
1909! endif
1910! return
1911
1912! end function erf
1913
1914!
1915!=======================================================================
1916!
1917 real*8 function gammln(xx)
1918!
1919!=======================================================================
1920!
1921!
1922 integer :: j
1923 real*8 :: cof(6),stp,half,one,fpf,x,tmp,ser,xx
1924!
1925 data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0, -
1926 &1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
1927 data half,one,fpf/0.5d0,1.0d0,5.5d0/
1928 x=xx-one
1929 tmp=x+fpf
1930 tmp=(x+half)*log(tmp)-tmp
1931 ser=one
1932 do 11 j=1,6
1933 x=x+one
1934 ser=ser+cof(j)/x
1935 11 continue
1936 gammln=tmp+log(stp*ser)
1937 return
1938 end function gammln
1939
1940
1941!
1942!=======================================================================
1943!
1944! real*8 function gammp(a,x)
1945!
1946!=======================================================================
1947!
1948! real*8 :: a,x,gln,gamser,gammcf
1949
1950! if(x.lt.0.d0.or.a.le.0.d0)pause
1951! if(x.lt.a+1.d0)then
1952! call gser(gamser,a,x,gln)
1953! gammp=gamser
1954! else
1955! call gcf(gammcf,a,x,gln)
1956! gammp=1.d0-gammcf
1957! endif
1958! return
1959! end function gammp
1960
1961
1962!
1963!=======================================================================
1964!
1965! subroutine gcf(gammcf,a,x,gln)
1966!
1967!=======================================================================
1968!
1969!
1970! integer :: n
1971! integer, parameter :: itmax=100
1972! real*8, parameter :: eps_par=3.e-7
1973! real*8 :: gln, gold,a0,a1,x,b0,b1,fac,an, &
1974! float,ana,a,anf,g,gammcf
1975! gln=gammln(a)
1976! gold=0.
1977! a0=1.
1978! a1=x
1979! b0=0.
1980! b1=1.
1981! fac=1.
1982! do 11 n=1,itmax
1983! an=float(n)
1984! ana=an-a
1985! a0=(a1+a0*ana)*fac
1986! b0=(b1+b0*ana)*fac
1987! anf=an*fac
1988! a1=x*a0+anf*a1
1989! b1=x*b0+anf*b1
1990! if(a1.ne.0.)then
1991! fac=1./a1
1992! g=b1*fac
1993! if(abs((g-gold)/g).lt.eps_par)go to 1
1994! gold=g
1995! endif
1996!1 continue
1997! pause 'a too large, itmax too small'
1998! gammcf=exp(-x+a*log(x)-gln)*g
1999! return
2000! end subroutine gcf
2001
2002
2003!
2004!=======================================================================
2005!sft
2006! subroutine gser(gamser,a,x,gln)
2007!
2008!=======================================================================
2009!
2010
2011! integer :: n
2012! integer, parameter :: itmax=100
2013! real*8, parameter :: eps_par=3.e-7
2014! real*8 :: gln,x,gamser,a,ap,sum,del,abs
2015
2016! gln=gammln(a)
2017! if(x.le.0.)then
2018! if(x.lt.0.)pause
2019! gamser=0.
2020! return
2021! endif
2022! ap=a
2023! sum=1./a
2024! del=sum
2025! do 11 n=1,itmax
2026! ap=ap+1.
2027! del=del*x/ap
2028! sum=sum+del
2029! if(abs(del).lt.abs(sum)*eps_par)go to 1
2030!1 continue
2031! pause 'a too large, itmax too small'
2032! gamser=sum*exp(-x+a*log(x)-gln)
2033! return
2034! end subroutine gser
2035
2036! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2037!
2038
2039!
2040
2041!*************************************************************************
2042!
2043! ICE NUCLEATION PARAMETERIZATION FILES START HERE
2044!
2045! ************************************************************************
2046!
2047!
2048!======================================================================
2049!
2050! Code Developer
2051! Donifan Barahona, GA TECH
2052! donifan@gatech.edu
2053! -------------------------------
2054! DESCRIPTION
2055!
2056
2057!***********************************************************
2058!** Parameterization of ICE crystal number concentration
2059!** for large scale models.
2060!** Donifan Barahona, Athanasios Nenes
2061! JGR, 111, D11211, 2008
2062! ACP, 9, 369-381, 2009
2063! ACP 9, 5933-5948, 2009
2064! Homogeneoeus and heterogeneous nucleation considered
2065!*** SI units unless otherwise specified.
2066!
2067!
2068! *** WRITTEN BY DONIFAN BARAHONA
2069!
2070!=======================================================================
2071
2072
2073
2074 SUBROUTINE iceparam (sigma_w, denice_ice,ddry_ice,np_ice,
2075 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
2076 & g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,sc_ice,
2077 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,
2078 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2079 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2080 & si0bc_ice,nbc_ice,one_over_tao_preex,
2081 & nhet, nice, smax, nlim,Nhet_dep,Nhet_dhf,fdust_dep
2082 & ,P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2083 & use_av_v,miuv_ice,vpresw_ice,vpresi_ice,denair_ice)
2084
2085 real*8, intent(in) :: sigma_w,denice_ice,ddry_ice,np_ice,
2086 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
2087 & g1_ice, g2_ice,gdoin_ice, z_ice, lambda_ice,
2088 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,
2089 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2090 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2091 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,p_ice, t_ice,
2092 & vpresw_ice,vpresi_ice,denair_ice
2093 real*8, dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
2094
2095
2096 real*8, intent(out) :: nice, nhet, smax, nlim,nhet_dep,nhet_dhf
2097 & ,fdust_dep
2098
2099 real*8 :: wpar_ice, sigmav_ice,vmax_ice,miuv_ice
2100 integer, intent(in)::nbindust_ice
2101 logical use_av_v
2102
2103
2104
2105
2106 vmin_ice=0.005d0
2107 vmax_ice=2.5d0
2108 sigmav_ice=sigma_w
2109 vmax_ice= max(min(miuv_ice+(4d0*sigmav_ice), vmax_ice),
2110 & vmin_ice +0.01)
2111
2112 if ((sigmav_ice .lt. 0.05) .or. (t_ice .gt. thom)) then
2113 use_av_v= .true.
2114 end if
2115
2116 if (vmax_ice .gt. vmin_ice + 0.01) then
2117 if (use_av_v) then
2118 wpar_ice=min(max(miuv_ice + sigma_w*0.8, vmin_ice), vmax_ice)
2119 CALL nice_param(wpar_ice, denice_ice,ddry_ice,np_ice,nin_ice,
2120 & alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2121 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2122 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2123 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2124 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2125 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
2126 & nice, smax, nhet, nlim,nhet_dep,nhet_dhf,fdust_dep,
2127 & p_ice, t_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2128 & vpresw_ice,vpresi_ice,denair_ice)
2129 else
2130 call nice_vdist(denice_ice,ddry_ice,np_ice,
2131 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2132 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2133 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2134 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2135 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2136 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
2137 & nice, smax, nhet, nlim,nhet_dep,nhet_dhf,fdust_dep,
2138 & p_ice, t_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice
2139 & ,vpresw_ice,vpresi_ice,denair_ice)
2140 end IF
2141 else
2142 nice = zero_par
2143 nhet = zero_par
2144 nlim = 1.0e10
2145 smax = zero_par
2146 end if
2147
2148
2149 return
2150 END subroutine iceparam
2151
2152
2153!*************************************************************
2154! Subroutine nice_Vdist. Calculates the ice crystal number concentration
2155! at the maximum supersaturation using a PDF of updraft using a
2156! sixth order Gauss-Legendre quadrature
2157! Inputs: T, and P all SI units)
2158! Output NC, smax, nhet, nlim (m-3)
2159! Barahona and Nenes, JGR, D11211 (2008) and ACPD, 15665-15698, (2008)
2160! Written by Donifan Barahona
2161!************************************************************
2162
2163
2164 subroutine nice_vdist(denice_ice,ddry_ice,np_ice,
2165 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2166 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2167 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2168 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2169 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2170 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
2171 & nice, smax, nhet, nlim,Nhet_dep,Nhet_dhf,fdust_dep,
2172 & P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2173 & vpresw_ice,vpresi_ice,denair_ice)
2174
2175 real*8 :: quadx(6), wpar, sum1, quadw(6), dp, sum2, sum3, sum4,
2176 & sum5, x1, x2
2177 real*8, intent(out) :: nice, smax, nhet, nlim,nhet_dep,nhet_dhf
2178 & ,fdust_dep
2179 real*8 wpar_icex,denice_ice,ddry_ice,np_ice,
2180 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2181 & normv_ice,sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2182 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2183 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,sc_ice,
2184 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2185 & si0bc_ice,nbc_ice,one_over_tao_preex,p_ice, t_ice,
2186 & vpresw_ice,vpresi_ice,denair_ice
2187 real*8, dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
2188
2189 INTEGER :: INDEX,nbindust_ice
2190
2191 DATA quadx/0.23861918d0, -0.23861918d0, 0.66120939d0, -
2192 &0.66120939d0, 0.93246951d0, -0.93246951d0/
2193
2194 DATA quadw/0.46791393d0, 0.46791393d0, 0.36076157d0,
2195 & 0.36076157d0, 0.17132449d0, 0.17132449d0/
2196
2197
2198
2199 x1=(vmin_ice-miuv_ice)/(sq2_par*sigmav_ice)
2200 x2=(vmax_ice-miuv_ice)/(sq2_par*sigmav_ice)
2201
2202 x2=max(x1 +0.01, x2)
2203 normv_ice=(erfapp(x2)-erfapp(x1))*0.5d0
2204
2205 sum1=0d0
2206 sum2=0d0
2207 sum3=0d0
2208 sum4=0d0
2209 sum5=0d0
2210
2211 DO index =1, 6
2212 wpar=max(0.5d0*(((vmax_ice-vmin_ice)*quadx(index)) +(vmax_ice+
2213 & vmin_ice)), 0.01)
2214
2215
2216 CALL nice_param(wpar, denice_ice,ddry_ice,np_ice,nin_ice,
2217 & alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2218 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2219 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2220 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2221 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2222 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
2223 & nice, smax, nhet, nlim,nhet_dep,nhet_dhf,fdust_dep,
2224 & p_ice, t_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2225 & vpresw_ice,vpresi_ice,denair_ice)
2226
2227 CALL gausspdf(wpar, dp, sigmav_ice,miuv_ice,normv_ice)
2228 sum1=sum1+(nice*dp*quadw(index))
2229 sum2=sum2+(smax*dp*quadw(index))
2230 sum3=sum3+(nhet*dp*quadw(index))
2231 sum4=sum4+(nlim*dp*quadw(index))
2232 sum5=sum5+(sc_ice*dp*quadw(index))
2233
2234
2235
2236 END DO
2237 nice=sum1*(vmax_ice-vmin_ice)*0.5d0
2238 smax=sum2*(vmax_ice-vmin_ice)*0.5d0
2239 nhet=sum3*(vmax_ice-vmin_ice)*0.5d0
2240 nlim=sum4*(vmax_ice-vmin_ice)*0.5d0
2241 sc_ice=sum5*(vmax_ice-vmin_ice)*0.5d0
2242 RETURN
2243
2244
2245 END subroutine nice_vdist
2246
2247!*************************************************************
2248
2249!*************************************************************
2250 real*8 function erfapp(x)
2251
2252 real*8, intent(in) :: x
2253 real*8 :: a
2254 a = x*x*(1.27324d0+(0.147d0*x*x))/(1d0+(0.147d0*x*x))
2255 erfapp = sqrt(1d0-exp(-a))
2256
2257 if (x .lt. 0.0) then
2258 erfapp = - erfapp
2259 end if
2260
2261 end function erfapp
2262
2263
2264
2265!*************************************************************
2266! Subroutine Het_freezing. Use only for mixed phase clouds . Inputs: Wpar, T, and P all SI units)
2267! Output Nc=Nhet (m-3)
2268! ! Written by Donifan Barahona
2269!************************************************************
2270
2271
2272! SUBROUTINE Het_freezing(nhet, nice, smax)
2273
2274! real*8, intent(out) :: nice, nhet, smax
2275! real*8, intent(in) :: np_ice
2276
2277! real*8 :: I, SX, NHET_, DSH
2278
2279! nhet=zero_par
2280! smax=zero_par
2281! nice=zero_par
2282
2283! SX=vpresw_ice/vpresi_ice
2284! call INSPEC_ice(SX-1.0, NHET_, DSH,np_ice)
2285! sc_ice=1.0
2286! nhet=max(NHET_, zero_par)
2287! nice=nhet
2288! smax=max(SX-1.0, zero_par)
2289
2290
2291! END SUBROUTINE Het_freezing
2292
2293
2294!*************************************************************
2295! Subroutine nice_param. This is the implementation of the Barahona and Nenes(2008, 2009)
2296! parameterization
2297! Written by Donifan Barahona
2298!************************************************************
2299
2300
2301 subroutine nice_param(wpar_icex,denice_ice,ddry_ice,np_ice,
2302 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2303 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2304 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,lambda_ice,
2305 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2306 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2307 & si0bc_ice,nbc_ice,sc_ice,one_over_tao_preex,
2308 & nice, smax, nhet, nlim_,Nhet_dep,Nhet_dhf,fdust_dep,
2309 & P_ice, T_ice,ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2310 & vpresw_ice,vpresi_ice,denair_ice)
2311
2312 real*8, intent(in) :: wpar_icex,denice_ice,np_ice,nin_ice,
2313 & alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,miuv_ice,
2314 & sigmav_ice,g1_ice, g2_ice,gdoin_ice, z_ice,vmax_ice,
2315 & norg_ice,sigorg_ice,dorg_ice, dbc_ice,sigbc_ice,
2316 & kdust_ice, kbc_ice, shdust_ice, shbc_ice, ddry_ice,
2317 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2318 & si0bc_ice,nbc_ice,one_over_tao_preex,p_ice, t_ice,
2319 & vpresw_ice,vpresi_ice,denair_ice
2320 real*8, intent(out) :: nice, smax, nhet, nlim_,nhet_dep,nhet_dhf
2321 & ,fdust_dep
2322
2323 real*8 :: aux1, aux2, g, dpmax, miu, monvol, fds, nlim, dlim,
2324 & dstar, ds, nstar, nhom, fc, phido, auxnc, sizecorr, dsh, nhet_,
2325 & f1, f2, saux, sup, slow, shalf, fhalf, dpmin, gam, wpar_ice,
2326 & preex_effect, swsat, sstep,sc_ice,lambda_ice
2327 real*8, dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
2328
2329 integer :: INDEX, maxiter_s,nbindust_ice
2330
2331
2332
2333
2334
2335 preex_effect=1.0-(one_over_tao_preex*(shom_ice)/alfa_ice/
2336 & ((shom_ice+1.0))/wpar_icex)
2337 swsat = vpresw_ice/vpresi_ice -1d0
2338
2339
2340 if (preex_effect .le. 0.0) then
2341 nhet=0d0
2342 nhom=0d0
2343 smax=shom_ice
2344 dsh =0.d0
2345 fds=1.d0
2346
2347 sc_ice = shom_ice + 1.d0
2348 nice = 0.d0
2349 nlim_=0d0
2350 return
2351 else
2352
2353 wpar_ice = wpar_icex*preex_effect
2354 end if
2355
2356
2357
2358
2359 if (np_ice .gt. 1.0) then
2360
2361 monvol=np_ice*1.0d-6*ddry_ice*ddry_ice*ddry_ice
2362 aux1=1.6397d-14*t_ice-3.1769d-12
2363 dpmax=aux1*(monvol**(-0.373d0))*(wpar_ice**(-0.05))
2364 IF (dpmax.GT.1.0d-4) THEN
2365 dpmax=1.0d-4
2366 END IF
2367
2368 else
2369 dpmax=1.0d-4
2370
2371 endif
2372
2373
2374 dpmin=dliq_ice+(0.02/sqrt(alfa_ice*wpar_ice*g1_ice))
2375 dpmax=max(dpmax,dpmin)
2376
2377 aux1=dpmax-dliq_ice
2378 aux2=dlog((g2_ice+(g1_ice*dpmax))/(g2_ice+(g1_ice*dliq_ice)))
2379 g=1.3346d0*((g1_ice*aux1)-(g2_ice*aux2))/(aux1*g1_ice*g1_ice)
2380 lambda_ice=lambda_ice/sqrt(wpar_ice)
2381 aux1 = g1_ice*alfa_ice*wpar_ice
2382 nstar=(aux1*sqrt(aux1))/beta_ice/z_ice/sq2_par
2383
2384 gam=g2_ice/g1_ice
2385
2386
2387
2388
2389 fds=1d0
2390 nhom=0d0
2391
2392
2393
2394 if (typeofspec_ice .ge. 0d0) then
2395
2396 call inspec_ice(shom_ice,nhet_,dsh,np_ice,norg_ice,sigorg_ice,
2397 & dorg_ice, dbc_ice,sigbc_ice,
2398 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2399 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2400 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2401 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice
2402 & ,vpresw_ice,vpresi_ice)
2403 sizecorr=exp(-2d0/lambda_ice/shom_ice)
2404 dstar=((4d0*dsh*dsh/3d0)+(2d0*dsh*(shom_ice-dsh))) /(shom_ice-
2405 &dsh+1d0)
2406
2407 if (dstar .gt. 0.0) then
2408 nlim=min(nstar*(shom_ice+1d0)/shom_ice/sqrt(dstar)/sizecorr,
2409 & 1.e10)
2410 else
2411 nlim=1.e10
2412 end if
2413
2414 else
2415 dsh=shom_ice-sh_ice
2416 dstar=((4d0*dsh*dsh/3d0)+(2d0*dsh*(shom_ice-dsh))) /(shom_ice-
2417 &dsh+1d0)
2418 dlim=-gam+sqrt((gam*gam)+(2d0*dstar/g1_ice/alfa_ice/wpar_ice))
2419 nlim=alfa_ice*wpar_ice*(shom_ice+1d0)/z_ice/beta_ice/shom_ice
2420 nlim=nlim*((g1_ice*dlim)+g2_ice)/dlim/dlim
2421 nhet_=nin_ice
2422 end if
2423
2424
2425 nlim_=min(nlim, 1d10)
2426 nlim_=max(nlim, 1d-6)
2427
2428 if (nhet_ .gt. 0.0) then
2429 aux1 =nhet_/nlim_
2430 fds=1d0-(aux1*sqrt(aux1))
2431 else
2432 fds = 1d0
2433 end if
2434
2435
2436
2437 if (purehom_ice) then
2438 fds=1.0d0
2439 end if
2440
2441 if ((purehet_ice) .or. (t_ice .GE. thom)) then
2442 fds = 0d0
2443 end if
2444
2445
2446 IF (fds .GE. 1.e-10) THEN
2447
2448 miu=fds*alfa_ice*(shom_ice+1d0)*wpar_ice*koft_ice/shom_ice
2449
2450 phido=sqrt(pi_ice*g/miu/2d0)*(g/miu)
2451 auxnc=2d0*denair_ice/beta_ice/koft_ice/denice_ice/pi_ice/np_ice
2452 fc=auxnc/phido
2453
2454
2455
2456 if (np_ice .gt.0d0) then
2457
2458 IF (fc .le. 0.6d0) then
2459 nhom=np_ice*exp(-fc)*(1.0d0-exp(-fc))
2460 else
2461 nhom=np_ice/(1d0+exp((9d0-2d0*fc)/7d0))
2462 end if
2463
2464 else
2465 nhom=0d0
2466 end if
2467
2468 smax=shom_ice
2469 nhet=nhet_
2470 if (purehom_ice) nhet_ = 0.d0
2471
2472
2473 ELSE
2474
2475
2476 nhom = 0d0
2477
2478
2479 smax=0d0
2480 nhet=0d0
2481 saux=0.01d0
2482
2483 if (typeofspec_ice .lt. 0d0) saux=sh_ice+0.00000000001d0
2484
2485 f1= findsmax(saux,dsh,
2486 & nhet_,np_ice,norg_ice,sigorg_ice,
2487 & dorg_ice, dbc_ice,sigbc_ice,
2488 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2489 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2490 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2491 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
2492 & gdoin_ice,alfa_ice,wpar_ice,gam,g1_ice,g2_ice,z_ice,
2493 & beta_ice,nin_ice,vpresw_ice,vpresi_ice,nstar)
2494 f2=1.0
2495 sstep = 0.05d0
2496 maxiter_s = int(1d0/sstep) + 1
2497
2498 do index =1, maxiter_s
2499
2500 if (saux .ge. swsat) then
2501 f2=f1
2502 saux = swsat
2503 exit
2504 end if
2505
2506
2507 saux=saux+sstep
2508 f2=findsmax(saux,dsh,
2509 & nhet_,np_ice,norg_ice,sigorg_ice,
2510 & dorg_ice, dbc_ice,sigbc_ice,
2511 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2512 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2513 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2514 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
2515 & gdoin_ice,alfa_ice,wpar_ice,gam,g1_ice,g2_ice,z_ice,
2516 & beta_ice,nin_ice,vpresw_ice,vpresi_ice,nstar)
2517 IF (f2*f1 .lt. 0d0) exit
2518
2519 end do
2520
2521 if (f2*f1 .gt.0d0) then
2522 nhet=0d0
2523 smax=saux
2524 else
2525
2526 if (saux .lt. swsat) then
2527
2528
2529 sup=saux
2530 slow=saux-(sstep + 0.001d0)
2531
2532 DO index=1,50
2533 shalf=0.5d0*(sup+slow)
2534 fhalf=findsmax(shalf,dsh,
2535 & nhet_,np_ice,norg_ice,sigorg_ice,
2536 & dorg_ice, dbc_ice,sigbc_ice,
2537 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2538 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2539 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2540 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
2541 & gdoin_ice,alfa_ice,wpar_ice,gam,g1_ice,g2_ice,z_ice,
2542 & beta_ice,nin_ice,vpresw_ice,vpresi_ice,nstar)
2543
2544 IF (sign(1.d0,f1)*sign(1.d0,fhalf) .LE. 0d0) THEN
2545 f2 = fhalf
2546 sup = shalf
2547 ELSE
2548 f1 = fhalf
2549 slow = shalf
2550 ENDIF
2551
2552
2553 IF (abs(slow-sup) .LE. 5d-3) exit
2554 END DO
2555 else
2556 shalf = swsat
2557 end if
2558
2559 smax=shalf
2560
2561 end if
2562
2563 if ((smax .gt. shom_ice) .and.(t_ice .LE. thom)) smax=shom_ice
2564
2565 if (typeofspec_ice .ge. 0d0) then
2566 call inspec_ice(smax, nhet,dsh,np_ice,norg_ice,sigorg_ice,
2567 & dorg_ice, dbc_ice,sigbc_ice,
2568 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2569 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2570 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2571 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2572 & vpresw_ice,vpresi_ice)
2573 else
2574 nhet=nin_ice
2575 end if
2576
2577 END IF
2578
2579
2580 nice=nhom+nhet
2581 sc_ice=max(smax+1.0-dsh, 1.0)
2582
2583
2584 if (.true.) then
2585
2586 if (fds .gt. 0.0) then
2587 sc_ice = (shom_ice+1.0)*fds + sc_ice*(1.0-fds)
2588 end if
2589 end if
2590
2591 sc_ice=min(shom_ice+1.0, sc_ice)
2592
2593
2594 return
2595
2596 END subroutine nice_param
2597!*************************************************************
2598 real*8 function findsmax(SX,DSH,
2599 & NHET_,np_ice,norg_ice,sigorg_ice,
2600 & dorg_ice, dbc_ice,sigbc_ice,
2601 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2602 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2603 & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
2604 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,lambda_ice,
2605 & gdoin_ice,alfa_ice,wpar_ice,GAM,g1_ice,g2_ice,z_ice,
2606 & beta_ice,nin_ice,vpresw_ice,vpresi_ice,NSTAR)
2607 real*8, intent(in) :: sx
2608 real*8 :: tao
2609 real*8 :: dstar, sizecorr, dsh, nstar,dlim
2610 integer nbindust_ice
2611 real*8, dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
2612 real*8 shom_ice,nhet_,np_ice,norg_ice,sigorg_ice,
2613 & dorg_ice, dbc_ice,sigbc_ice,
2614 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2615 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2616 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2617 & lambda_ice,
2618 & gdoin_ice,alfa_ice,wpar_ice,gam,g1_ice,g2_ice,z_ice,
2619 & beta_ice,nin_ice,vpresw_ice,vpresi_ice
2620
2621 if (typeofspec_ice .ge. 0d0) then
2622
2623 call inspec_ice(sx,nhet_,dsh,np_ice,norg_ice,sigorg_ice,
2624 & dorg_ice, dbc_ice,sigbc_ice,
2625 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2626 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2627 & si0bc_ice,nbc_ice,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice,
2628 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
2629 & vpresw_ice,vpresi_ice)
2630 sizecorr=exp(-2d0/lambda_ice/sx)
2631 dstar=((4d0*dsh*dsh/3d0)+(2d0*dsh*(sx-dsh)))/(sx-dsh+1d0)
2632 dstar=dstar+(gdoin_ice*alfa_ice*wpar_ice)
2633
2634
2635 tao=nhet_*sizecorr*sx*sqrt(dstar)/(sx+1d0)/nstar
2636
2637
2638 else
2639 dsh=sx-sh_ice
2640 dstar=((4d0*dsh*dsh/3d0)+(2d0*dsh*(sx-dsh))) /(sx-dsh+1d0)
2641 dlim=-gam+sqrt((gam*gam)+(2d0*dstar/g1_ice/alfa_ice/wpar_ice))
2642 tao=alfa_ice*wpar_ice*(sx+1d0)/z_ice/beta_ice/sx
2643 tao=tao*((g1_ice*dlim)+g2_ice)/dlim/dlim/nin_ice
2644
2645 end if
2646
2647
2648 findsmax=1d0-tao
2649
2650 end function findsmax
2651
2652
2653!*************************************************************
2654! Function VPRESWATER. Calculates the saturated vapor pressure
2655! of water (Pa) according to Murphy & Koop (2005)
2656! T in K (173.15-373.15)
2657!************************************************************
2658
2659 real*8 function vpreswater_ice(T)
2660
2661 real*8, intent(in) :: t
2662 real*8 :: a(0:9)
2663
2664 DATA a/54.842763d0, -6763.22d0, -4.21d0, 0.000367d0, 0.0415d0,
2665 & 218.8d0, 53.878d0, -1331.22d0, -9.44523d0, 0.014025d0/
2666
2667
2668 vpreswater_ice = a(0)+(a(1)/t)+(a(2)*log(t))+(a(3)*t)+
2669 & (tanh(a(4)*(t-a(5)))*((a(6)+(a(7)/t))+ (a(8)*log(t))+ (a(9)*t)))
2670
2671 vpreswater_ice=exp(vpreswater_ice)
2672
2673 return
2674 END function vpreswater_ice
2675
2676!*************************************************************
2677! Function VPRESICE. Calculates the saturated vapor pressure
2678! of ice (pa) according to Murphy & Koop (2005)
2679! T in K (>110)
2680!************************************************************
2681
2682 real*8 function vpresice(T)
2683
2684 real*8, intent(in) :: t
2685 real*8 :: a(0:3)
2686
2687 DATA a/9.550426d0, -5723.265d0, 3.53068d0, -0.00728332d0/
2688
2689
2690 vpresice = a(0)+(a(1)/t)+(a(2)*log(t))+(a(3)*t)
2691 vpresice=exp(vpresice)
2692
2693 return
2694 END function vpresice
2695
2696!*************************************************************
2697! Function DHSUB. Calculates the latent heat of sublimation
2698! of ice (J/Kg) according to Murphy & Koop (2005)
2699! T in K (>30)
2700!************************************************************
2701
2702 real*8 function dhsub_ice(T)
2703
2704 real*8, intent(in) :: t
2705 real*8 :: a(0:4)
2706
2707
2708 DATA a/46782.5d0, 35.8925d0, -0.07414d0, 541.5d0, 123.75d0/
2709
2710 dhsub_ice = a(0) + (a(1) * t) + (a(2)*t*t) + (a(3) * exp(-((t/
2711 & a(4))**2)))
2712
2713 dhsub_ice=1000d0*dhsub_ice/18d0
2714 return
2715 END function dhsub_ice
2716
2717!*************************************************************
2718! Function ICEDENSITY. Calculates the DENSITY OF ICE
2719! of ice (Kg/m3) according to PK97
2720! T in K (>30)
2721!************************************************************
2722
2723 real*8 function densityice(T)
2724
2725 real*8, intent(in) :: t
2726 real*8 :: a(0:2), ttemp
2727
2728 DATA a/0.9167d0, -1.75d-4, -5.0d-7/
2729
2730 ttemp=t-273d0
2731
2732 densityice= 1000d0*(a(0)+(a(1)*ttemp)+(a(2)*ttemp*ttemp))
2733 return
2734 END function densityice
2735
2736!*************************************************************
2737! Function WATDENSITY. Calculates the DENSITY OF ICE
2738! of liquid water (Kg/m3) according to PK97
2739! T in K (>240)
2740!************************************************************
2741
2742 real*8 function watdensity_ice(T)
2743
2744 real*8, intent(in) :: t
2745 real*8 :: a(0:6), ttemp, watdensity
2746 INTEGER :: i
2747
2748
2749 DATA a/0.99986d0, 6.690d-5, -8.486d-6, 1.518d-7, -6.9984d-9, -
2750 &3.6449d-10, -7.497d-12 /
2751
2752 ttemp=t-273d0
2753
2754 IF (ttemp .le. -40d0) THEN
2755 ttemp=-40d0
2756 END IF
2757
2758 watdensity=a(6)*ttemp
2759
2760 IF (t .GE. 240.0) THEN
2761 DO i=5,1, -1
2762 watdensity= (watdensity+a(i))*(ttemp)
2763 ENDDO
2764 watdensity=watdensity + a(0)
2765 ELSE
2766 watdensity=0.979d0
2767 END IF
2768
2769 watdensity=watdensity*1000d0
2770 watdensity_ice=watdensity
2771 return
2772 END function watdensity_ice
2773
2774
2775!*************************************************************
2776! Subroutine PROPERTIES. Set physical an thermodynamic
2777! properties at T and P for ice param
2778!************************************************************
2779
2780
2781 SUBROUTINE prop_ice(T, P, denice_ice,ddry_ice,
2782 & nin_ice,alfa_ice,beta_ice,shom_ice, koft_ice, dliq_ice,
2783 & g1_ice, g2_ice,gdoin_ice, z_ice,lambda_ice,
2784 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2785 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2786 & si0bc_ice,nbc_ice,D_preex, N_preex,one_over_tao_preex,
2787 & P_ice, T_ice,act_param,ndust_ice,vpresw_ice,vpresi_ice,
2788 & denair_ice)
2789
2790 real*8, intent(in) :: t, p,ddry_ice,nbc_ice,d_preex, n_preex
2791 real*8, intent(out) :: alfa_ice,beta_ice,shom_ice,
2792 & koft_ice, dliq_ice,g1_ice, g2_ice,gdoin_ice, z_ice,
2793 & lambda_ice,kdust_ice, kbc_ice, shdust_ice, shbc_ice,
2794 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
2795 & si0bc_ice,one_over_tao_preex,p_ice, t_ice
2796
2797 real*8 :: aux, aux1, aux2, sw, fice, mice, tc, hdust, hbc, b0,
2798 & b1, b2, b3, x, t0bc, t0dust, gam, gamma,nin_ice,
2799 & tr, vw, den_m, eact, toact, dact, acc, n1, siw, rgo,
2800 & ngo, mw_molec,dhs_ice,denwat_ice,denice_ice,vpresw_ice,
2801 & vpresi_ice,denair_ice,diff_ice,aircond_ice
2802
2803 real*8,dimension(:) :: ndust_ice
2804 integer, intent(in) :: act_param
2805
2806
2807 t_ice = min(max(t, tmin_ice), to_ice)
2808 p_ice = max(p, pmin_ice)
2809 dliq_ice = 1.0d-6
2810
2811
2812! rv_ice = rgas_ice/wmw_ice
2813 dhs_ice = dhsub_ice(t_ice)
2814 vpresw_ice = vpreswater_ice(t_ice)
2815 vpresi_ice = vpresice(t_ice)
2816 denice_ice = densityice(t_ice)
2817 denwat_ice = watdensity_ice(t_ice)
2818 denair_ice = p_ice*amw_ice/rgas_ice/t_ice
2819
2820
2821 diff_ice=(0.211d0*101325d0/p_ice)*((t_ice/273d0)**1.94d0)*1.0d-4
2822 aux1=1.0e-3*(4.39d0+0.071d0*t_ice)
2823
2824
2825 aux2=(2d0*aux1/(thaccom_ice*1.0d-6*denair_ice*cpa_ice)) *((58.0d-
2826 &3*pi_ice/(rgas_ice*t_ice))**0.5d0)
2827
2828 aircond_ice=aux1/(1.d0+aux2)
2829
2830
2831
2832
2833 aux1=grav_ice*dhs_ice/rv_ice/t_ice/t_ice/cpa_ice
2834 aux2=grav_ice*amw_ice/rgas_ice/t_ice
2835 alfa_ice=aux1-aux2
2836 beta_ice=amw_ice*p_ice/wmw_ice/vpresi_ice
2837 gamma=1.5d0*dhs_ice*dhs_ice/rv_ice/t_ice/t_ice/cpa_ice
2838
2839 beta_ice=beta_ice+gamma
2840
2841
2842
2843 shom_ice=2.349d0-(t_ice/259d0)
2844 sw=shom_ice*vpresi_ice/vpresw_ice
2845 shom_ice=shom_ice-1d0
2846 koft_ice=(0.0240d0*t_ice*t_ice)-(8.035d0*t_ice)+934.0d0
2847
2848
2849
2850
2851
2852 if (sw .lt. 0.99) then
2853 aux1=(1d0/(1d0-sw))-1.1764d0
2854 else
2855 aux1=(1d0/0.01)-1.1764d0
2856 end if
2857 dliq_ice=ddry_ice*0.9344d0*(aux1**0.333)
2858
2859
2860
2861 aux1=denice_ice*rv_ice*t_ice/vpresi_ice/diff_ice
2862 aux2=dhs_ice*denice_ice/aircond_ice/t_ice
2863 aux2=aux2*((dhs_ice/rv_ice/t_ice)-1.0d0)
2864 g1_ice=(aux1+aux2)/4.0d0
2865
2866
2867 one_over_tao_preex = beta_ice*denice_ice*pi_ice*0.5* d_preex*
2868 &n_preex/g1_ice/denair_ice
2869
2870
2871
2872 g2_ice=denice_ice*rv_ice*t_ice/2.0d0/vpresi_ice/depcoef_ice
2873 g2_ice=g2_ice*((2.0d0*pi_ice/rv_ice/t_ice)**0.5d0)
2874
2875 doin_ice=1.0d-6
2876 gdoin_ice=(g1_ice*0.5d0*doin_ice*doin_ice)+(g2_ice*doin_ice)
2877 z_ice=denice_ice*pi_ice/2.0d0/denair_ice
2878
2879 gam=g2_ice/g1_ice
2880 lambda_ice=1d0/sqrt(alfa_ice*g1_ice*gam*gam)
2881
2882
2883
2884
2885
2886 if (typeofspec_ice .lt. 0) then
2887 sh_ice=0.3d0
2888 nin_ice=(sum(ndust_ice)+nbc_ice)*0.05d0
2889 elseif (typeofspec_ice .eq. 3) then
2890
2891 shdust_ice = 0.2d0
2892 effdust_ice=0.6d0
2893 shbc_ice = 0.35d0
2894 effbc_ice=0.05d0
2895 mice = 0.96d0
2896 fice=0.25d0*((mice*mice*mice)-(3d0*mice)+2d0)
2897 kdust_ice=koft_ice*fice
2898 mice = 0.76d0
2899 fice=0.25d0*((mice*mice*mice)-(3d0*mice)+2d0)
2900 kbc_ice=koft_ice*fice
2901 elseif (typeofspec_ice .eq. 4) then
2902
2903 tc=t_ice-273.15d0
2904 hdust=0.15d0
2905 t0dust=-40d0
2906 b0=-1.0261d0; b1=3.1656d-3; b2=5.3938d-4; b3=8.2584d-6
2907 x=b0+(b1*tc)+(b2*tc*tc)+(b3*tc*tc*tc)
2908 si0dust_ice=1d0+(10d0**x)
2909 del1dust_ice=cubicint_ice(tc, t0dust, t0dust+5d0, 1d0, hdust)
2910 hbc=0d0
2911 t0bc=-50d0
2912 b0=0.5652d0; b1=1.085d-2; b2=-3.118d-5
2913 si0bc_ice=b0+(b1*t_ice)+(b2*t_ice*t_ice)-0.1d0
2914 del1bc_ice=cubicint_ice(tc, t0bc, t0bc+5d0, 1d0, hbc)
2915 end if
2916
2917 RETURN
2918
2919 END SUBROUTINE prop_ice
2920
2921!*************************************************************
2922! Subroutine gauspdf (normalized width of the updraft distribution).
2923!************************************************************
2924
2925 SUBROUTINE gausspdf(x, dp, sigmav_ice,miuv_ice,normv_ice)
2926
2927
2928 real*8, intent(in) :: x
2929 real*8 sigmav_ice,miuv_ice,normv_ice
2930 real*8, intent(out) :: dp
2931
2932 sigmav_ice =max(sigmav_ice, 0.01)
2933 normv_ice =max(normv_ice, 0.01)
2934
2935
2936 dp=exp(-0.5d0*(x-miuv_ice)*(x-miuv_ice)/sigmav_ice/sigmav_ice) /
2937 &sigmav_ice/sq2pi_par/(normv_ice + 0.001)
2938
2939
2940 RETURN
2941
2942
2943 END SUBROUTINE gausspdf
2944
2945
2946!*************************************************************
2947! Function cubicint_ice (cubic interpolation between y1 and y2 within a and b).
2948!************************************************************
2949
2950 real*8 function cubicint_ice(y, y1, y2, a, b)
2951
2952 real*8, intent(in) :: y, y1, y2, a, b
2953 real*8 :: a_, b_, a0, a1, a2, a3, d, aux
2954
2955 if (y .le. y1) then
2956 d=a
2957 goto 5065
2958 end if
2959
2960 if (y .ge. y2) then
2961 d=b
2962 goto 5065
2963 end if
2964
2965
2966 aux=y2-y1
2967 a_=6d0*(a-b)/(aux*aux*aux)
2968 b_=a+(a_*(y1*y1*y1)/6d0)-(a_*(y1*y1)*y2*0.5d0)
2969
2970 a0=b_
2971 a1=a_*y1*y2
2972 a2=-a_*(y1+y2)*0.5d0
2973 a3=a_/3d0
2974 d=a0+(a1*y)+(a2*y*y)+(a3*y*y*y)
2975
2976
2977 5065 cubicint_ice=d
2978
2979
2980 end function cubicint_ice
2981
2982!*************************************************************
2983! Function dcubicint_ice (used in the PDA08 spectrum).
2984!************************************************************
2985
2986 real*8 function dcubicint_ice(y, y1, y2, a, b)
2987
2988 real*8, intent(in) :: y, y1, y2, a, b
2989 real*8 :: a_, a0, a1, a2, a3, d, aux
2990
2991 if (y .le. y1) then
2992 d=0
2993 goto 5065
2994 end if
2995
2996 if (y .ge. y2) then
2997 d=0
2998 goto 5065
2999 end if
3000
3001
3002 aux=y2-y1
3003 a_=6d0*(a-b)/(aux*aux*aux)
3004
3005 a1=a_*y1*y2
3006 a2=-a_*(y1+y2)*0.5d0
3007 a3=a_/3d0
3008 d=(a1)+(2d0*a2*y)+(3d0*a3*y*y)
3009
3010
3011 5065 dcubicint_ice=d
3012
3013
3014 end function dcubicint_ice
3015
3016!*************************************************************
3017! Function PDG07 (simplified ice nucleation
3018! spectra according to Phillips et. al. 2007).
3019! si is supersaturation wrt ice and T is in K
3020!************************************************************
3021
3022 real*8 function pdg07_ice(si, T)
3023
3024 real*8, intent(in) :: si, t
3025 real*8 :: n
3026
3027 if (t .le. 243d0)then
3028 n=1000d0*exp(-0.388d0)*(exp(3.88d0*si)-1d0)
3029 else
3030 n=60d0*exp(-0.639d0)*(exp(12.96d0*si)-1d0)
3031 end if
3032
3033 pdg07_ice=n
3034
3035 end function pdg07_ice
3036
3037
3038
3039
3040!*************************************************************
3041! Subroutine INSPEC_ice
3042! Provides the Ice Nuclei concentration (m-3)
3043! and the chracteristic freezing threeshold, DSh (Barahona & Nenes 2009), at given
3044! si and T. The variable typeofspec_ice (integer) has the values
3045! 1 Meyers et. al. 1992
3046! 2 Phillips et. al. 2007
3047! 3 Barahona 2011
3048! 4 Phillips et. al. 2008 (simplifed)
3049! 5 Phillips et. al. 2013 (simplifed)
3050! si is supersaturation wrt ice and T is in K
3051
3052! Written by Donifan Barahona
3053! donifan.o.barahona@nasa.gov
3054
3055!************************************************************
3056
3057 subroutine inspec_ice(six, N, Dsh,np_ice,norg_ice, sigorg_ice,
3058 & dorg_ice, dbc_ice,sigbc_ice,
3059 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
3060 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
3061 & si0bc_ice,nbc_ice,Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice,
3062 & ndust_ice, sigdust_ice,ddust_ice,nbindust_ice,
3063 & vpresw_ice,vpresi_ice)
3064
3065 real*8, intent(in) :: six,np_ice,norg_ice, sigorg_ice,
3066 & dorg_ice, dbc_ice,sigbc_ice,
3067 & kdust_ice, kbc_ice, shdust_ice, shbc_ice,
3068 & effdust_ice, effbc_ice, del1dust_ice, si0dust_ice, del1bc_ice,
3069 & si0bc_ice,nbc_ice,p_ice, t_ice
3070 real*8, intent(out) :: n, dsh,nhet_dep,nhet_dhf,fdust_dep
3071 real*8 :: nd, nbc, aux, si_, sw, del0, ddel0, fc, delw0, ddelw0,
3072 & sw0, hdust, hbc, nbase, dnd, dnbc, dnbase, dh, dfc, ndaux,
3073 & dndaux, dnorg, norg, ndustaux, frac, aux2, dx2, fdep, ndep, ndhf,
3074 & dndep, dndhf, si, dfrac, ncdep_, ncdhf_ , fglassy, nglassy,
3075 & dnglassy, siw, d_grid_bio, n_grid_bio,vpresw_ice,vpresi_ice
3076
3077 real*8, dimension(3) :: sig_array, the_array, frac_array
3078 real*8, dimension(1:nbindust_ice) :: ndust_ice, sigdust_ice,
3079 & ddust_ice
3080
3081 real :: n_iw, DSh_s , nbc_s, dbc_s, Asolo
3082 real, dimension (nbindust_ice) :: ndust_s, ddust_s
3083 integer :: index, kindex, mode,nbindust_ice
3084
3085 si=six
3086 si_=si+1d0
3087 sw=si_*vpresi_ice/vpresw_ice
3088 nd = zero_par
3089 nbc = zero_par
3090 norg = zero_par
3091 nglassy = zero_par
3092
3093 dnd = zero_par
3094 dnbc = zero_par
3095 dnorg = zero_par
3096 dnglassy = zero_par
3097
3098 if ((six .lt. 0.02) .or. (t_ice .gt. 270.0)) then
3099 n=0d0
3100 dsh=si
3101 return
3102 end if
3103
3104 if (sw .ge. 1.0) then
3105 sw=1.0
3106 si_=vpresw_ice/vpresi_ice
3107 si=si_-1.0
3108 end if
3109 siw=vpresw_ice/vpresi_ice
3110
3111 sig_array = 0.0
3112 the_array = 0.0
3113 frac_array = 0.0
3114
3115
3116 select case (typeofspec_ice)
3117
3118 case(1)
3119 n=1000d0*exp(-0.639d0)*(exp(12.96d0*si)-1d0)
3120 dsh=1d0/12.96d0
3121
3122 case(2)
3123 n=pdg07_ice(si, t_ice)
3124 if (t_ice .le. 243d0)then
3125 dsh=1d0/3.88d0
3126 else
3127 dsh=1d0/12.96d0
3128 end if
3129
3130 case(3)
3131
3132 ndustaux=0.0d0
3133 DO index=1, nbindust_ice
3134
3135 ndustaux=ndustaux+ndust_ice(index)
3136 end do
3137
3138
3139 if (si .le.shdust_ice) then
3140 nd=(si/shdust_ice)*ndustaux*effdust_ice* exp(-kdust_ice*
3141 &(shdust_ice-si))
3142 dnd=nd*((1d0/si)+kdust_ice)
3143
3144 else
3145 nd=ndustaux*effdust_ice
3146 dnd=0d0
3147 end if
3148
3149
3150 if (si .le.shbc_ice) then
3151 nbc=(si/shbc_ice)*nbc_ice*effbc_ice* exp(-kbc_ice*(shbc_ice-si))
3152 dnbc=nbc*((1d0/si)+kbc_ice)
3153 else
3154 nbc=nbc_ice*effbc_ice
3155 dnbc=0d0
3156 end if
3157
3158 n=nd+nbc
3159 if (((dnd+dnbc) .gt. 0d0) .and. (n .gt. 0.0)) then
3160 dsh=n/(dnd+dnbc)
3161 else
3162 dsh=0.0
3163 end if
3164
3165 case(4)
3166
3167
3168
3169 sw0=0.97d0
3170 delw0=cubicint_ice(sw, sw0, 1d0, 0d0, 1d0)
3171 ddelw0=dcubicint_ice(sw, sw0, 1d0, 0d0, 1d0)
3172
3173 nbase=pdg07_ice(si, t_ice)
3174
3175
3176 if (t_ice .le. 243d0)then
3177 dnbase=3.88d0*nbase
3178 else
3179 dnbase=12.96d0*nbase
3180 end if
3181
3182
3183 del0=cubicint_ice(si_, si0dust_ice, si0dust_ice+0.1d0, 0d0, 1d0)
3184 ddel0=dcubicint_ice(si_, si0dust_ice, si0dust_ice+0.1d0, 0d0,
3185 & 1d0)
3186
3187 fc=0.5d0*del1dust_ice*del0
3188 dfc=0.5d0*del1dust_ice*ddel0
3189
3190 hdust=fc+((1d0-fc)*delw0)
3191 dh=(dfc*(1d0-delw0))+(ddelw0*(1d0-fc))
3192
3193 if (hdust .gt. 1d0) then
3194 hdust=1d0
3195 dh=0d0
3196 end if
3197
3198
3199 aux=(2d0/3d0)*hdust*(nbase/0.76d0)*pi_ice/5.0d-7/4d0
3200 aux2=(2d0/3d0)*pi_ice/0.76d0/5.0d-7/4d0
3201
3202 nd=0d0
3203 dnd=0d0
3204
3205 DO index =1, nbindust_ice
3206
3207
3208 dx2=ddust_ice(index)*ddust_ice(index)*ddust_ice(index)*0.52*
3209 &acorr_dust
3210
3211 frac=0.5d0*(1d0-erfapp(-log(ddust_ice(index)/0.1e-6) /
3212 &sigdust_ice(index)/sq2_par))
3213
3214 ndaux=frac*ndust_ice(index)*(1d0-exp(-aux*dx2))
3215
3216 nd=nd+ndaux
3217 ndaux=(frac*ndust_ice(index)-ndaux)
3218 dndaux=ndaux*((dh*nbase)+(hdust*dnbase))*aux2*dx2
3219
3220 dnd=dnd+dndaux
3221
3222 END DO
3223
3224
3225
3226
3227 del0=cubicint_ice(si_, si0bc_ice, si0bc_ice+0.1d0, 0d0, 1d0)
3228 ddel0=dcubicint_ice(si_, si0bc_ice, si0bc_ice+0.1d0, 0d0, 1d0)
3229
3230 fc=0.5d0*del1bc_ice*del0
3231 hbc=fc+((1d0-fc)*delw0)
3232 dfc=0.5d0*del1bc_ice*ddel0
3233 dh=(dfc*(1d0-delw0))+(ddelw0*(1d0-fc))
3234
3235
3236 if (hbc .gt. 1d0) then
3237 hbc=1d0
3238 dh=0d0
3239 end if
3240
3241 frac=0.5d0*(1d0 -erfapp(-log(dbc_ice/0.1e-6) /sigbc_ice/
3242 &sq2_par))
3243
3244
3245 dx2=dbc_ice*dbc_ice*dbc_ice*0.52*acorr_bc
3246
3247 aux=((1d0/3d0)-0.06d0)*hbc*(nbase/0.76d0)*pi_ice/2.7d-7
3248 aux2=((1d0/3d0)-0.06d0)*pi_ice/0.76d0/2.7d-7
3249
3250 nbc=nbc_ice*frac*(1d0-exp(-aux*dx2))
3251 dnbc=(nbc_ice*frac-nbc)*((dh*nbase)+(hbc*dnbase))*aux2*dx2
3252
3253
3254
3255
3256 frac=0.5d0*(1d0-erfapp(-log(dorg_ice/0.1e-6) /sigorg_ice/
3257 &sq2_par))
3258
3259 dx2=dorg_ice*dorg_ice
3260
3261 aux=0.06d0*hbc*(nbase/0.76d0)*pi_ice/9.1d-7
3262 aux2=0.06d0*pi_ice/0.76d0/9.1d-7
3263
3264
3265 norg=norg_ice*frac*(1d0-exp(-aux*dx2))
3266 dnorg=(norg_ice*frac-norg)*((dh*nbase)+(hbc*dnbase))*aux2*dx2
3267
3268
3269
3270
3271 n=nd+nbc+norg
3272
3273 fdust_dep = nd
3274 nhet_dep = n
3275
3276
3277 if (.false.) then
3278
3279 if (t_ice .lt. 210.0) then
3280 nglassy= min(0.01 +0.0045*(210.0 -t_ice), 0.1)
3281 fglassy= 7.7211*1e-3*si_ - 9.2688*1e-3
3282 fglassy=min(fglassy, 3.3587e-3)
3283 fglassy=max(fglassy, 0.0)
3284 else
3285 nglassy = 0.0
3286 fglassy = 0.0
3287 end if
3288 nglassy = np_ice*nglassy*fglassy
3289 dnglassy =nglassy*7.7211*1e-3
3290 n=n+nglassy
3291
3292 end if
3293
3294
3295
3296 if ((dnd+dnbc+dnorg+dnglassy) .gt. 0d0) then
3297 dsh=n/(dnd+dnbc+dnorg+dnglassy)
3298 else
3299 dsh=0.0
3300 end if
3301
3302
3303 case (5)
3304
3305
3306
3307 d_grid_bio =dorg_ice
3308
3309
3310
3311 n_grid_bio = 0.0
3312
3313 if (is_gocart) then
3314 ndust_s = sngl(frac_dust*ndust_ice)
3315 nbc_s = sngl(frac_bc*nbc_ice)
3316 asolo = sngl(0.25d0*frac_org*norg_ice*pi_par*dorg_ice*dorg_ice)
3317 else
3318
3319
3320 DO index =1, nbindust_ice
3321 frac=0.5d0*(1d0-erfapp(log(0.1e-6/ddust_ice(index)) /
3322 &sigdust_ice(index)/sq2_par))
3323 ndust_s(index) = sngl(frac*ndust_ice(index))
3324 end do
3325
3326 frac=0.5d0*(1d0-erfapp(log(0.1e-6/dbc_ice) /sigbc_ice/sq2_par))
3327
3328 nbc_s = sngl(frac*nbc_ice)
3329
3330 frac=0.5d0*(1d0-erfapp(log(0.1e-6/dorg_ice) /sigorg_ice/
3331 &sq2_par))*0.25d0
3332
3333 asolo = sngl(frac*norg_ice*pi_par*dorg_ice*dorg_ice)
3334 end if
3335
3336 ddust_s=sngl(ddust_ice)
3337 dbc_s = sngl(dbc_ice*1.0)
3338
3339 call empirical_param_phillips(sngl(si_), sngl(siw), sngl(sw), (/
3340 &ddust_s/), (/ndust_s/), 5, (/dbc_s/), (/nbc_s/), 1, (/
3341 &sngl(d_grid_bio)/), (/sngl(n_grid_bio)/), 1, asolo, n_iw, dsh_s
3342 & ,nhet_dep,nhet_dhf,fdust_dep,p_ice, t_ice)
3343
3344 n=dble(n_iw)
3345 dsh=dble(dsh_s)
3346
3347 case default
3348 n=zero_par
3349 dsh=0.0
3350 end select
3351
3352 if (dsh .ge. si) then
3353 dsh=si
3354 end if
3355
3356 if (t_ice .gt. 255.0) then
3357 n=zero_par
3358 dsh=zero_par
3359 end if
3360
3361 end subroutine inspec_ice
3362
3363
3364!*************************************************************
3365! Subroutine INimmersion
3366! Provides the Immersion IN (for activated droplets) concentration at given T(K) according to Barahona et al. GMD, 2014.
3367! Written by Donifan Barahona
3368! donifan.o.barahona@nasa.gov
3369!========================================
3370
3371 subroutine inimmersion(INconc, dINconc, wparcel,dbc_ice,sigbc_ice
3372 & ,nbc_ice,fdust_imm,fdrop_dust,fdrop_bc,ndust_ice, sigdust_ice,
3373 & ddust_ice,nbindust_ice,vpresw_ice,vpresi_ice,T_ice)
3374
3375 real*8, intent (in) :: wparcel,dbc_ice,sigbc_ice,nbc_ice
3376 & ,fdrop_dust,fdrop_bc
3377 real*8, intent (out) :: inconc, dinconc,fdust_imm
3378 real*8 :: nd, nd_unc, nd_coa, nbc, ahet, frac, dfrac, naux,dnaux,
3379 & tx, nssoot, nsdust, ninbkg, sx , dnsd, dnss, dnd, dnbc, coolr,
3380 & min_ns_dust, min_ns_soot, dninbkg,vpresw_ice,vpresi_ice,t_ice
3381
3382 real*8, dimension(3) :: sig_array, the_array, frac_array
3383 real*8, dimension(:) :: ndust_ice, sigdust_ice,ddust_ice
3384 integer :: index, kindex, nbindust_ice
3385
3386 if (t_ice .lt. thom) then
3387 inconc = zero_par
3388 dinconc = zero_par
3389 return
3390 end if
3391
3392 if (t_ice .ge. to_ice) then
3393 inconc = zero_par
3394 dinconc = zero_par
3395 return
3396 end if
3397
3398 tx=t_ice -273.16
3399 coolr=5.0e-3*wparcel
3400 min_ns_dust= 3.75e6
3401 min_ns_soot= 3.75e9
3402
3403 sx=(vpresw_ice/vpresi_ice)-1.0
3404
3405
3406 ninbkg = 0.0
3407 dninbkg=0.0
3408
3409
3410
3411
3412
3413 nsdust= max(exp(-0.517*tx + 8.934)-min_ns_dust, 0.0)
3414 dnsd = max(0.517*nsdust, 0.0)
3415
3416 nssoot= max(1.0e4*exp(-0.0101*tx*tx - 0.8525*tx + 0.7667)-
3417 &min_ns_soot, 0.0)
3418 dnss = max(-(-2.0*0.0101*tx -0.8525)*nssoot, 0.0)
3419
3420 naux=zero_par
3421 dnaux=zero_par
3422
3423
3424 if (is_gocart) then
3425
3426 DO index =1,nbindust_ice
3427
3428 ahet= ahet_dust(index)
3429 naux=(1.0-exp(-nsdust*ahet))*ndust_ice(index)+naux
3430 dnaux = exp(-nsdust*ahet)*ndust_ice(index)*dnsd*coolr*ahet+dnaux
3431 END DO
3432 ahet = ahet_bc
3433 else
3434
3435 DO index =1,nbindust_ice
3436
3437 ahet= ddust_ice(index)*ddust_ice(index)*ddust_ice(index)*0.52*
3438 &acorr_dust* exp(4.5*sigdust_ice(index)*sigdust_ice(index))
3439 naux=(1.0-exp(-nsdust*ahet))*ndust_ice(index)+naux
3440 dnaux = exp(-nsdust*ahet)*ndust_ice(index)*dnsd*coolr*ahet+dnaux
3441 END DO
3442
3443 ahet=dbc_ice*dbc_ice*dbc_ice*0.52*acorr_bc* exp(4.5*sigbc_ice*
3444 &sigbc_ice)
3445 end if
3446
3447
3448 nd=naux*fdrop_dust
3449 dnd=dnaux*fdrop_dust
3450
3451 nbc=nbc_ice*(1.0-exp(-nssoot*ahet))*fdrop_bc
3452 dnbc= nbc_ice*exp(-nssoot*ahet)*fdrop_bc*dnss*coolr*ahet
3453
3454
3455
3456
3457 inconc=nbc+ nd + ninbkg
3458 dinconc=dnbc+dnd + dninbkg
3459 fdust_imm = nd
3460
3461
3462 end subroutine inimmersion
3463
3464!
3465!=======================================================================================
3466!=======================================================================================
3467!=======================================================================================
3468!!====================================================================================
3469! EMPIRICAL PARAMETERISATION (Phillips et al. 2013, JAS)
3470! Code contributed by Vaughan Phillips, University of Leeds
3471! Implementation: Donifan Barahona donifan.o.barahona@nasa.gov
3472!====================================================================================
3473
3474
3475 SUBROUTINE empirical_param_phillips(SI, SIW, SW, D_grid_dust,
3476 & n_grid_dust, ijstop_dust, D_grid_soot, n_grid_soot, ijstop_soot,
3477 & D_grid_bio, n_grid_bio, ijstop_bio, A_solo, n_iw, DSH,
3478 & Nhet_dep,Nhet_dhf,fdust_dep,P_ice, T_ice)
3479 implicit none
3480 real, intent(IN):: SI, SIW, SW, A_solo
3481 real*8, intent(IN):: p_ice, t_ice
3482 real, dimension(:), intent(IN):: D_grid_dust, n_grid_dust,
3483 & d_grid_soot, n_grid_soot, d_grid_bio, n_grid_bio
3484 integer, intent(IN):: ijstop_dust, ijstop_soot, ijstop_bio
3485! integer ijstop_dust, ijstop_soot, ijstop_bio
3486
3487 real :: nin_a_nuc_dust, nin_a_nuc_soot, nin_a_nuc_bio,
3488 & nin_a_nuc_solo, num_ic_dust_imm, num_ic_soot_imm, num_ic_bio_imm,
3489 & num_ic_solo_imm
3490
3491 real, intent (inout) :: DSH, n_iw
3492 real*8, intent (out) :: nhet_dep,nhet_dhf,fdust_dep
3493
3494 real :: dn_in_dust, dn_in_soot, dn_in_bio, dn_in_solo, dNall,
3495 & dnaux, naux, ss_w, dh_frac_dust, dh_frac_soot, dh_frac_solo, aux,
3496 & dfdep, temperature_k, p_sat, ahet_aux
3497
3498
3499
3500 REAL :: RHO_CFDC, BASE_DUST_OMEGA, BASE_SOOT_PHILIC_OMEGA,
3501 & base_bio_omega, alpha_dust, alpha_soot, alpha_bio,
3502 & fraction_depnucl_warm_dust, pie, base_solo_omega,
3503 & temp_max_dust_degc, temp_max_soot_degc, temp_max_bio_degc,
3504 & glass_frac
3505 parameter( base_dust_omega = 2.0e-6, base_soot_philic_omega = 1.e-
3506 &7, base_bio_omega = 0.89e-6, base_solo_omega = 5.6e-5,
3507 & glass_frac = 0.5, alpha_dust = 2./3., alpha_soot = 1./3. - 0.03,
3508 & alpha_bio = 0.03, rho_cfdc = 50000./(287.*228.15),
3509 & fraction_depnucl_warm_dust = 0.15, pie = 3.1415926,
3510 & temp_max_dust_degc = -10., temp_max_soot_degc = -15.,
3511 & temp_max_bio_degc = -2.)
3512
3513 real :: FAC_CORRECT_RH = 2., rho_aida
3514 real:: H_frac_dust, n_in, n_in_dust, n_in_ultra, n_in_dust_ultra,
3515 & cihenc_dust, esw, esi, ss_i, n_in_soot_ultra, h_frac_soot,
3516 & h_frac_bio, n_in_soot, n_in_bio, n_in_bio_ultra, cihenc_soot,
3517 & cihenc_bio, delta_si, delta_t, delta_sw, n_in_max, ss_iw, rho
3518
3519 real :: H_frac_solO, RHI, n_in_solO, n_in_solO_star, CIHENC_solO,
3520 & psi_solo
3521
3522 real :: mu, S_i_0, RH_crit, S_i_w_warm, S_i_w_cold, S_i_w,
3523 & tc_hm_degc
3524 real :: S_w_0, dep_frac, n_in_hat, n_in_tilde
3525 real*4 :: dh1smooth
3526 real :: EPS = 0.622
3527 integer :: ij
3528!intrinsic :: exp, DEXP, SIZE, DBLE
3529
3530
3531
3532!print *, SIZE(n_grid_dust(:))
3533 if(ijstop_dust .ne. SIZE(n_grid_dust)) stop 6366
3534 if(ijstop_soot .ne. SIZE(n_grid_soot)) stop 6366
3535 if(ijstop_bio .ne. SIZE(n_grid_bio)) stop 6366
3536! ijstop_dust = SIZE(n_grid_dust)
3537! ijstop_soot = SIZE(n_grid_soot)
3538! ijstop_bio = SIZE(n_grid_bio)
3539
3540
3541
3542!Naux=12.96 !default
3543 naux =0.0
3544! Nall =dNaux
3545! Sh=0.0
3546 n_iw=0.0
3547 nin_a_nuc_dust=0.0; nin_a_nuc_soot=0.0; nin_a_nuc_bio=
3548 &0.0; nin_a_nuc_solo=0.0
3549 num_ic_dust_imm=0.0; num_ic_soot_imm=0.0; num_ic_bio_imm=
3550 &0.0; num_ic_solo_imm=0.0
3551! AnningC uncom following four lines
3552 n_in_dust=0.0; dn_in_dust=0.0;n_in_soot=0.0
3553 dn_in_soot=0.0; dn_in_bio=0.0; n_in_bio=0.0;
3554 dn_in_solo=0.0;n_in_solo=0.0;n_in_dust_ultra=0.0
3555 n_in_soot_ultra=0.0;n_in_bio_ultra=0.0;
3556
3557
3558 h_frac_dust = 0.0
3559 h_frac_soot = 0.0
3560 h_frac_solo = 0.0
3561! H1smooth=0.0
3562 aux=0.0
3563
3564 temperature_k=sngl(t_ice)
3565 p_sat =sngl(p_ice)
3566
3567!A_solo = 1e-7 !m2 kg-1
3568
3569
3570!====================================================================================
3571! COMPUTATION BLOCK
3572!
3573!====================================================================================
3574!
3575 rho_aida = 90000./(287.*205.)
3576
3577 rho = p_sat/(287.*temperature_k)
3578
3579 psi_solo = a_solo/base_solo_omega
3580 ss_i = min(max(si-1.0, 0.0), 1.0)
3581 ss_w = min(max(sw-1.0, -1.0), 1.0)
3582 ss_iw = min(max(siw - 1.0, 0.0), 1.0)
3583
3584
3585 if(ss_i > 0.0) then
3586 if(temperature_k < 273.15 .and. temperature_k > 273.15 -
3587 & 90. ) then
3588
3589
3590 if(ss_w > 0.) then
3591 ss_i = ss_iw
3592 ss_w = 0.0
3593 end if
3594
3595
3596! S_i_zero = 1.15 !this is taken care of
3597
3598
3599 delta_si = h_1_smooth(ss_i + 1, 1.1, 1.2, 0.0, 1.,dh1smooth);
3600 delta_t = h_1_smooth(-(temperature_k-273.15), 35., 40.,
3601 & fraction_depnucl_warm_dust, 1.,dh1smooth);
3602 delta_sw = h_1_smooth(ss_w + 1.0, 0.97, 1., 0., 1.,dh1smooth);
3603
3604 tc_hm_degc = temperature_k - 273.15
3605
3606
3607 s_i_0 = 1. + 10.**(8.2584e-6*tc_hm_degc*tc_hm_degc*tc_hm_degc +
3608 & 5.3938e-4*tc_hm_degc*tc_hm_degc + 3.1656e-3*tc_hm_degc - 1.0261)
3609
3610
3611 s_w_0 = 0.97
3612
3613 aux =h_1_smooth(-(temperature_k-273.15), 35., 40.,
3614 & fraction_depnucl_warm_dust, 1.,dh1smooth)/fac_correct_rh
3615
3616 dep_frac = h_1_smooth(ss_i + 1, s_i_0, s_i_0 + 0.1, 0.,1.,
3617 & dh1smooth)* aux
3618 dfdep=dh1smooth*aux
3619
3620 aux= h_1_smooth(ss_w + 1.0, s_w_0, 1., 0.,1.,dh1smooth)
3621
3622 h_frac_dust = dep_frac + (1. - dep_frac)*aux
3623
3624 dh_frac_dust = dfdep + (siw*(1. - dep_frac)*dh1smooth)- aux*
3625 &dfdep
3626
3627 if(h_frac_dust > 1.) h_frac_dust = 1.
3628
3629 if ((h_frac_dust .gt. 1.0e-6) .and. (h_frac_dust .lt. 1.)) then
3630 dh_frac_dust = dh_frac_dust/h_frac_dust
3631 else
3632 dh_frac_dust =0.0
3633 end if
3634
3635
3636 s_i_0 = 1.2
3637
3638 aux =h_1_smooth(-(temperature_k-273.15), 65., 75., 0.,1.,
3639 & dh1smooth)
3640 dep_frac = h_1_smooth(ss_i + 1, s_i_0, s_i_0+0.1, 0.,1.,
3641 & dh1smooth)*aux
3642 h_frac_solo = dep_frac
3643 dh_frac_solo=0.0
3644 if ((h_frac_solo .gt. 1.0e-6) .and. (h_frac_solo .lt. 1.)) then
3645 dh_frac_solo = dh1smooth/h_frac_solo
3646 end if
3647 if(h_frac_solo > 1.) h_frac_solo = 1.
3648
3649
3650 s_w_0 = 0.97
3651
3652 s_i_0 = 1.3
3653!
3654 aux = h_1_smooth(-(temperature_k-273.15), 40., 50., 0.,1.,
3655 & dh1smooth) / fac_correct_rh
3656 dep_frac = h_1_smooth(ss_i + 1, s_i_0, s_i_0+0.1, 0.,1.,
3657 & dh1smooth)* aux
3658
3659 dfdep= dh1smooth*aux
3660
3661 aux = h_1_smooth(ss_w + 1.0, s_w_0, 1., 0.,1.,dh1smooth)
3662 h_frac_soot = dep_frac + (1. - dep_frac)*aux
3663 if(h_frac_soot > 1.) h_frac_soot = 1.
3664
3665 dh_frac_soot = dfdep + (siw*(1. - dep_frac)*dh1smooth)- aux*
3666 &dfdep
3667 if ((h_frac_soot .gt. 1.0e-6) .and. (h_frac_soot .lt. 1.)) then
3668 dh_frac_soot = dh_frac_soot/h_frac_soot
3669 else
3670 dh_frac_soot =0.0
3671 end if
3672
3673
3674 h_frac_bio = h_frac_soot
3675
3676 if(temperature_k < 273.15 .and. temperature_k >= 273.15 -
3677 & 35.) then
3678 n_in = 1.e3* (exp(12.96*ss_i - 0.639)/rho_cfdc) *0.0587*
3679 &fac_correct_rh
3680 if( temperature_k > 273.15 -5. .and. temperature_k < 273.15 -
3681 & 2. ) then
3682 n_in = n_in*h_1_smooth(-(temperature_k-273.15), 2., 5., 0., 1.,
3683 & dh1smooth)
3684 endif
3685 if(temperature_k >= 273.15 - 2. ) n_in = 0.
3686
3687
3688 if(temperature_k < 273.15 -25. ) then
3689 n_in_tilde = 1000.*(exp(0.1296*(ss_i*100.-10.))**0.3)*
3690 &fac_correct_rh/rho_cfdc
3691 n_in_hat = n_in
3692
3693 if(temperature_k >= 273.15 - 30.) n_in_max = 1.e3* (exp(12.96*
3694 &ss_iw - 0.639)/rho_cfdc) *0.0587*fac_correct_rh
3695 if(temperature_k < 273.15 - 30.) n_in_max = 1000.*(exp(0.1296*
3696 &(ss_iw*100.-10.))**0.3)*fac_correct_rh/rho_cfdc
3697
3698 if(n_in_hat > n_in_max) n_in_hat = n_in_max
3699 if(n_in_tilde > n_in_max) n_in_tilde = n_in_max
3700
3701
3702
3703 n_in = n_in_hat * ((n_in_tilde/n_in_hat)**(h_1_smooth(-
3704 &(temperature_k-273.15), 25., 35., 0., 1.,dh1smooth)))
3705
3706
3707 if(n_in > n_in_max) n_in = n_in_max
3708
3709 endif
3710 n_in_dust = 0.
3711 dn_in_dust = 0.
3712
3713
3714 if(temperature_k < 273.15 - 30.) then
3715 dnaux = 3.88
3716 else
3717 dnaux = 12.96
3718 end if
3719
3720
3721 naux=0.0
3722 do ij = 1, ijstop_dust
3723 if (is_gocart) then
3724 mu = n_in*alpha_dust*h_frac_dust*ahet_dust(ij)/base_dust_omega
3725 else
3726 ahet_aux = pie*d_grid_dust(ij)*d_grid_dust(ij)*d_grid_dust(ij)*
3727 &4.73*acorr_dust/6.0
3728 mu = n_in*alpha_dust*h_frac_dust*ahet_aux/base_dust_omega
3729 end if
3730 naux = (1. - exp(-mu))*n_grid_dust(ij)
3731 n_in_dust = n_in_dust + naux
3732 dn_in_dust = max(mu*(n_grid_dust(ij)-naux)*(dnaux +
3733 & dh_frac_dust), 0.0) + dn_in_dust
3734
3735 enddo
3736
3737 if( temperature_k > 273.15 +temp_max_dust_degc -
3738 & 20. .and. temperature_k < 273.15 + temp_max_dust_degc) then
3739 n_in_dust = n_in_dust*h_1_smooth(-(temperature_k-273.15),-
3740 &temp_max_dust_degc,-temp_max_dust_degc+20., 0., 1.,dh1smooth)
3741 endif
3742 if(temperature_k >= 273.15 + temp_max_dust_degc) n_in_dust = 0.
3743
3744
3745 n_in_soot = 0.
3746 dn_in_soot = 0.
3747 do ij = 1, ijstop_soot
3748
3749 if (is_gocart) then
3750 mu = n_in*alpha_soot*h_frac_soot*ahet_bc/base_soot_philic_omega
3751 else
3752 ahet_aux = pie*d_grid_soot(ij)*d_grid_soot(ij)*d_grid_soot(ij)*
3753 &16.4*acorr_bc/6.0
3754 mu = n_in*alpha_soot*h_frac_soot*ahet_aux/base_soot_philic_omega
3755 end if
3756 naux = (1. - exp(-mu))*n_grid_soot(ij)
3757 n_in_soot = n_in_soot + naux
3758 dn_in_soot = max(mu*(n_grid_soot(ij)-naux)*(dnaux+dh_frac_soot),
3759 & 0.0) + dn_in_soot
3760
3761 enddo
3762
3763 if( temperature_k > 273.15 + temp_max_soot_degc -
3764 & 10. .and. temperature_k < 273.15 + temp_max_soot_degc) then
3765 n_in_soot = n_in_soot*h_1_smooth(-(temperature_k-273.15),-
3766 &temp_max_soot_degc,-temp_max_soot_degc+10., 0., 1.,dh1smooth)
3767
3768 endif
3769 if(temperature_k >= 273.15 + temp_max_soot_degc) n_in_soot = 0.
3770
3771 n_in_bio = 0.
3772 dn_in_bio = 0.
3773
3774
3775 do ij = 1, ijstop_bio
3776 mu = n_in*alpha_bio*h_frac_bio*pie*(d_grid_bio(ij)**2.) /
3777 &base_bio_omega
3778
3779 mu = n_in*alpha_bio*h_frac_bio
3780 naux = (1. - exp(-mu))*n_grid_bio(ij)
3781
3782
3783
3784
3785 n_in_bio = n_in_bio + naux
3786 dn_in_bio = max(mu*(n_grid_bio(ij)-naux)*dnaux, 0.0) + dn_in_bio
3787
3788
3789 enddo
3790
3791
3792 if( temperature_k > 273.15 + temp_max_bio_degc -
3793 & 3. .and. temperature_k < 273.15 + temp_max_bio_degc) then
3794 n_in_bio = n_in_bio*h_1_smooth(-(temperature_k-273.15),-
3795 &temp_max_bio_degc,-temp_max_bio_degc+3., 0., 1.,dh1smooth)
3796
3797 endif
3798 if(temperature_k >= 273.15 + temp_max_bio_degc ) n_in_bio = 0.
3799
3800
3801
3802 else
3803 n_in = 0.; n_in_ultra = 0.; n_in_dust = 0.; n_in_soot =
3804 & 0.; n_in_bio = 0.;
3805 endif
3806
3807 if(temperature_k < 273.15 - 35.) then
3808 n_in_ultra = 1000.*(exp(0.1296*(ss_i*100.-10.))**0.3)*
3809 &fac_correct_rh/rho_cfdc
3810 dnaux = 3.88
3811 naux=0.0
3812
3813
3814 rhi = (ss_i+1.)*100.
3815 if(rhi < 0.) rhi = 0.
3816 n_in_solo_star = 1000.e6*(7.7211e-5 * rhi - 9.2688e-3)/rho_aida
3817
3818
3819 n_in_dust_ultra = 0.;
3820 dn_in_dust = 0.0
3821 do ij = 1, ijstop_dust
3822
3823 if (is_gocart) then
3824 mu = n_in_ultra*alpha_dust*h_frac_dust*ahet_dust(ij)/
3825 &base_dust_omega
3826 else
3827 ahet_aux = pie*d_grid_dust(ij)*d_grid_dust(ij)*d_grid_dust(ij)*
3828 &4.73*acorr_dust/6.0
3829 mu = n_in_ultra*alpha_dust*h_frac_dust*ahet_aux/base_dust_omega
3830 end if
3831
3832
3833
3834 naux = (1. - exp(-mu))*n_grid_dust(ij)
3835 n_in_dust_ultra = n_in_dust_ultra + naux
3836 dn_in_dust = max(mu*(n_grid_dust(ij)-naux)*(dnaux +dh_frac_dust),
3837 & 0.0) + dn_in_dust
3838
3839
3840
3841 enddo
3842
3843
3844
3845 n_in_soot_ultra = 0.0
3846 dn_in_soot = 0.0
3847 do ij = 1, ijstop_soot
3848
3849 if (is_gocart) then
3850 mu = n_in_ultra*alpha_soot*h_frac_soot*ahet_bc/
3851 &base_soot_philic_omega
3852 else
3853 ahet_aux = pie*d_grid_soot(ij)*d_grid_soot(ij)*d_grid_soot(ij)*
3854 &16.4*acorr_bc/6.0
3855 mu = n_in_ultra*alpha_soot*h_frac_soot*ahet_aux/
3856 &base_soot_philic_omega
3857 end if
3858
3859
3860
3861 naux = (1. - exp(-mu))*n_grid_soot(ij)
3862 n_in_soot_ultra = n_in_soot_ultra + naux
3863 dn_in_soot = max(mu*(n_grid_soot(ij)-naux)*(dnaux +dh_frac_soot),
3864 & 0.0) + dn_in_soot
3865
3866 enddo
3867
3868
3869 n_in_bio_ultra = 0.
3870 dn_in_bio = 0.0
3871
3872
3873
3874
3875
3876
3877
3878
3879
3880
3881 n_in_solo = psi_solo*glass_frac*h_frac_solo*n_in_solo_star
3882 dn_in_solo =max(psi_solo*glass_frac* (h_frac_solo*77211.0*100.0/
3883 &rho_aida + n_in_solo_star*dh_frac_solo), 0.0)
3884
3885
3886 else
3887 n_in_ultra = 0.; n_in_dust_ultra = 0.; n_in_soot_ultra =
3888 & 0.; n_in_solo = 0.; n_in_bio_ultra = 0.;
3889 endif
3890
3891
3892
3893 n_in_dust = n_in_dust + n_in_dust_ultra;
3894 n_in_soot = n_in_soot + n_in_soot_ultra;
3895 n_in_bio = n_in_bio + n_in_bio_ultra;
3896
3897
3898
3899
3900! PROBLEM: how to ensure that the frozen fraction does not exceed 1 ?
3901
3902 if (.false.) then
3903 if(n_in_dust + n_in_bio + n_in_soot + n_in_solo > 0.) then
3904
3905 cihenc_dust = n_in_dust - nin_a_nuc_dust
3906 if(cihenc_dust < 0.) cihenc_dust = 0.
3907
3908 cihenc_soot = n_in_soot - nin_a_nuc_soot
3909 if(cihenc_soot < 0.) cihenc_soot = 0.
3910
3911 cihenc_bio = n_in_bio - nin_a_nuc_bio
3912 if(cihenc_bio < 0.) cihenc_bio = 0.
3913
3914 cihenc_solo = n_in_solo - nin_a_nuc_solo
3915 if(cihenc_solo < 0.) cihenc_solo = 0.
3916
3917
3918 n_iw = n_iw + cihenc_dust
3919 nin_a_nuc_dust = nin_a_nuc_dust + cihenc_dust
3920 num_ic_dust_imm = num_ic_dust_imm + cihenc_dust
3921
3922 n_iw = n_iw + cihenc_soot
3923 nin_a_nuc_soot = nin_a_nuc_soot + cihenc_soot
3924 num_ic_soot_imm = num_ic_soot_imm + cihenc_soot
3925
3926 n_iw = n_iw + cihenc_bio
3927 nin_a_nuc_bio = nin_a_nuc_bio + cihenc_bio
3928 num_ic_bio_imm = num_ic_bio_imm + cihenc_bio
3929
3930 n_iw = n_iw + cihenc_solo
3931 nin_a_nuc_solo = nin_a_nuc_solo + cihenc_solo
3932 num_ic_solo_imm = num_ic_solo_imm + cihenc_solo
3933
3934
3935 endif
3936 end if
3937 endif
3938 endif
3939
3940 n_iw = n_in_dust + n_in_bio + n_in_soot + n_in_solo
3941 dnall = dn_in_dust + dn_in_bio + dn_in_soot + dn_in_solo
3942
3943 if (n_iw .gt. zero_par) then
3944 fdust_dep = dble(n_in_dust/n_iw)
3945 end if
3946 nhet_dep = dble(n_in_dust)
3947 nhet_dhf = dble(n_in_solo)
3948
3949
3950
3951
3952 if (( dnall > 0.) .and. (n_iw .gt. 0.0)) then
3953 dsh=max(min(n_iw/dnall, ss_i), 0.005)
3954 else
3955 dsh=0.005
3956 end if
3957
3958
3959 END SUBROUTINE empirical_param_phillips
3960
3961 real function h_1(x, x_1, x_2, hlo)
3962 real, intent(in) :: hlo, x, x_1, x_2
3963
3964 if(x >= x_2) h_1 = 1
3965 if(x <= x_1) h_1 = hlo
3966 if(x > x_1 .and. x < x_2) h_1 = (x - x_1)/(x_2 - x_1)
3967
3968 if( x_2 <= x_1) stop 91919
3969
3970 return
3971 end function
3972
3973
3974 real function h_1_smooth(x, x_1, x_2, hlo, hhi,dh1smooth)
3975 real, intent(in) :: hlo, hhi, x, x_1, x_2
3976 real*4,intent(out) :: dh1smooth
3977 real :: a_0, a_1, a_2, a_3, a, b
3978
3979 if(x >= x_2) h_1_smooth = hhi
3980 if(x <= x_1) h_1_smooth = hlo
3981
3982 if(x >= x_2) dh1smooth = 0.0
3983 if(x <= x_1) dh1smooth = 0.0
3984
3985 if(x > x_1 .and. x < x_2) then
3986 a = 6.*(hlo - hhi)/(x_2**3. - x_1**3. + 3.*(x_2*x_1*x_1 - x_1*
3987 &x_2*x_2) )
3988 a_3 = (a/3.)
3989 a_2 = -(a/2.)*(x_1 + x_2)
3990 a_1 = a*x_2*x_1
3991 b = hlo + a*(x_1**3.)/6. - a*x_1*x_1*x_2/2.
3992 a_0 = b
3993 h_1_smooth = a_0 + a_1*x + a_2*x*x + a_3*x*x*x
3994 dh1smooth = a_1 + 2.0*a_2*x + 3.0*a_3*x*x
3995 endif
3996
3997!H1smooth =min(dH1smooth , 1000000.0)
3998!H1smooth =max(dH1smooth , 0.0)
3999
4000 if( x_2 <= x_1) stop 91919
4001
4002 return
4003 end function
4004
4005
4006
4007
4008
4009! END ICE PARAMETERIZATION DONIF
4010!
4011!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4012
4013
4014
4015 END MODULE aer_cloud
subroutine sstep(nsoil, sh2oin, rhsct, dt, smcmax, cmcmax, zsoil, sice, cmc, rhstt, ai, bi, ci, sh2oout, runoff3, smc)
This subroutine calculates/updates soil moisture content values and canopy moisture content values.
Definition sflx.f:5296
subroutine copy_aer(a, b)
This subroutine handles aer structure.
Definition aer_cloud.F:1145
subroutine ccn_at_super(super, ccn_at_s, nmodes, sig_par, sg_par, tp_par)
Definition aer_cloud.F:1289
subroutine, public aer_cloud_init()
This subroutine calculates.
Definition aer_cloud.F:116
subroutine, public vertical_vel_variance(omeg, lc_turb, tm_gw, pm_gw,
This subroutine calculates subgrid scale distribution of vertical velocity.
Definition aer_cloud.F:1027
subroutine, public init_aer(aerout)
This subroutine initialize aerosol properties in MG sheme.
Definition aer_cloud.F:1182
subroutine, public aerconversion1(aer_mass, aerpr)
This subroutine sets the properties of the aerosol distributions.
Definition aer_cloud.F:930
subroutine, public aerconversion(aer_mass, aerpr, kappa, sulfate, org, bcarbon, dust, seasalt)
This subroutine sets the properties of the aerosol distributions Mass-number conversion based on Bara...
Definition aer_cloud.F:823
subroutine, public getinsubset(typ, aerin, aerout)
This subroutine extracts aerosol props with INactive = typ.
Definition aer_cloud.F:1102
subroutine, public aerosol_activate(tparc_in, pparc_in, sigwparc_in, wparc_ls, aer_props, npre_in, dpre_in, ccn_diagr8, ndropr8, cdncr8, smaxliqr8, incr8, smaxicer8, nheticer8, inimmr8, dinimmr8, ncdepr8, ncdhfr8, sc_icer8, fdust_immr8, fdust_depr8, fdust_dhfr8, nlimr8, use_average_v, ccn_param, in_param, fd_dust, fd_soot, pfrz_inc_r8, sigma_nuc, rhi_cell, nccn)
This subroutine sets the variables needed for the activation subroutines and return the activated dro...
Definition aer_cloud.F:191
subroutine copy_mode(a_out, a_in, mode_in, mode_out)
This subroutine.
Definition aer_cloud.F:1164
subroutine aerconversion_base()
This subrotine sets basic properties of the aerosol size distributions when using GOCART aerosol Mass...
Definition aer_cloud.F:637
subroutine arg_activ(wparc, sigw, nact, smax, nmodes, tp_par, dpg_par, kappa_par, sig_par, temp_par, pres_par)
This subroutine finds the activated droplet number following Abdul-Razzak and Ghan (2000) Abdul_Razza...
Definition aer_cloud.F:1209
subroutine pdfactiv(wparc, sigw, nact, smax, nmodes,
This subroutine calculates the ccn activation fraction according to the nenes and seinfeld (2003) par...
Definition aer_cloud.F:1406
This module contains some of the most frequently used math and physics constants for GCM models.
Definition physcons.F90:36