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 Koop et al., Nature 406, 611-614, 2000
9 : !! 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 Eric Jensen, Chuck Bardeen
17 : !! @version Dec-2003, Apr-2010
18 213271518 : subroutine freezaerl_koop2000(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_KOOP_2000) /= 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 : ! Calculate approximate critical saturation needed for homogeneous freezing
118 : ! of sulfate aerosols (see Jensen and Toon, GRL, 1994).
119 0 : sifreeze = 0.3_f
120 :
121 : ! Homogeneous freezing of sulfate aerosols should only occur if SL < Scrit
122 : ! and SI > <sifreeze>.
123 0 : if (ssi > sifreeze) then
124 :
125 : ! Koop et al. nucleation rate parameterization
126 0 : td = t(iz)
127 0 : rlnt = log(td)
128 : ! eqn 2, potential difference [J mol-1]
129 0 : dmy = 210368._f + 131.438_f * td - (3.32373e6_f / td) - 41729.1_f * rlnt
130 0 : rsi = RGAS / 1.e7_f ! gas constant [J mol-1 K-1]
131 : ! Notes (p: ambient vs. at pressure) ?
132 0 : awi = exp(dmy / (rsi * td))
133 :
134 : ! eqn 4
135 0 : vw0 = -230.76_f - 0.1478_f * td + (4099.2_f / td) + 48.8341_f * rlnt
136 : ! eqn 5
137 0 : vi = 19.43_f - 2.2e-3_f * td + 1.08e-5_f * td * td
138 :
139 0 : pp = 1.e-10_f * p(iz) ! pressure [GPa]
140 0 : pp2 = pp * pp * 0.5_f
141 0 : pp3 = pp2 * pp / 3._f
142 0 : riv = vw0 * (pp - kt0 * pp2 - dkt0dp * pp3) - vi * (pp - kti * pp2 - dktidp * pp3) ! eqn 3
143 :
144 0 : riv = riv * 1.e3_f ! [GPa cm3 mol-1] to [Pa m3 mol-1]
145 :
146 : ! NOTE: The wieght percent can become negative from this parameterization,
147 : ! which is not physicsal. With small supersaturations, the water activity
148 : ! becomes postive (>1.013) the weight percent becomes negative. Don't allow
149 : ! the the supsatl to be greater than 0.
150 0 : ssl = max(-1.0_f, min(0._f, ssl))
151 :
152 : ! Water activity
153 0 : aw = 1._f + ssl ! ?
154 :
155 : ! Kelvin effect on water activity
156 0 : fkelv = exp(akelvin(iz,igas) / r(ibin,igroup)) ! ?
157 0 : aw = aw / fkelv
158 :
159 : ! Nucleation rate
160 : !
161 : ! NOTE: This formulation is only valid for daw in the range of
162 : ! 0.26 < daw < 0.34, so limit daw to that range.
163 0 : daw = aw * exp(riv / (rsi*td)) - awi ! eqn 6
164 0 : daw = min(0.34_f, max(daw, 0.26_f)) ! eqn 7
165 :
166 0 : rlogj = ((29180._f * daw - 26924._f) * daw + 8502._f) * daw - 906.7_f ! eqn 7
167 0 : rlogj = min(rlogj, POWMAX*0.3_f)
168 0 : rjj = 10._f**(rlogj) ! [cm-3 s-1]
169 :
170 :
171 : ! Calculate volume ratio of wet/dry aerosols
172 0 : if (aw < 0.05_f) then
173 0 : CONTL = 12.37208932_f * (aw**(-0.16125516114_f)) - 30.490657554_f * aw - 2.1133114241_f
174 0 : CONTH = 13.455394705_f * (aw**(-0.1921312255_f)) - 34.285174604_f * aw - 1.7620073078_f
175 0 : elseif (aw <= 0.85_f) then
176 0 : CONTL = 11.820654354_f * (aw**(-0.20786404244_f)) - 4.807306373_f * aw - 5.1727540348_f
177 0 : CONTH = 12.891938068_f * (aw**(-0.23233847708_f)) - 6.4261237757_f * aw - 4.9005471319_f
178 : else
179 0 : CONTL = -180.06541028_f * (aw**(-0.38601102592_f)) - 93.317846778_f * aw + 273.88132245_f
180 0 : CONTH = -176.95814097_f * (aw**(-0.36257048154_f)) - 90.469744201_f * aw + 267.45509988_f
181 : endif
182 :
183 0 : H2SO4m = CONTL + ((CONTH - CONTL) * (t(iz) - 190._f) / 70._f)
184 0 : WT = (98.0_f * H2SO4m) / (1000._f + 98._f * H2SO4m)
185 0 : WT = max(0._f, min(1._f, WT))
186 0 : WT = 100._f * WT
187 :
188 : ! Volume ratio of wet/dry aerosols.
189 0 : if (WT <= 0._f) then
190 : volrat = 1.e10_f
191 : else
192 0 : volrat = rhosol(isol) / RHO_W * ((100._f - WT) / WT) + 1._f
193 : endif
194 :
195 : ! [s-1]
196 0 : rnuclg(ibin,igroup,ignucto) = rnuclg(ibin,igroup,ignucto) + rjj * volrat * vol(ibin,igroup)
197 : endif ! ssi > sifreeze .and. target droplets not evaporating
198 : enddo ! ibin = 1,NBIN
199 : endif ! inucproc(iepart,ienucto) .eq. I_DROPACT
200 : endif
201 : enddo ! inuc = 1,nnuc2elem(iepart)
202 : endif ! pconmax .gt. FEW_PC
203 : endif ! (igas = inucgas(igroup) .ne. 0)
204 : enddo ! igroup = 1,NGROUP
205 : endif
206 :
207 : ! Return to caller with particle loss rates due to nucleation evaluated.
208 213271518 : return
209 213271518 : end subroutine
|
