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 : module wetr
6 : use carma_precision_mod
7 :
8 : contains
9 :
10 : !! This routine calculates the wet radius for hydrophilic particles that are
11 : !! assumed to grow in size based upon the realtive humidity.
12 : !!
13 : !! Parameterizations based upon Fitzgerald [1975] and Gerber [1985] are support and the
14 : !! particles are assumed to be spherical.
15 : !!
16 : !! Hygroscopicity routine MUST be called prior to getwetr to update kappa for
17 : !! Yu et al. [JAMES, 2015] parameterization (irhswell = I_PETTERS)
18 : !!
19 : !! @author Chuck Bardeen, Pete Colarco
20 : !! @version May-2009 from Nov-2000
21 38135321593 : subroutine getwetr(carma, igroup, rh, rdry, rwet, rhopdry, rhopwet, rc, h2o_mass, h2o_vp, temp, kappa, wgtpct)
22 :
23 : ! types
24 : use carma_precision_mod
25 : use carma_enums_mod
26 : use carma_constants_mod
27 : use carma_types_mod
28 : use carmastate_mod
29 : use carma_mod
30 : use sulfate_utils
31 :
32 : implicit none
33 :
34 : type(carma_type), intent(in) :: carma !! the carma object
35 : integer, intent(in) :: igroup !! group index
36 : real(kind=f), intent(in) :: rh !! relative humidity
37 : real(kind=f), intent(in) :: rdry !! dry radius [cm]
38 : real(kind=f), intent(out) :: rwet !! wet radius [cm]
39 : real(kind=f), intent(in) :: rhopdry !! dry radius [cm]
40 : real(kind=f), intent(out) :: rhopwet !! wet radius [cm]
41 : integer, intent(inout) :: rc !! return code, negative indicates failure
42 : real(kind=f), intent(in), optional :: h2o_mass!! water vapor mass concentration (g/cm3)
43 : real(kind=f), intent(in), optional :: h2o_vp !! water eq. vaper pressure (dynes/cm2)
44 : real(kind=f), intent(in), optional :: temp !! temperature [K]
45 : real(kind=f), intent(in), optional :: kappa !! hygroscopicity parameter (Petters & Kreidenweis, ACP, 2007)
46 : real(kind=f), intent(in), optional :: wgtpct !! weight percent h2so4 (overrides rh)
47 :
48 : ! Local declarations
49 : real(kind=f) :: humidity
50 : real(kind=f) :: r_ratio
51 : real(kind=f) :: wtpkelv, den1, den2, drho_dwt
52 : real(kind=f) :: sigkelv, sig1, sig2, dsigma_dwt
53 : real(kind=f) :: rkelvinH2O_a, rkelvinH2O_b, rkelvinH2O, h2o_kelv
54 : real(kind=f) :: rh190
55 :
56 :
57 : ! The following parameters relate to the swelling of seasalt like particles
58 : ! following Fitzgerald, Journal of Applied Meteorology, [1975].
59 : !
60 : ! Question - Should epsilon be 1._f? It means alpharat is 1 by definition.
61 : real(kind=f), parameter :: epsilon_ = 1._f ! soluble fraction of deliquescing particle
62 : real(kind=f) :: alphaComp
63 : real(kind=f) :: alpha
64 : real(kind=f) :: alpha1
65 : real(kind=f) :: alpharat
66 : real(kind=f) :: beta
67 : real(kind=f) :: theta
68 : real(kind=f) :: f1
69 : real(kind=f) :: f2
70 :
71 : ! Parameters from Gerber [1985]
72 : real(kind=f) :: c1
73 : real(kind=f) :: c2
74 : real(kind=f) :: c3
75 : real(kind=f) :: c4
76 :
77 : ! Define formats
78 : 1 format(/,'Non-spherical particles specified for group ',i3, ' (ishape=',i3,') but spheres assumed in wetr.f.'/)
79 :
80 : ! If humidty affects the particle, then determine the equilbirium
81 : ! radius and density based upon the relative humidity.
82 38135321593 : if (irhswell(igroup) == I_NO_SWELLING) then
83 :
84 : ! No swelling, just use the dry values.
85 0 : rwet = rdry
86 0 : rhopwet = rhopdry
87 :
88 : else ! irhswell(igroup) /= I_NO_SWELLING
89 :
90 : ! Warning message for non-spherical particles!
91 38135321593 : if( ishape(igroup) .ne. I_SPHERE )then
92 0 : if (do_print) write(LUNOPRT,1) igroup, ishape(igroup)
93 0 : rc = RC_ERROR
94 0 : return
95 : endif
96 :
97 : ! The Parameterizations don't handle relative humidities of 0, and
98 : ! behave poorly when RH > 0.995, so cap the relative humidity
99 : ! used to these values.
100 38135321593 : humidity = min(max(rh,tiny(1.0_f)), 0.995_f)
101 :
102 : ! Fitzgerald Parameterization
103 38135321593 : if (irhswell(igroup) == I_FITZGERALD) then
104 :
105 : ! Calculate the alpha and beta parameters for the wet particle
106 : ! relative to amonium sulfate
107 0 : beta = exp((0.00077_f * humidity) / (1.009_f - humidity))
108 0 : if (humidity .le. 0.97_f) then
109 : theta = 1.058_f
110 : else
111 0 : theta = 1.058_f - (0.0155_f * (humidity - 0.97_f)) / (1.02_f - humidity**1.4_f)
112 : endif
113 :
114 0 : alpha1 = 1.2_f * exp((0.066_f * humidity) / (theta - humidity))
115 0 : f1 = 10.2_f - 23.7_f * humidity + 14.5_f * humidity**2
116 0 : f2 = -6.7_f + 15.5_f * humidity - 9.2_f * humidity**2
117 0 : alpharat = 1._f - f1 * (1._f - epsilon_) - f2 * (1._f - epsilon_**2)
118 :
119 : ! Scale the size based on the composition of the particle.
120 0 : select case(irhswcomp(igroup))
121 :
122 : case (I_SWF_NH42SO4)
123 0 : alphaComp = 1.00_f
124 :
125 : case(I_SWF_NH4NO3)
126 0 : alphaComp = 1.06_f
127 :
128 : case(I_SWF_NANO3)
129 0 : alphaComp = 1.17_f
130 :
131 : case(I_SWF_NH4CL)
132 0 : alphaComp = 1.23_f
133 :
134 : case(I_SWF_CACL2)
135 0 : alphaComp = 1.29_f
136 :
137 : case(I_SWF_NABR)
138 0 : alphaComp = 1.32_f
139 :
140 : case(I_SWF_NACL)
141 0 : alphaComp = 1.35_f
142 :
143 : case(I_SWF_MGCL2)
144 0 : alphaComp = 1.41_f
145 :
146 : case(I_SWF_LICL)
147 0 : alphaComp = 1.54_f
148 :
149 : case default
150 0 : if (do_print) write(LUNOPRT,*) "wetr:: ERROR - Unknown composition type (", irhswcomp(igroup), &
151 0 : ") for Fitzgerald."
152 0 : rc = RC_ERROR
153 0 : return
154 : end select
155 :
156 0 : alpha = alphaComp * (alpha1 * alpharat)
157 :
158 : ! Determine the wet radius.
159 : !
160 : ! NOTE: Fitgerald's equations assume r in [um], so scale the cgs units
161 : ! appropriately.
162 0 : rwet = (alpha * (rdry * 1e4_f)**beta) * (1e-4_f)
163 :
164 : ! Determine the wet density from the wet radius.
165 0 : r_ratio = (rdry / rwet)**3._f
166 0 : rhopwet = r_ratio * rhopdry + (1._f - r_ratio) * RHO_W
167 : end if
168 :
169 :
170 : ! Gerber Paremeterization
171 38135321593 : if (irhswell(igroup) == I_GERBER) then
172 :
173 : ! Scale the size based on the composition of the particle.
174 0 : select case(irhswcomp(igroup))
175 :
176 : case (I_SWG_NH42SO4)
177 : c1 = 0.4809_f
178 : c2 = 3.082_f
179 : c3 = 3.110e-11_f
180 0 : c4 = -1.428_f
181 :
182 : case(I_SWG_URBAN)
183 0 : c1 = 0.3926_f
184 0 : c2 = 3.101_f
185 0 : c3 = 4.190e-11_f
186 0 : c4 = -1.404_f
187 :
188 : case(I_SWG_RURAL)
189 0 : c1 = 0.2789_f
190 0 : c2 = 3.115_f
191 0 : c3 = 5.415e-11_f
192 0 : c4 = -1.399_f
193 :
194 : case(I_SWG_SEA_SALT)
195 0 : c1 = 0.7674_f
196 0 : c2 = 3.079_f
197 0 : c3 = 2.572e-11_f
198 0 : c4 = -1.424_f
199 :
200 : case default
201 0 : if (do_print) write(LUNOPRT,*) "wetr:: ERROR - Unknown composition type (", irhswcomp(igroup), &
202 0 : ") for Gerber."
203 0 : rc = RC_ERROR
204 0 : return
205 : end select
206 :
207 0 : rwet = ((c1 * rdry**c2 / (c3 * rdry**c4 - log10(humidity))) + rdry**3)**(1._f / 3._f)
208 :
209 : ! Determine the wet density from the wet radius.
210 0 : r_ratio = (rdry / rwet)**3
211 0 : rhopwet = r_ratio * rhopdry + (1._f - r_ratio) * RHO_W
212 :
213 : end if ! irhswell(igroup) == I_GERBER
214 :
215 : ! Mixed aerosol paremeterization (Pengfei Yu et al., JAMES, 2015) based on
216 : ! Petters and Kreidenweis (ACP, 2007) hygroscopicity parameter kappa
217 38135321593 : if (irhswell(igroup) == I_PETTERS) then
218 :
219 19508283799 : if (.not.( present(temp) .and. &
220 : present(kappa) )) then
221 0 : if (do_print) write(LUNOPRT,*) "wetr:: ERROR - h2o_mass,temp,kappa for PETTERS"
222 0 : rc = RC_ERROR
223 0 : return
224 : endif
225 :
226 19508283799 : if (temp .le. 190._f) then
227 1176090213 : rh190 = rh * pvap_h2o(temp) / pvap_h2o(190._f)
228 1176090213 : rh190 = min(max(rh190,tiny(1.0_f)), 0.995_f)
229 1176090213 : rwet = rdry * (1._f + rh190*kappa/(1._f-rh190))**(1._f/3._f)
230 : else ! temp > 190
231 18332193586 : rwet = rdry * (1._f + humidity*kappa/(1._f-humidity))**(1._f/3._f)
232 : end if
233 19508283799 : r_ratio = (rdry / rwet)**3._f
234 19508283799 : rhopwet = r_ratio * rhopdry + (1._f - r_ratio) * RHO_W
235 : end if ! irhswell(igroup) == I_PETTERS
236 :
237 :
238 : ! Sulfate Aerosol, using weight percent.
239 38135321593 : if (irhswell(igroup) == I_WTPCT_H2SO4) then
240 :
241 : ! Has the weight percent already been specified, or do we need to determine
242 : ! in based upon the water and temperature.
243 : !
244 : ! NOTE: This is used when generating optical properties files.
245 18627037794 : if (present(wgtpct) .and. present(temp)) then
246 0 : rhopwet = sulfate_density(carma, wgtpct, temp, rc)
247 0 : rwet = rdry * (100._f * rhopdry / wgtpct / rhopwet)**(1._f / 3._f)
248 18627037794 : else if (.not.( present(h2o_mass) .and. &
249 : present(h2o_vp) .and. &
250 : present(temp) )) then
251 0 : if (do_print) write(LUNOPRT,*) "wetr:: ERROR - h2o_mass,h2o_vp,temp for WTPCT_H2SO4"
252 0 : rc = RC_ERROR
253 0 : return
254 : else
255 : ! Adjust calculation for the Kelvin effect of H2O:
256 18627037794 : wtpkelv = 80._f ! start with assumption of 80 wt % H2SO4
257 18627037794 : den1 = 2.00151_f - 0.000974043_f * temp ! density at 79 wt %
258 18627037794 : den2 = 2.01703_f - 0.000988264_f * temp ! density at 80 wt %
259 18627037794 : drho_dwt = den2-den1 ! change in density for change in 1 wt %
260 :
261 18627037794 : sig1 = 79.3556_f - 0.0267212_f * temp ! surface tension at 79.432 wt %
262 18627037794 : sig2 = 75.608_f - 0.0269204_f * temp ! surface tension at 85.9195 wt %
263 18627037794 : dsigma_dwt = (sig2-sig1) / (85.9195_f - 79.432_f) ! change in density for change in 1 wt %
264 18627037794 : sigkelv = sig1 + dsigma_dwt * (80.0_f - 79.432_f)
265 :
266 18627037794 : rwet = rdry * (100._f * rhopdry / wtpkelv / den2)**(1._f / 3._f)
267 :
268 : rkelvinH2O_b = 1._f + wtpkelv * drho_dwt / den2 - 3._f * wtpkelv &
269 18627037794 : * dsigma_dwt / (2._f*sigkelv)
270 :
271 18627037794 : rkelvinH2O_a = 2._f * gwtmol(igash2so4) * sigkelv / (den1 * RGAS * temp * rwet)
272 :
273 18627037794 : rkelvinH2O = exp (rkelvinH2O_a*rkelvinH2O_b)
274 :
275 18627037794 : h2o_kelv = h2o_mass / rkelvinH2O
276 18627037794 : wtpkelv = wtpct_tabaz(carma, temp, h2o_kelv, h2o_vp, rc)
277 18627037794 : rhopwet = sulfate_density(carma, wtpkelv, temp, rc)
278 18627037794 : rwet = rdry * (100._f * rhopdry / wtpkelv / rhopwet)**(1._f / 3._f)
279 : endif
280 : end if ! irhswell(igroup) == I_WTPCT_H2SO4
281 : end if ! irhswell(igroup) /= I_NO_SWELLING
282 :
283 : ! Return to caller with wet radius evaluated.
284 38135321593 : return
285 :
286 : contains
287 :
288 : !! This function is needed for generating wet radii for optics when using the
289 : !! PETTERS scheme, and should not be used generally within CARMA.
290 : !!
291 : !! The vaporp_h2o_murphy2005 equation to calculate the vapor pressure at 190 K
292 : !! for liquid water
293 2352180426 : pure function pvap_h2o(temp) result(pvap)
294 :
295 : real(kind=f), intent(in) :: temp
296 : real(kind=f) :: pvap
297 :
298 : pvap = temp / ( 10.0_f * exp(54.842763_f - (6763.22_f / temp) &
299 : - (4.210_f * log(temp)) + (0.000367_f * temp) + (tanh(0.0415_f * (temp - 218.8_f)) &
300 2352180426 : * (53.878_f - (1331.22_f / temp) - (9.44523_f * log(temp)) + 0.014025_f * temp))) )
301 :
302 38135321593 : end function pvap_h2o
303 :
304 : end subroutine getwetr
305 :
306 : end module wetr
307 :
|
