Line data Source code
1 : ! Include shortname defintions, so that the F77 code does not have to be modified to
2 : ! reference the CARMA structure.
3 : #include "carma_globaer.h"
4 :
5 : !! This routine evaluates particle loss rates due to nucleation <rnuclg>:
6 : !! aerosol freezing only.
7 : !!
8 : !! The parameterization described by Mohler et al., presented at the AMS
9 : !! Cloud physics workshop (2010) is used.
10 : !!
11 : !! The loss rates for all particle elements in a particle group are equal.
12 : !!
13 : !! To avoid nucleation into an evaporating bin, this subroutine must
14 : !! be called after growp, which evaluates evaporation loss rates <evaplg>.
15 : !!
16 : !! @author Chuck Bardeen
17 : !! @version Aug-2010
18 213271518 : subroutine freezaerl_mohler2010(carma, cstate, iz, rc)
19 :
20 : ! types
21 : use carma_precision_mod
22 : use carma_enums_mod
23 : use carma_constants_mod
24 : use carma_types_mod
25 : use carmastate_mod
26 : use carma_mod
27 :
28 : implicit none
29 :
30 : type(carma_type), intent(in) :: carma !! the carma object
31 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
32 : integer, intent(in) :: iz !! z index
33 : integer, intent(inout) :: rc !! return code, negative indicates failure
34 :
35 : ! Local declarations
36 : real(kind=f), parameter :: prenuc = 2.075e33_f * RHO_W / RHO_I
37 : real(kind=f), parameter :: kt0 = 1.6e0_f
38 : real(kind=f), parameter :: dkt0dp = -8.8e0_f
39 : real(kind=f), parameter :: kti = 0.22e0_f
40 : real(kind=f), parameter :: dktidp = -0.17e0_f
41 :
42 : logical :: evapfrom_nucto
43 : integer :: igas ! gas index
44 : integer :: igroup ! group index
45 : integer :: ibin ! bin index
46 : integer :: iepart ! element for condensing group index
47 : integer :: inuc ! nucleating element index
48 : integer :: isol ! solute index of freezing particle
49 : integer :: ienucto ! index of target nucleation element
50 : integer :: ignucto ! index of target nucleation group
51 : integer :: inucto ! index of target nucleation bin
52 : real(kind=f) :: sifreeze
53 : real(kind=f) :: aw
54 : real(kind=f) :: CONTL
55 : real(kind=f) :: CONTH
56 : real(kind=f) :: H2SO4m
57 : real(kind=f) :: WT
58 : real(kind=f) :: volrat
59 : real(kind=f) :: ssi
60 : real(kind=f) :: ssl
61 : real(kind=f) :: rjj
62 : real(kind=f) :: rlogj
63 : real(kind=f) :: daw
64 : real(kind=f) :: riv
65 : real(kind=f) :: vw0
66 : real(kind=f) :: awi
67 : real(kind=f) :: rsi
68 : real(kind=f) :: dmy
69 : real(kind=f) :: rlnt
70 : real(kind=f) :: td
71 : real(kind=f) :: pp
72 : real(kind=f) :: pp2
73 : real(kind=f) :: pp3
74 : real(kind=f) :: vi
75 : real(kind=f) :: fkelv
76 : real(kind=f) :: fkelvi
77 :
78 :
79 213271518 : rc = RC_OK
80 :
81 : ! Aerosol freezing limited to T < 240K
82 213271518 : if (t(iz) <= 240._f) then
83 :
84 : ! Loop over particle groups.
85 371230947 : do igroup = 1,NGROUP
86 :
87 247487298 : igas = inucgas(igroup)
88 247487298 : iepart = ienconc(igroup)
89 247487298 : isol = isolelem(iepart)
90 :
91 371230947 : if (igas .ne. 0) then
92 :
93 : ! Bypass calculation if few particles are present
94 123743649 : if (pconmax(iz,igroup) .gt. FEW_PC) then
95 :
96 : ! Calculate nucleation loss rates. Do not allow nucleation into
97 : ! an evaporating bin.
98 241022720 : do inuc = 1, nnuc2elem(iepart)
99 :
100 120511360 : ienucto = inuc2elem(inuc,iepart)
101 241022720 : if (ienucto /= 0) then
102 120511360 : ignucto = igelem( ienucto )
103 :
104 : ! Only compute nucleation rate for aerosol freezing
105 : !
106 : ! NOTE: If heterogeneous nucleation of glassy aerosols is being used
107 : ! as a nucleation mechanism, then both the heterogeneous nucleation and
108 : ! the homogeneous freezing need to be considered.
109 120511360 : if (iand(inucproc(iepart,ienucto), I_AF_MOHLER_2010) /= 0) then
110 :
111 : ! Loop over particle bins.
112 0 : do ibin = 1, NBIN
113 :
114 0 : ssi = supsati(iz,igas)
115 0 : ssl = supsatl(iz,igas)
116 :
117 : ! Adjust ssi for the Kelvin effect.
118 0 : fkelvi = exp(akelvini(iz,igas) / r(ibin,igroup))
119 0 : ssi = ssi / fkelvi
120 :
121 : ! Calculate approximate critical saturation needed for homogeneous freezing
122 : ! of sulfate aerosols (see Jensen and Toon, GRL, 1994).
123 0 : sifreeze = 0.3_f
124 :
125 : ! Homogeneous freezing of sulfate aerosols should only occur if SL < Scrit
126 : ! and SI > <sifreeze>.
127 0 : if (ssi > sifreeze) then
128 :
129 : ! Mohler et al. 2010? nucleation rate parameterization
130 0 : rlogj = 97.973292_f - 154.67476_f * (ssi + 1._f) - 0.84952712_f * t(iz) + 1.0049467_f * (ssi + 1._f) * t(iz)
131 0 : rjj = 10._f**(rlogj) ! [cm-3 s-1]
132 :
133 : ! NOTE: The weight percent can become negative from this parameterization,
134 : ! which is not physicsal. With small supersaturations, the water activity
135 : ! becomes postive (>1.013) the weight percent becomes negative. Don't allow
136 : ! the the supsatl to be greater than 0.
137 0 : ssl = max(-1.0_f, min(0._f, ssl))
138 :
139 : ! Kelvin effect on water activity
140 0 : aw = 1._f + ssl ! ?
141 0 : fkelv = exp(akelvin(iz,igas) / r(ibin,igroup))
142 0 : aw = aw / fkelv
143 :
144 : ! Calculate volume ratio of wet/dry aerosols
145 0 : if (aw < 0.05_f) then
146 0 : CONTL = 12.37208932_f * (aw**(-0.16125516114_f)) - 30.490657554_f * aw - 2.1133114241_f
147 0 : CONTH = 13.455394705_f * (aw**(-0.1921312255_f)) - 34.285174604_f * aw - 1.7620073078_f
148 0 : elseif (aw <= 0.85_f) then
149 0 : CONTL = 11.820654354_f * (aw**(-0.20786404244_f)) - 4.807306373_f * aw - 5.1727540348_f
150 0 : CONTH = 12.891938068_f * (aw**(-0.23233847708_f)) - 6.4261237757_f * aw - 4.9005471319_f
151 : else
152 0 : CONTL = -180.06541028_f * (aw**(-0.38601102592_f)) - 93.317846778_f * aw + 273.88132245_f
153 0 : CONTH = -176.95814097_f * (aw**(-0.36257048154_f)) - 90.469744201_f * aw + 267.45509988_f
154 : endif
155 :
156 0 : H2SO4m = CONTL + ((CONTH - CONTL) * (t(iz) - 190._f) / 70._f)
157 0 : WT = (98.0_f * H2SO4m) / (1000._f + 98._f * H2SO4m)
158 0 : WT = max(0._f, min(1._f, WT))
159 0 : WT = 100._f * WT
160 :
161 : ! Volume ratio of wet/dry aerosols.
162 0 : if (WT <= 0._f) then
163 : volrat = 1.e10_f
164 : else
165 0 : volrat = rhosol(isol) / RHO_W * ((100._f - WT) / WT) + 1._f
166 : endif
167 :
168 : ! NOTE: Limit the rate for stability.
169 : ! [s-1]
170 0 : rnuclg(ibin,igroup,ignucto) = rnuclg(ibin,igroup,ignucto) + min(1e20_f, rjj * volrat * vol(ibin,igroup))
171 : endif ! ssi > sifreeze .and. target droplets not evaporating
172 : enddo ! ibin = 1,NBIN
173 : endif ! inucproc(iepart,ienucto) .eq. I_DROPACT
174 : endif
175 : enddo ! inuc = 1,nnuc2elem(iepart)
176 : endif ! pconmax .gt. FEW_PC
177 : endif ! (igas = inucgas(igroup) .ne. 0)
178 : enddo ! igroup = 1,NGROUP
179 : endif
180 :
181 : ! Return to caller with particle loss rates due to nucleation evaluated.
182 213271518 : return
183 213271518 : end subroutine
|
