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 : !! Method: first use Stokes flow (with Fuchs' size corrections,
11 : !! valid only for Stokes flow) to estimate fall velocity, then calculate
12 : !! Reynolds' number (Re) (for spheres, Stokes drag coefficient is 24/Re).
13 : !! Then for Re > 1, correct drag coefficient (Cd) for turbulent boundary
14 : !! layer through standard trick to solving the drag problem:
15 : !! fit y = log( Re ) as a function of x = log( Cd Re^2 ).
16 : !! We use the data for rigid spheres taken from Figure 10-6 of
17 : !! Pruppacher and Klett (1978):
18 : !!
19 : !! Re Cd
20 : !! ----- ------
21 : !! 1 24
22 : !! 10 4.3
23 : !! 100 1.1
24 : !! 1000 0.45
25 : !!
26 : !! Note that we ignore the "drag crisis" at Re > 200,000
27 : !! (as discussed on p. 341 and shown in Fig 10-36 of P&K 1978), where
28 : !! Cd drops dramatically to 0.2 for smooth, rigid spheres, and instead
29 : !! assume Cd = 0.45 for Re > 1,000
30 : !!
31 : !! Note that we also ignore hydrodynamic deformation of liquid droplets
32 : !! as well as any breakup due to Rayleigh-Taylor instability.
33 : !!
34 : !! This routine requires that vertical profiles of temperature <t>,
35 : !! air density <rhoa>, and viscosity <rmu> are defined (i.e., initatm.f
36 : !! must be called before this). The vertical profile with ix = iy = 1
37 : !! is used.
38 : !!
39 : !! We assume spherical particles -- call setupvf_std_shape() to use legacy
40 : !! code from old Toon model for non-spherical effects -- use (better
41 : !! yet, fix) at own risk.
42 : !!
43 : !! Added support for the particle radius being dependent on the relative
44 : !! humidity according to the parameterizations of Gerber [1995] and
45 : !! Fitzgerald [1975]. The fall velocity is then based upon the wet radius
46 : !! rather than the dry radius. For particles that are not subject to
47 : !! swelling, the wet and dry radii are the same.
48 : !!
49 : !! @author Chuck Bardeen, Pete Colarco from Andy Ackerman
50 : !! @version Mar-2010 from Nov-2000
51 2101248 : subroutine setupvf_std(carma, cstate, j, rc)
52 :
53 : ! types
54 : use carma_precision_mod
55 : use carma_enums_mod
56 : use carma_constants_mod
57 : use carma_types_mod
58 : use carmastate_mod
59 : use carma_mod
60 :
61 : implicit none
62 :
63 : type(carma_type), intent(in) :: carma !! the carma object
64 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
65 : integer, intent(in) :: j !! group index
66 : integer, intent(inout) :: rc !! return code, negative indicates failure
67 :
68 : ! Local declarations
69 : integer :: i, k
70 : real(kind=f) :: x, y, cdrag
71 : real(kind=f) :: rhoa_cgs, vg, rmfp, rkn, expon
72 :
73 : ! Define formats
74 : 1 format(/,'Non-spherical particles specified for group ',i3, &
75 : ' (ishape=',i3,') but spheres assumed in I_FALLRTN_STD.', &
76 : ' Suggest using non-spherical code in I_FALLRTN_STD_SHAPE.')
77 :
78 : ! Warning message for non-spherical particles!
79 2101248 : if( ishape(j) .ne. 1 )then
80 0 : if (do_print) write(LUNOPRT,1) j, ishape(j)
81 : endif
82 :
83 : ! Loop over all atltitudes.
84 149188608 : do k = 1, NZ
85 :
86 : ! This is <rhoa> in cartesian coordinates (good old cgs units)
87 147087360 : rhoa_cgs = rhoa(k) / zmet(k)
88 :
89 : ! <vg> is mean thermal velocity of air molecules [cm/s]
90 147087360 : vg = sqrt(8._f / PI * R_AIR * t(k))
91 :
92 : ! <rmfp> is mean free path of air molecules [cm]
93 147087360 : rmfp = 2._f * rmu(k) / (rhoa_cgs * vg)
94 :
95 : ! Loop over particle size bins.
96 3090935808 : do i = 1,NBIN
97 :
98 : ! <rkn> is knudsen number
99 2941747200 : rkn = rmfp / (r_wet(k,i,j) * rrat(i,j))
100 :
101 : ! <bpm> is the slip correction factor, the correction term for
102 : ! non-continuum effects. Also used to calculate coagulation kernels
103 : ! and diffusion coefficients.
104 2941747200 : expon = -.87_f / rkn
105 2941747200 : expon = max(-POWMAX, expon)
106 2941747200 : bpm(k,i,j) = 1._f + (1.246_f*rkn + 0.42_f*rkn*exp(expon))
107 :
108 : ! Stokes fall velocity and Reynolds' number
109 2941747200 : vf(k,i,j) = (ONE * 2._f / 9._f) * rhop_wet(k,i,j) * r_wet(k,i,j)**2 * GRAV * bpm(k,i,j) / rmu(k) / rprat(i,j)
110 2941747200 : re(k,i,j) = 2._f * rhoa_cgs * r_wet(k,i,j) * rprat(i,j) * vf(k,i,j) / rmu(k)
111 :
112 3088834560 : if (re(k,i,j) .ge. 1._f) then
113 :
114 : ! Correct drag coefficient for turbulence
115 1480198 : x = log(re(k,i,j) / bpm(k,i,j))
116 1480198 : y = x*(0.83_f - 0.013_f*x)
117 :
118 1480198 : re(k,i,j) = exp(y) * bpm(k,i,j)
119 :
120 1480198 : if (re(k,i,j) .le. 1.e3_f) then
121 :
122 : ! drag coefficient from quadratic fit y(x) when Re < 1,000
123 1480198 : vf(k,i,j) = re(k,i,j) * rmu(k) / (2._f * r_wet(k,i,j) * rprat(i,j) * rhoa_cgs)
124 : else
125 :
126 : ! drag coefficient = 0.45 independent of Reynolds number when Re > 1,000
127 0 : cdrag = 0.45_f
128 0 : vf(k,i,j) = bpm(k,i,j) * &
129 0 : sqrt( 8._f * rhop_wet(k,i,j) * r_wet(k,i,j) * GRAV / &
130 0 : (3._f * cdrag * rhoa_cgs * rprat(i,j)**2._f) )
131 : endif
132 : endif
133 : enddo ! <i=1,NBIN>
134 : enddo ! <k=1,NZ>
135 :
136 : ! Return to caller with particle fall velocities evaluated.
137 2101248 : return
138 2101248 : end
|
