forked from ESCOMP/CARMA_base
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwetr.F90
305 lines (242 loc) · 11.1 KB
/
wetr.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
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
! Include shortname defintions, so that the F77 code does not have to be modified to
! reference the CARMA structure.
#include "carma_globaer.h"
module wetr
use carma_precision_mod
contains
!! This routine calculates the wet radius for hydrophilic particles that are
!! assumed to grow in size based upon the realtive humidity.
!!
!! Parameterizations based upon Fitzgerald [1975] and Gerber [1985] are support and the
!! particles are assumed to be spherical.
!!
!! Hygroscopicity routine MUST be called prior to getwetr to update kappa for
!! Yu et al. [JAMES, 2015] parameterization (irhswell = I_PETTERS)
!!
!! @author Chuck Bardeen, Pete Colarco
!! @version May-2009 from Nov-2000
subroutine getwetr(carma, igroup, rh, rdry, rwet, rhopdry, rhopwet, rc, h2o_mass, h2o_vp, temp, kappa, wgtpct)
! types
use carma_precision_mod
use carma_enums_mod
use carma_constants_mod
use carma_types_mod
use carmastate_mod
use carma_mod
use sulfate_utils
implicit none
type(carma_type), intent(in) :: carma !! the carma object
integer, intent(in) :: igroup !! group index
real(kind=f), intent(in) :: rh !! relative humidity
real(kind=f), intent(in) :: rdry !! dry radius [cm]
real(kind=f), intent(out) :: rwet !! wet radius [cm]
real(kind=f), intent(in) :: rhopdry !! dry radius [cm]
real(kind=f), intent(out) :: rhopwet !! wet radius [cm]
integer, intent(inout) :: rc !! return code, negative indicates failure
real(kind=f), intent(in), optional :: h2o_mass!! water vapor mass concentration (g/cm3)
real(kind=f), intent(in), optional :: h2o_vp !! water eq. vaper pressure (dynes/cm2)
real(kind=f), intent(in), optional :: temp !! temperature [K]
real(kind=f), intent(in), optional :: kappa !! hygroscopicity parameter (Petters & Kreidenweis, ACP, 2007)
real(kind=f), intent(in), optional :: wgtpct !! weight percent h2so4 (overrides rh)
! Local declarations
real(kind=f) :: humidity
real(kind=f) :: r_ratio
real(kind=f) :: wtpkelv, den1, den2, drho_dwt
real(kind=f) :: sigkelv, sig1, sig2, dsigma_dwt
real(kind=f) :: rkelvinH2O_a, rkelvinH2O_b, rkelvinH2O, h2o_kelv
real(kind=f) :: rh190
! The following parameters relate to the swelling of seasalt like particles
! following Fitzgerald, Journal of Applied Meteorology, [1975].
!
! Question - Should epsilon be 1._f? It means alpharat is 1 by definition.
real(kind=f), parameter :: epsilon_ = 1._f ! soluble fraction of deliquescing particle
real(kind=f) :: alphaComp
real(kind=f) :: alpha
real(kind=f) :: alpha1
real(kind=f) :: alpharat
real(kind=f) :: beta
real(kind=f) :: theta
real(kind=f) :: f1
real(kind=f) :: f2
! Parameters from Gerber [1985]
real(kind=f) :: c1
real(kind=f) :: c2
real(kind=f) :: c3
real(kind=f) :: c4
! Define formats
1 format(/,'Non-spherical particles specified for group ',i3, ' (ishape=',i3,') but spheres assumed in wetr.f.'/)
! If humidty affects the particle, then determine the equilbirium
! radius and density based upon the relative humidity.
if (irhswell(igroup) == I_NO_SWELLING) then
! No swelling, just use the dry values.
rwet = rdry
rhopwet = rhopdry
else ! irhswell(igroup) /= I_NO_SWELLING
! Warning message for non-spherical particles!
if( ishape(igroup) .ne. I_SPHERE )then
if (do_print) write(LUNOPRT,1) igroup, ishape(igroup)
rc = RC_ERROR
return
endif
! The Parameterizations don't handle relative humidities of 0, and
! behave poorly when RH > 0.995, so cap the relative humidity
! used to these values.
humidity = min(max(rh,tiny(1.0_f)), 0.995_f)
! Fitzgerald Parameterization
if (irhswell(igroup) == I_FITZGERALD) then
! Calculate the alpha and beta parameters for the wet particle
! relative to amonium sulfate
beta = exp((0.00077_f * humidity) / (1.009_f - humidity))
if (humidity .le. 0.97_f) then
theta = 1.058_f
else
theta = 1.058_f - (0.0155_f * (humidity - 0.97_f)) / (1.02_f - humidity**1.4_f)
endif
alpha1 = 1.2_f * exp((0.066_f * humidity) / (theta - humidity))
f1 = 10.2_f - 23.7_f * humidity + 14.5_f * humidity**2
f2 = -6.7_f + 15.5_f * humidity - 9.2_f * humidity**2
alpharat = 1._f - f1 * (1._f - epsilon_) - f2 * (1._f - epsilon_**2)
! Scale the size based on the composition of the particle.
select case(irhswcomp(igroup))
case (I_SWF_NH42SO4)
alphaComp = 1.00_f
case(I_SWF_NH4NO3)
alphaComp = 1.06_f
case(I_SWF_NANO3)
alphaComp = 1.17_f
case(I_SWF_NH4CL)
alphaComp = 1.23_f
case(I_SWF_CACL2)
alphaComp = 1.29_f
case(I_SWF_NABR)
alphaComp = 1.32_f
case(I_SWF_NACL)
alphaComp = 1.35_f
case(I_SWF_MGCL2)
alphaComp = 1.41_f
case(I_SWF_LICL)
alphaComp = 1.54_f
case default
if (do_print) write(LUNOPRT,*) "wetr:: ERROR - Unknown composition type (", irhswcomp(igroup), &
") for Fitzgerald."
rc = RC_ERROR
return
end select
alpha = alphaComp * (alpha1 * alpharat)
! Determine the wet radius.
!
! NOTE: Fitgerald's equations assume r in [um], so scale the cgs units
! appropriately.
rwet = (alpha * (rdry * 1e4_f)**beta) * (1e-4_f)
! Determine the wet density from the wet radius.
r_ratio = (rdry / rwet)**3._f
rhopwet = r_ratio * rhopdry + (1._f - r_ratio) * RHO_W
end if
! Gerber Paremeterization
if (irhswell(igroup) == I_GERBER) then
! Scale the size based on the composition of the particle.
select case(irhswcomp(igroup))
case (I_SWG_NH42SO4)
c1 = 0.4809_f
c2 = 3.082_f
c3 = 3.110e-11_f
c4 = -1.428_f
case(I_SWG_URBAN)
c1 = 0.3926_f
c2 = 3.101_f
c3 = 4.190e-11_f
c4 = -1.404_f
case(I_SWG_RURAL)
c1 = 0.2789_f
c2 = 3.115_f
c3 = 5.415e-11_f
c4 = -1.399_f
case(I_SWG_SEA_SALT)
c1 = 0.7674_f
c2 = 3.079_f
c3 = 2.572e-11_f
c4 = -1.424_f
case default
if (do_print) write(LUNOPRT,*) "wetr:: ERROR - Unknown composition type (", irhswcomp(igroup), &
") for Gerber."
rc = RC_ERROR
return
end select
rwet = ((c1 * rdry**c2 / (c3 * rdry**c4 - log10(humidity))) + rdry**3)**(1._f / 3._f)
! Determine the wet density from the wet radius.
r_ratio = (rdry / rwet)**3
rhopwet = r_ratio * rhopdry + (1._f - r_ratio) * RHO_W
end if ! irhswell(igroup) == I_GERBER
! Mixed aerosol paremeterization (Pengfei Yu et al., JAMES, 2015) based on
! Petters and Kreidenweis (ACP, 2007) hygroscopicity parameter kappa
if (irhswell(igroup) == I_PETTERS) then
if (.not.( present(temp) .and. &
present(kappa) )) then
if (do_print) write(LUNOPRT,*) "wetr:: ERROR - h2o_mass,temp,kappa for PETTERS"
rc = RC_ERROR
return
endif
if (temp .le. 190._f) then
rh190 = rh * pvap_h2o(temp) / pvap_h2o(190._f)
rh190 = min(max(rh190,tiny(1.0_f)), 0.995_f)
rwet = rdry * (1._f + rh190*kappa/(1._f-rh190))**(1._f/3._f)
else ! temp > 190
rwet = rdry * (1._f + humidity*kappa/(1._f-humidity))**(1._f/3._f)
end if
r_ratio = (rdry / rwet)**3._f
rhopwet = r_ratio * rhopdry + (1._f - r_ratio) * RHO_W
end if ! irhswell(igroup) == I_PETTERS
! Sulfate Aerosol, using weight percent.
if (irhswell(igroup) == I_WTPCT_H2SO4) then
! Has the weight percent already been specified, or do we need to determine
! in based upon the water and temperature.
!
! NOTE: This is used when generating optical properties files.
if (present(wgtpct) .and. present(temp)) then
rhopwet = sulfate_density(carma, wgtpct, temp, rc)
rwet = rdry * (100._f * rhopdry / wgtpct / rhopwet)**(1._f / 3._f)
else if (.not.( present(h2o_mass) .and. &
present(h2o_vp) .and. &
present(temp) )) then
if (do_print) write(LUNOPRT,*) "wetr:: ERROR - h2o_mass,h2o_vp,temp for WTPCT_H2SO4"
rc = RC_ERROR
return
else
! Adjust calculation for the Kelvin effect of H2O:
wtpkelv = 80._f ! start with assumption of 80 wt % H2SO4
den1 = 2.00151_f - 0.000974043_f * temp ! density at 79 wt %
den2 = 2.01703_f - 0.000988264_f * temp ! density at 80 wt %
drho_dwt = den2-den1 ! change in density for change in 1 wt %
sig1 = 79.3556_f - 0.0267212_f * temp ! surface tension at 79.432 wt %
sig2 = 75.608_f - 0.0269204_f * temp ! surface tension at 85.9195 wt %
dsigma_dwt = (sig2-sig1) / (85.9195_f - 79.432_f) ! change in density for change in 1 wt %
sigkelv = sig1 + dsigma_dwt * (80.0_f - 79.432_f)
rwet = rdry * (100._f * rhopdry / wtpkelv / den2)**(1._f / 3._f)
rkelvinH2O_b = 1._f + wtpkelv * drho_dwt / den2 - 3._f * wtpkelv &
* dsigma_dwt / (2._f*sigkelv)
rkelvinH2O_a = 2._f * gwtmol(igash2so4) * sigkelv / (den1 * RGAS * temp * rwet)
rkelvinH2O = exp (rkelvinH2O_a*rkelvinH2O_b)
h2o_kelv = h2o_mass / rkelvinH2O
wtpkelv = wtpct_tabaz(carma, temp, h2o_kelv, h2o_vp, rc)
rhopwet = sulfate_density(carma, wtpkelv, temp, rc)
rwet = rdry * (100._f * rhopdry / wtpkelv / rhopwet)**(1._f / 3._f)
endif
end if ! irhswell(igroup) == I_WTPCT_H2SO4
end if ! irhswell(igroup) /= I_NO_SWELLING
! Return to caller with wet radius evaluated.
return
contains
!! This function is needed for generating wet radii for optics when using the
!! PETTERS scheme, and should not be used generally within CARMA.
!!
!! The vaporp_h2o_murphy2005 equation to calculate the vapor pressure at 190 K
!! for liquid water
pure function pvap_h2o(temp) result(pvap)
real(kind=f), intent(in) :: temp
real(kind=f) :: pvap
pvap = temp / ( 10.0_f * exp(54.842763_f - (6763.22_f / temp) &
- (4.210_f * log(temp)) + (0.000367_f * temp) + (tanh(0.0415_f * (temp - 218.8_f)) &
* (53.878_f - (1331.22_f / temp) - (9.44523_f * log(temp)) + 0.014025_f * temp))) )
end function pvap_h2o
end subroutine
end module