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 fall velocities, vf(k) [cm s^-1]
6 : !! and reynolds' numbers based on fall velocities, re(j,i,k) [dimensionless].
7 : !! indices correspond to vertical level <k>, bin index <i>, and aerosol
8 : !! group <j>.
9 : !!
10 : !! Non-spherical particles are treated through shape factors <ishape>
11 : !! and <eshape>.
12 : !!
13 : !! General method is to first use Stokes' flow to estimate fall
14 : !!! velocity, then calculate reynolds' number, then use "y function"
15 : !! (defined in Pruppacher and Klett) to reevaluate reynolds' number,
16 : !! from which the fall velocity is finally obtained.
17 : !!
18 : !! This routine requires that vertical profiles of temperature <t>,
19 : !! air density <rhoa>, and viscosity <rmu> are defined (i.e., initatm.f
20 : !! must be called before this).
21 : !!
22 : !! @author Chuck Bardeen, Pete Colarco from Andy Ackerman
23 : !! @version Mar-2010 from Oct-1995
24 :
25 :
26 0 : subroutine setupvf_std_shape(carma, cstate, j, rc)
27 :
28 : ! types
29 : use carma_precision_mod
30 : use carma_enums_mod
31 : use carma_constants_mod
32 : use carma_types_mod
33 : use carmastate_mod
34 : use carma_mod
35 :
36 : implicit none
37 :
38 : type(carma_type), intent(in) :: carma !! the carma object
39 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
40 : integer, intent(in) :: j !! group index
41 : integer, intent(inout) :: rc !! return code, negative indicates failure
42 :
43 : ! Local declarations
44 : integer :: i, k, ilast
45 : real(kind=f) :: x, y
46 : real(kind=f) :: rhoa_cgs, vg, rmfp, rkn, expon
47 : real(kind=f) :: f1, f2, f3, ex, exx, exy, xcc, xa, bxx, r_shape, rfix, b0, bb1, bb2, bb3, z
48 :
49 : ! Define formats
50 : 1 format('setupvfall::ERROR - ishape != 1, no fall velocity algorithm')
51 :
52 :
53 : ! First evaluate factors that depend upon particle shape (used in correction
54 : ! factor <bpm> below).
55 0 : if (ishape(j) .eq. I_SPHERE) then
56 :
57 : ! Spheres
58 : f1 = 1.0_f
59 : f2 = 1.0_f
60 :
61 0 : else if (ishape(j) .eq. I_HEXAGON) then
62 :
63 : ! Hexagons: taken from Turco et al (Planet. Space Sci. Rev. 30, 1147-1181, 1982)
64 : ! with diffuse reflection of air molecules assumed
65 0 : f2 = (PI / 9._f / tan(PI / 6._f))**(ONE/3._f) * eshape(j)**(ONE/6._f)
66 :
67 0 : else if (ishape(j) .eq. I_CYLINDER)then
68 :
69 : ! Spheroids: also from Turco et al. [1982]
70 0 : f2 = (2._f / 3._f)**(ONE/3._f) * eshape(j)**(ONE/6._f)
71 : endif
72 :
73 : ! (following statement yields <f3> = 1.0 for <eshape> = I_SPHERE)
74 0 : f3 = 1.39_f / sqrt((1.14_f + 0.25_f / eshape(j)) * (0.89_f + eshape(j) / 2._f))
75 0 : f2 = f2 * f3
76 :
77 0 : if (eshape(j) .gt. 1._f) then
78 :
79 : ! For Stokes regime there is no separate data for hexagonal plates or columns,
80 : ! so we use prolate spheroids. This is from Fuchs' book.
81 0 : exx = eshape(j)**2 - 1._f
82 0 : exy = sqrt(exx)
83 0 : xcc = 1.333_f * exx / ((2._f * eshape(j)**2 - 1._f) * log(eshape(j) + exy) / exy-eshape(j))
84 0 : xa = 2.666_f * exx / ((2._f * eshape(j)**2 - 3._f) * log(eshape(j) + exy) / exy+eshape(j))
85 : ! f1 = eshape(j)**(-ONE/3._f) * (xcc + 2._f*xa) / 3._f
86 0 : f1 = eshape(j)**(-2._f/3._f) * (xcc + 2._f*xa) / 3._f
87 :
88 0 : elseif (eshape(j) .lt. 1._f) then
89 :
90 : ! Use oblate spheroids for disks (eshape < 1.). Also from Fuchs' book.
91 0 : bxx = 1._f / eshape(j)
92 0 : exx = bxx**2 - 1._f
93 0 : exy = sqrt(exx)
94 0 : xcc = 1.333_f * exx / (bxx * (bxx**2 - 2._f) * atan(exy) / exy + bxx)
95 0 : xa = 2.666_f * exx / (bxx * (3._f * bxx**2 - 2._f) * atan(exy) / exy - bxx)
96 0 : f1 = bxx**(ONE/3._f) * (xcc + 2._f * xa) / 3._f
97 : endif
98 :
99 :
100 : ! Loop over column with ixy = 1
101 0 : do k = 1,NZ
102 :
103 : ! This is <rhoa> in cartesian coordinates (good old cgs units)
104 0 : rhoa_cgs = rhoa(k) / zmet(k)
105 :
106 : ! <vg> is mean thermal velocity of air molecules [cm/s]
107 0 : vg = sqrt(8._f / PI * R_AIR * t(k))
108 :
109 : ! <rmfp> is mean free path of air molecules [cm]
110 0 : rmfp = 2._f * rmu(k) / (rhoa_cgs * vg)
111 :
112 : ! Loop over particle size bins.
113 0 : do i = 1,NBIN
114 :
115 : ! <r_shape> is radius of particle used to calculate <re>.
116 0 : if (ishape(j) .eq. I_SPHERE) then
117 0 : r_shape = r_wet(k,i,j)
118 0 : else if (ishape(j) .eq. I_HEXAGON) then
119 0 : r_shape = r_wet(k,i,j) * 0.8456_f * eshape(j)**(-ONE/3._f)
120 0 : else if(ishape(j) .eq. I_CYLINDER) then
121 : ! r_shape = r_wet(k,i,j) * eshape(j)**(-ONE/3._f)
122 :
123 : ! Shouldn't this have a factor related to being a cylinder vs a
124 : ! sphere in addition to the aspect ratio factor?
125 0 : r_shape = r_wet(k,i,j) * 0.8736_f * eshape(j)**(-ONE/3._f)
126 : endif
127 :
128 : ! <rkn> is knudsen number
129 0 : rkn = rmfp / r_wet(k,i,j)
130 :
131 : ! <bpm> is the slip correction factor, the correction term for
132 : ! non-continuum effects. Also used to calculate coagulation kernels
133 : ! and diffusion coefficients.
134 0 : expon = -.87_f / rkn
135 0 : expon = max(-POWMAX, expon)
136 0 : bpm(k,i,j) = 1._f + f1*f2*(1.246_f*rkn + 0.42_f*rkn*exp(expon))
137 :
138 : ! These are first guesses for fall velocity and Reynolds' number,
139 : ! valid for Reynolds' number < 0.01
140 : !
141 : ! This is "regime 1" in Pruppacher and Klett (chap. 10, pg 416).
142 0 : vf(k,i,j) = (2._f / 9._f) * rhop_wet(k,i,j) *(r_wet(k,i,j)**2) * GRAV * bpm(k,i,j) / (f1 * rmu(k))
143 0 : re(k,i,j) = 2._f * rhoa_cgs * r_shape * vf(k,i,j) / rmu(k)
144 :
145 :
146 : ! <rfix> is used in drag coefficient.
147 0 : rfix = vol(i,j) * rhop_wet(k,i,j) * GRAV * rhoa_cgs / rmu(k)**2
148 :
149 0 : if ((re(k,i,j) .ge. 0.01_f) .and. (re(k,i,j) .le. 300._f)) then
150 :
151 : ! This is "regime 2" in Pruppacher and Klett (chap. 10, pg 417).
152 : !
153 : ! NOTE: This sphere case is not the same solution used when
154 : ! interpolating other shape factors. This seems potentially inconsistent.
155 0 : if (ishape(j) .eq. I_SPHERE) then
156 :
157 0 : x = log(24._f * re(k,i,j) / bpm(k,i,j))
158 : y = -0.3318657e1_f + x * 0.992696_f - x**2 * 0.153193e-2_f - &
159 : x**3 * 0.987059e-3_f - x**4 * 0.578878e-3_f + &
160 0 : x**5 * 0.855176E-04_f - x**6 * 0.327815E-05_f
161 :
162 0 : if (y .lt. -675._f) y = -675._f
163 0 : if (y .ge. 741._f) y = 741._f
164 :
165 0 : re(k,i,j) = exp(y) * bpm(k,i,j)
166 :
167 0 : else if (eshape(j) .le. 1._f) then
168 :
169 : ! P&K pg. 427
170 0 : if (ishape(j) .eq. I_HEXAGON) then
171 0 : x = log10(16._f * rfix / (3._f * sqrt(3._f)))
172 0 : else if (ishape(j) .eq. I_CYLINDER) then
173 0 : x = log10(8._f * rfix / PI)
174 : endif
175 :
176 0 : if (eshape(j) .le. 0.2_f) then
177 :
178 : ! P&K, page 424-427
179 : b0 = -1.33_f
180 : bb1 = 1.0217_f
181 : bb2 = -0.049018_f
182 : bb3 = 0.0_f
183 0 : else if (eshape(j) .le. 0.5_f) then
184 :
185 : ! NOTE: This interpolation/extrapolation method is
186 : ! not discussed in P&K; although, the solution for
187 : ! eshape = 0.5 is shown. Does this really work?
188 0 : ex = (eshape(j) - 0.2_f) / 0.3_f
189 0 : b0 = -1.33_f + ex * (-1.3247_f + 1.33_f)
190 0 : bb1 = 1.0217_f + ex * (1.0396_f - 1.0217_f)
191 0 : bb2 = -0.049018_f + ex * (-0.047556_f + 0.049018_f)
192 0 : bb3 = ex * (-0.002327_f)
193 : else
194 :
195 : ! Extrapolating to cylinder cases on 436.
196 0 : ex = (eshape(j) - 0.5_f) / 0.5_f
197 0 : b0 = -1.3247_f + ex * (-1.310_f + 1.3247_f)
198 0 : bb1 = 1.0396_f + ex * (0.98968_f - 1.0396_f)
199 0 : bb2 = -0.047556_f + ex * (-0.042379_f + 0.047556_f)
200 0 : bb3 = -0.002327_f + ex * ( 0.002327_f)
201 : endif
202 :
203 0 : y = b0 + x * bb1 + x**2 * bb2 + x**3 * bb3
204 0 : re(k,i,j) = 10._f**y * bpm(k,i,j)
205 :
206 0 : else if (eshape(j) .gt. 1._f) then
207 : ! Why is this so different from the oblate case?
208 : ! This seems wrong.
209 : ! x = log10(2._f * rfix / eshape(j))
210 0 : if (ishape(j) .eq. I_CYLINDER) then
211 0 : x = log10(8._f * rfix / PI)
212 : endif
213 :
214 : ! P&K pg 430
215 0 : if( eshape(j) .le. 2._f )then
216 0 : ex = eshape(j) - 1._f
217 0 : b0 = -1.310_f + ex * (-1.11812_f + 1.310_f)
218 0 : bb1 = 0.98968_f + ex * (0.97084_f - 0.98968_f)
219 0 : bb2 = -0.042379_f + ex * (-0.058810_f + 0.042379_f)
220 0 : bb3 = ex * (0.002159_f)
221 0 : else if (eshape(j) .le. 10._f) then
222 0 : ex = (eshape(j) - 2._f) / 8.0_f
223 0 : b0 = -1.11812_f + ex * (-0.90629_f + 1.11812_f)
224 0 : bb1 = 0.97084_f + ex * (0.90412_f - 0.97084_f)
225 0 : bb2 = -0.058810_f + ex * (-0.059312_f + 0.058810_f)
226 0 : bb3 = 0.002159_f + ex * (0.0029941_f - 0.002159_f)
227 : else
228 :
229 : ! This is interpolating to a solution for an infinite
230 : ! cylinder, so it may not be the greatest estimate.
231 0 : ex = 10._f / eshape(j)
232 0 : b0 = -0.79888_f + ex * (-0.90629_f + 0.79888_f)
233 0 : bb1 = 0.80817_f + ex * (0.90412_f - 0.80817_f)
234 0 : bb2 = -0.030528_f + ex * (-0.059312_f + 0.030528_f)
235 0 : bb3 = ex * (0.0029941_f)
236 : endif
237 :
238 0 : y = b0 + x * bb1 + x**2 * bb2 + x**3 * bb3
239 0 : re(k,i,j) = 10._f**y * bpm(k,i,j)
240 :
241 : endif
242 :
243 : ! Adjust <vf> for non-sphericicity.
244 0 : vf(k,i,j) = re(k,i,j) * rmu(k) / (2._f * r_shape * rhoa_cgs)
245 :
246 : endif
247 :
248 0 : if (re(k,i,j) .gt. 300._f) then
249 :
250 : ! This is "regime 3" in Pruppacher and Klett (chap. 10, pg 418).
251 :
252 : ! if ((do_print) .and. (ishape(j) .ne. I_SPHERE)) write(LUNOPRT,1)
253 : ! if ((do_print) .and. (ishape(j) .ne. I_SPHERE)) write(LUNOPRT,*) "setupvfall:", j, i, k, re(k,i,j)
254 : ! rc = RC_ERROR
255 : ! return
256 :
257 0 : z = ((1.e6_f * rhoa_cgs**2) / (GRAV * rhop_wet(k,i,j) * rmu(k)**4))**(ONE/6._f)
258 0 : b0 = (24._f * vf(k,i,j) * rmu(k)) / 100._f
259 0 : x = log(z * b0)
260 : y = -5.00015_f + x * (5.23778_f - x * (2.04914_f - x * (0.475294_f - &
261 0 : x * (0.0542819_f - x * 0.00238449_f))))
262 :
263 0 : if (y .lt. -675._f) y = -675.0_f
264 0 : if (y .ge. 741._f) y = 741.0_f
265 :
266 0 : re(k,i,j) = z * exp(y) * bpm(k,i,j)
267 0 : vf(k,i,j) = re(k,i,j) * rmu(k) / ( 2._f * r_wet(k,i,j) * rhoa_cgs)
268 :
269 : ! Values should not decrease with diameter, but instead should
270 : ! reach a limiting velocity that is independent of size (see
271 : ! Figure 10-25 of Pruppacher and Klett, 1997)
272 0 : ilast = max(1,i-1)
273 0 : if ((vf(k,i,j) .lt. vf(k,ilast,j)) .or. (re(k,i,j) .gt. 4000._f)) then
274 0 : vf(k,i,j) = vf(k,ilast,j)
275 : endif
276 : endif
277 : enddo ! <i=1,NBIN>
278 : enddo ! <k=1,NZ>
279 :
280 : ! Return to caller with particle fall velocities evaluated.
281 0 : return
282 0 : end
|
