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