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 Tabazadeh et al. [GRL, 27, 1111, 2000.] is
9 : !! used.
10 : !!
11 : !! The loss rates for all particle elements in a particle group are equal.
12 : !!
13 : !! @author Eric Jensen, Chuck Bardeen
14 : !! @version Mar-1995, Nov-2009
15 1418703438 : subroutine freezaerl_tabazadeh2000(carma, cstate, iz, rc)
16 :
17 : ! types
18 : use carma_precision_mod
19 : use carma_enums_mod
20 : use carma_constants_mod
21 : use carma_types_mod
22 : use carmastate_mod
23 : use carma_mod
24 :
25 : implicit none
26 :
27 : type(carma_type), intent(in) :: carma !! the carma object
28 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
29 : integer, intent(in) :: iz !! z index
30 : integer, intent(inout) :: rc !! return code, negative indicates failure
31 :
32 : ! Local declarations
33 : ! Define parameters needed for freezing nucleation calculations.
34 : ! real(kind=f), parameter :: adelf = 1.29e-12_f
35 : ! real(kind=f), parameter :: bdelf = 0.05_f
36 : real(kind=f), parameter :: prenuc = 2.075e33_f * RHO_W / RHO_I
37 : ! real(kind=f), parameter :: rmiv = 0.6_f
38 :
39 : integer :: igas !! gas index
40 : integer :: igroup !! group index
41 : integer :: ibin !! bin index
42 : integer :: iepart !! element for condensing group index
43 : integer :: inuc !! nucleating element index
44 : integer :: ienucto !! index of target nucleation element
45 : integer :: ignucto !! index of target nucleation group
46 : integer :: isol
47 : real(kind=f) :: A0, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10
48 : real(kind=f) :: c0, C1, C2, C3, C4, c5
49 : real(kind=f) :: d0, d1, d2, d3, d4, d5
50 : real(kind=f) :: e0, e1, e2, e3, e4, e5
51 : real(kind=f) :: sifreeze
52 : real(kind=f) :: rhoibar
53 : real(kind=f) :: rlhbar
54 : real(kind=f) :: act
55 : real(kind=f) :: CONTL
56 : real(kind=f) :: CONTH
57 : real(kind=f) :: H2SO4m
58 : real(kind=f) :: WT
59 : real(kind=f) :: vrat
60 : real(kind=f) :: wtfrac
61 : real(kind=f) :: den
62 : real(kind=f) :: diffact
63 : real(kind=f) :: S260, S220, S180
64 : real(kind=f) :: sigma
65 : real(kind=f) :: sigsula
66 : real(kind=f) :: sigicea
67 : real(kind=f) :: sigsulice
68 : real(kind=f) :: ag
69 : real(kind=f) :: delfg
70 : real(kind=f) :: expon
71 : real(kind=f) :: ssl
72 : real(kind=f) :: fkelv
73 :
74 :
75 : ! Loop over particle groups.
76 4256110314 : do igroup = 1,NGROUP
77 :
78 2837406876 : igas = inucgas(igroup)
79 2837406876 : iepart = ienconc(igroup)
80 2837406876 : isol = isolelem(iepart)
81 :
82 4256110314 : if( igas .ne. 0 )then
83 :
84 : ! Calculate nucleation loss rates. Do not allow nucleation into
85 : ! an evaporating bin.
86 : ! if( nnuc2elem(iepart) .gt. 1 )then
87 2837406876 : do inuc = 1,nnuc2elem(iepart)
88 :
89 1418703438 : ienucto = inuc2elem(inuc,iepart)
90 2837406876 : if( ienucto .ne. 0 )then
91 1418703438 : ignucto = igelem( ienucto )
92 :
93 : ! Only compute nucleation rate for aerosol freezing.
94 : !
95 : ! NOTE: If heterogeneous nucleation of glassy aerosols is being used
96 : ! as a nucleation mechanism, then both the heterogeneous nucleation and
97 : ! the homogeneous freezing need to be considered.
98 1418703438 : if ((iand(inucproc(iepart,ienucto), I_AF_TABAZADEH_2000) /= 0)) then
99 :
100 : ! Loop over particle bins. Loop from largest to smallest for
101 : ! evaluation of index of smallest bin nucleated during time step <inucstep>.
102 0 : do ibin =NBIN,1,-1
103 :
104 : ! Bypass calculation if few particles are present
105 0 : if( pconmax(iz,igroup) .gt. FEW_PC )then
106 :
107 : ! Calculate approximate critical saturation needed for homogeneous freezing
108 : ! of sulfate aerosols (see Jensen and Toon, GRL, 1994).
109 0 : sifreeze = 0.3_f
110 :
111 : ! NOTE: The wieght percent can become negative from this parameterization,
112 : ! which is not physicsal. With small supersaturations, the water activity
113 : ! becomes postive (>1.013) the weight percent becomes negative. Don't allow
114 : ! the the supsatl to be greater than 0.
115 0 : ssl = max(-1.0_f, min(0._f, supsatl(iz,igas)))
116 :
117 :
118 : ! Homogeneous freezing of sulfate aerosols should only occur of SL < Scrit
119 : ! and SI > <sifreeze>.
120 0 : if( supsati(iz,igas) .gt. sifreeze)then
121 :
122 : ! Calculate mean ice density and latent heat of freezing over temperature
123 : ! interval [T0,T]
124 :
125 0 : rhoibar = ( 0.916_f * (t(iz)-T0) - &
126 : 1.75e-4_f/2._f * ((t(iz)-T0)**2) - &
127 0 : 5.e-7_f * ((t(iz)-T0)**3)/3._f ) / (t(iz)-T0)
128 :
129 : rlhbar = ( 79.7_f * (t(iz)-T0) + &
130 : 0.485_f/2._f * (t(iz)-T0)**2 - &
131 : 2.5e-3_f/3._f * (t(iz)-T0)**3 ) &
132 0 : / (t(iz)-T0) * 4.186e7_f*18._f
133 :
134 : ! Equilibrium H2SO4 weight percent for fixed water activity
135 0 : act = min(1.0_f, ssl + 1._f)
136 :
137 : ! Kelvin effect on water activity
138 0 : fkelv = exp(akelvin(iz,igas) / r(ibin,igroup)) ! ?
139 0 : act = act / fkelv
140 :
141 0 : IF(act .LT. 0.05_f) THEN
142 : CONTL = 12.37208932_f * (act**(-0.16125516114_f)) - &
143 0 : 30.490657554_f * act - 2.1133114241_f
144 : CONTH = 13.455394705_f * (act**(-0.1921312255_f)) - &
145 0 : 34.285174604_f * act - 1.7620073078_f
146 : END IF
147 0 : IF(act .GE. 0.05_f .and. act .LE. 0.85_f) THEN
148 : CONTL = 11.820654354_f * (act**(-0.20786404244_f)) - &
149 0 : 4.807306373_f * act - 5.1727540348_f
150 : CONTH = 12.891938068_f * (act**(-0.23233847708_f)) - &
151 0 : 6.4261237757_f * act - 4.9005471319_f
152 : END IF
153 0 : IF(act .GT. 0.85_f) THEN
154 : CONTL = -180.06541028_f * (act**(-0.38601102592_f)) - &
155 0 : 93.317846778_f * act + 273.88132245_f
156 : CONTH = -176.95814097_f * (act**(-0.36257048154_f)) - &
157 0 : 90.469744201_f * act + 267.45509988_f
158 : END IF
159 0 : H2SO4m = CONTL + ((CONTH - CONTL) * (t(iz) -190._f)/70._f)
160 0 : WT = (98.0_f * H2SO4m)/(1000._f + 98._f * H2SO4m)
161 0 : WT = 100._f * WT
162 :
163 : ! Volume ratio of wet/dry aerosols.
164 0 : vrat = rhosol(isol)/RHO_W * ((100._f-wt)/wt) + 1._f
165 :
166 : ! Calculation sulfate solution density from Myhre et al. (1998).
167 0 : wtfrac = WT/100._f
168 0 : C1 = t(iz) - 273.15_f
169 0 : C2 = C1**2
170 0 : C3 = C1**3
171 0 : C4 = C1**4
172 0 : A0 = 999.8426_f + 334.5402e-4_f*C1 - 569.1304e-5_f*C2
173 : A1 = 547.2659_f - 530.0445e-2_f*C1 + 118.7671e-4_f*C2 &
174 0 : + 599.0008e-6_f*C3
175 : A2 = 526.295e+1_f + 372.0445e-1_f*C1 + 120.1909e-3_f*C2 &
176 0 : - 414.8594e-5_f*C3 + 119.7973e-7_f*C4
177 : A3 = -621.3958e+2_f - 287.7670_f*C1 - 406.4638e-3_f*C2 &
178 0 : + 111.9488e-4_f*C3 + 360.7768e-7_f*C4
179 : A4 = 409.0293e+3_f + 127.0854e+1_f*C1 + 326.9710e-3_f*C2 &
180 0 : - 137.7435e-4_f*C3 - 263.3585e-7_f*C4
181 : A5 = -159.6989e+4_f - 306.2836e+1_f*C1 + 136.6499e-3_f*C2 &
182 0 : + 637.3031e-5_f*C3
183 0 : A6 = 385.7411e+4_f + 408.3717e+1_f*C1 - 192.7785e-3_f*C2
184 0 : A7 = -580.8064e+4_f - 284.4401e+1_f*C1
185 0 : A8 = 530.1976e+4_f + 809.1053_f*C1
186 0 : A9 = -268.2616e+4_f
187 0 : A10 = 576.4288e+3_f
188 : den = A0 + wtfrac*A1 + wtfrac**2 * A2 + &
189 0 : wtfrac**3 * A3 + wtfrac**4 * A4
190 : den = den + wtfrac**5 * A5 + wtfrac**6 * A6 + &
191 0 : wtfrac**7 * A7
192 : den = den + wtfrac**8 * A8 + wtfrac**9 * A9 + &
193 0 : wtfrac**10 * A10
194 :
195 : ! Activation energy is based on Koop's lab data.
196 0 : IF(t(iz) .GT. 220._f) then
197 : A0 = 104525.93058_f
198 : A1 = -1103.7644651_f
199 : A2 = 1.070332702_f
200 : A3 = 0.017386254322_f
201 : A4 = -1.5506854268e-06_f
202 : A5 = -3.2661912497e-07_f
203 : A6 = 6.467954459e-10_f
204 : ELSE
205 0 : A0 = -17459.516183_f
206 0 : A1 = 458.45827551_f
207 0 : A2 = -4.8492831317_f
208 0 : A3 = 0.026003658878_f
209 0 : A4 = -7.1991577798e-05_f
210 0 : A5 = 8.9049094618e-08_f
211 0 : A6 = -2.4932257419e-11_f
212 : END IF
213 :
214 : diffact = ( A0 + A1*t(iz) + A2*t(iz)**2 + &
215 : A3*t(iz)**3 + A4*t(iz)**4 + &
216 0 : A5*t(iz)**5 + A6*t(iz)**6 ) * 1.0e-13_f
217 :
218 : ! Surface energy
219 :
220 : ! Weight percent function for T = 260 K
221 0 : c0 = 77.40682664_f
222 0 : c1 = -0.006963123274_f
223 0 : c2 = -0.009682499074_f
224 0 : c3 = 0.00088797988_f
225 0 : c4 = -2.384669516e-05_f
226 0 : c5 = 2.095358048e-07_f
227 : S260 = c0 + c1*wt + c2*wt**2 + c3*wt**3 + &
228 0 : c4*wt**4 + c5*wt**5
229 :
230 : ! Weight percent function for T = 220 K
231 0 : d0 = 82.01197792_f
232 0 : d1 = 0.5312072092_f
233 0 : d2 = -0.1050692123_f
234 0 : d3 = 0.005415260617_f
235 0 : d4 = -0.0001145573827_f
236 0 : d5 = 8.969257061e-07_f
237 : S220 = d0 + d1*wt + d2*wt**2 + d3*wt**3 + &
238 0 : d4*wt**4 + d5*wt**5
239 :
240 : ! Weight percent function for T = 180K
241 0 : e0 = 85.75507114_f
242 0 : e1 = 0.09541966318_f
243 0 : e2 = -0.1103647657_f
244 0 : e3 = 0.007485866933_f
245 0 : e4 = -0.0001912224154_f
246 0 : e5 = 1.736789787e-06_f
247 : S180 = e0 + e1*wt + e2*wt**2 + e3*wt**3 + &
248 0 : e4*wt**4 + e5*wt**5
249 :
250 0 : if( t(iz) .GE. 220._f ) then
251 0 : sigma = S260 + ((260._f-t(iz))*(S220-S260))/40._f
252 : else
253 0 : sigma = S220 + ((220._f-t(iz))*(S180-S220))/40._f
254 : endif
255 :
256 0 : sigsula = sigma
257 0 : sigicea = 105._f
258 0 : sigsulice = abs( sigsula - sigicea )
259 :
260 : ! Critical ice germ radius formed in the sulfate solution
261 0 : ag = 2._f*gwtmol(igas)*sigsulice / &
262 : ( rlhbar * rhoibar * log(T0/t(iz)) + &
263 : rhoibar * rgas * 0.5_f * (T0+t(iz)) * &
264 0 : log(ssl+1._f) )
265 :
266 0 : if( ag .lt. 0._f ) ag = 1.e10_f
267 :
268 : ! Gibbs free energy of ice germ formation in the ice/sulfate solution
269 0 : delfg = 4._f/3._f*PI * sigsulice * (ag**2)
270 :
271 : ! Ice nucleation rate in a 0.2 micron aerosol (/sec)
272 0 : expon = ( -diffact - delfg ) / BK / t(iz)
273 0 : expon = max( -100._f*ONE, expon )
274 0 : rnuclg(ibin,igroup,ignucto) = prenuc * &
275 : sqrt(sigsulice*t(iz)) * &
276 0 : vrat*vol(ibin,igroup) * exp( expon )
277 :
278 : ! This parameterizations has problems that sometimes yield negative nucleation
279 : ! rates. It would be best to fix the parameterization, but at least keep negative
280 : ! values from being return.
281 0 : if (rnuclg(ibin,igroup,ignucto) < 0._f) then
282 0 : rnuclg(ibin,igroup,ignucto) = 0._f
283 : end if
284 :
285 :
286 : ! xh = 0.1 * r(ibin,igroup) / ag
287 : ! phih = sqrt( 1. - 2.*rmiv*xh + xh**2 )
288 : ! rath = (xh-rmiv) / phih
289 : ! fv3h = xh**3 * ( 2.*ONE - 3.*rath + rath**3 )
290 : ! fv4h = 3. * rmiv * xh**2 * (rath-1.)
291 : ! if( abs(rath) .gt. 1.e0-1.e-8 ) fv3h = 0.
292 : ! if( abs(rath) .gt. 1.e0-1.e-10 ) fv4h = 0.
293 : !
294 : ! fh = 0.5 * ( ONE + ((ONE-rmiv*xh)/phih)**3 +
295 : ! $ fv3h + fv4h )
296 : !
297 : ! expon = ( -delfwat2ice - delfg ) / BK / t3(ixyz)
298 : ! expon = max( -POWMAX, expon )
299 : endif
300 : endif ! pconmax(ixyz,igroup) .gt. FEW_PC
301 : enddo ! ibin = 1,NBIN
302 : endif ! inucproc(iepart,ienucto) .eq. I_DROPACT
303 : endif
304 : enddo ! inuc = 1,nnuc2elem(iepart)
305 : ! endif ! (nnuc2elem(iepart) .gt. 1)
306 : endif ! (igas = inucgas(igroup) .ne. 0)
307 : enddo ! igroup = 1,NGROUP
308 :
309 : ! Return to caller with particle loss rates due to nucleation evaluated.
310 1418703438 : return
311 1418703438 : end
|
