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 defines radius-dependent but time-independent parameters
6 : !! used to calculate condensational growth of particles. Growth rates
7 : !! are calculated at bin boundaries: the parameters calculated here
8 : !! ( <gro>, <gro1>, <gro2>, and <akelvin> )
9 : !! are defined at lower bin boundaries through the growth rate expression
10 : !! (for one particle) used in growevapl.f:
11 : !!>
12 : !! dm = gro*pvap*( S + 1 - Ak*As - gro1*gro2*qrad )
13 : !! -- -------------------------------------------
14 : !! dt 1 + gro*gro1*pvap
15 : !!
16 : !! where
17 : !!
18 : !! S = supersaturation
19 : !! Ak = exp(akelvin/r)
20 : !! As = exp(-sol_ions * solute_mass/solwtmol * gwtmol/condensate_mass)
21 : !! pvap = saturation vapor pressure [dyne cm**-2]
22 : !! qrad = radiative energy absorbed
23 : !!<
24 : !! This routine requires that vertical profiles of temperature <T>,
25 : !! and pressure <p> are defined.
26 : !!
27 : !! This routine also requires that particle Reynolds' numbers are
28 : !! defined (setupvfall.f must be called before this).
29 : !!
30 : !! @author Andy Ackerman
31 : !! @version Dec-1995
32 1050624 : subroutine setupgkern(carma, cstate, rc)
33 :
34 : ! types
35 : use carma_precision_mod
36 : use carma_enums_mod
37 : use carma_constants_mod
38 : use carma_types_mod
39 : use carmastate_mod
40 : use carma_mod
41 : use sulfate_utils
42 :
43 : implicit none
44 :
45 : type(carma_type), intent(in) :: carma !! the carma object
46 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
47 : integer, intent(inout) :: rc !! return code, negative indicates failure
48 :
49 : ! Local declarations
50 : integer :: igas !! gas index
51 : integer :: ielem !! element index
52 : integer :: k !! z index
53 : integer :: igroup !! group index
54 : integer :: i
55 : real(kind=f) :: gstick
56 : real(kind=f) :: cor
57 : real(kind=f) :: phish
58 : real(kind=f) :: esh1
59 : real(kind=f) :: a1
60 : real(kind=f) :: br
61 : real(kind=f) :: rknudn
62 : real(kind=f) :: rknudnt
63 : real(kind=f) :: rlam
64 : real(kind=f) :: rlamt
65 2101248 : real(kind=f) :: rhoa_cgs(NZ, NGAS)
66 2101248 : real(kind=f) :: freep(NZ, NGAS)
67 2101248 : real(kind=f) :: freept(NZ, NGAS)
68 : real(kind=f) :: rlh
69 : real(kind=f) :: diffus1
70 : real(kind=f) :: thcond1
71 : real(kind=f) :: reyn_shape
72 : real(kind=f) :: schn
73 : real(kind=f) :: prnum
74 : real(kind=f) :: x1
75 : real(kind=f) :: x2
76 : real(kind=f) :: fv
77 : real(kind=f) :: surf_tens ! surface tension of H2SO4 particle
78 : real(kind=f) :: rho_H2SO4 ! wet density of H2SO4 particle
79 :
80 :
81 : ! Calculate gas properties for all of the gases. Better to do them all once, than to
82 : ! repeat this for multiple groups.
83 3151872 : do igas = 1, NGAS
84 :
85 : ! Radius-independent parameters for condensing gas
86 : !
87 : ! This is <rhoa> in cgs units.
88 : !
89 149188608 : rhoa_cgs(:, igas) = rhoa(:) / zmet(:)
90 :
91 2101248 : if (igas .eq. igash2o) then
92 :
93 : ! Condensing gas is water vapor
94 : !
95 : ! <surfctwa> is surface tension of water-air interface (valid from 0 to 40 C)
96 : ! from Pruppacher and Klett (eq. 5-12).
97 74594304 : surfctwa(:) = 76.10_f - 0.155_f*( t(:) - 273.16_f )
98 :
99 : ! <surfctiw> is surface tension of water-ice interface
100 : ! from Pruppacher and Klett (eq. 5-48).!
101 74594304 : surfctiw(:) = 28.5_f + 0.25_f*( t(:) - 273.16_f )
102 :
103 : ! <surfctiw> is surface tension of water-ice interface
104 : ! from Hale and Plummer [J. Chem. Phys., 61, 1974].
105 74594304 : surfctia(:) = 141._f - 0.15_f * t(:)
106 :
107 : ! <akelvin> is argument of exponential in kelvin curvature term.
108 0 : akelvin(:,igas) = 2._f*gwtmol(igas)*surfctwa(:) &
109 74594304 : / ( t(:)*RHO_W*RGAS )
110 :
111 0 : akelvini(:,igas) = 2._f*gwtmol(igas)*surfctia(:) &
112 74594304 : / ( t(:)*RHO_W*RGAS )
113 :
114 : ! condensing gas is H2SO4
115 1050624 : else if (igas .eq. igash2so4) then
116 :
117 : ! Calculate Kelvin curvature factor for H2SO4 interactively with temperature:
118 74594304 : do k = 1, NZ
119 73543680 : surf_tens = sulfate_surf_tens(carma, wtpct(k), t(k), rc)
120 73543680 : rho_H2SO4 = sulfate_density(carma, wtpct(k), t(k), rc)
121 73543680 : akelvin(k, igas) = 2._f * gwtmol(igas) * surf_tens / (t(k) * rho_H2SO4 * RGAS)
122 :
123 : ! Not doing condensation of h2So4 on ice, so just set it to the value
124 : ! for water vapor.
125 74594304 : akelvini(k, igas) = akelvini(k, igash2o)
126 : end do
127 : else
128 :
129 : ! Condensing gas is not yet configured.
130 0 : if (do_print) write(LUNOPRT,*) 'setupgkern::ERROR - invalid igas'
131 0 : rc = RC_ERROR
132 0 : return
133 : endif
134 :
135 : ! Molecular free path of condensing gas
136 2101248 : freep(:,igas) = 3._f*diffus(:,igas) &
137 151289856 : * sqrt( ( PI*gwtmol(igas) ) / ( 8._f*RGAS*t(:) ) )
138 :
139 : ! Thermal free path of condensing gas
140 0 : freept(:,igas) = freep(:,igas)*thcond(:) / &
141 0 : ( diffus(:,igas) * rhoa_cgs(:, igas) &
142 150239232 : * ( CP - RGAS/( 2._f*WTMOL_AIR ) ) )
143 : end do
144 :
145 :
146 : ! Loop over aerosol groups only (no radius, gas, or spatial dependence).
147 3151872 : do igroup = 1, NGROUP
148 :
149 : ! Use gstickl or gsticki, depending on whether group is ice or not
150 2101248 : if( is_grp_ice(igroup) ) then
151 0 : gstick = gsticki
152 : else
153 2101248 : gstick = gstickl
154 : endif
155 :
156 : ! Non-spherical corrections (need a reference for these)
157 2101248 : if( ishape(igroup) .eq. I_SPHERE )then
158 :
159 : ! Spheres
160 : cor = 1._f
161 : phish = 1._f
162 : else
163 :
164 0 : if( ishape(igroup) .eq. I_HEXAGON )then
165 :
166 : ! Hexagons
167 : phish = 6._f/PI*tan(PI/6._f)*( eshape(igroup) + 0.5_f ) &
168 0 : * ( PI / ( 9._f*eshape(igroup)*tan(PI/6._f) ) )**(2._f/3._f)
169 :
170 0 : else if( ishape(igroup) .eq. I_CYLINDER )then
171 :
172 : ! Spheroids
173 : phish = ( eshape(igroup) + 0.5_f ) &
174 0 : * ( 2._f / ( 3._f*eshape(igroup) ) )**(2._f/3._f)
175 : endif
176 :
177 0 : if( eshape(igroup) .lt. 1._f )then
178 :
179 : ! Oblate spheroids
180 0 : esh1 = 1._f / eshape(igroup)
181 0 : a1 = sqrt(esh1**2 - 1._f)
182 0 : cor = a1 / asin( a1 / esh1 ) / esh1**(2._f/3._f)
183 : else
184 :
185 : ! Prolate spheroids
186 0 : a1 = sqrt( eshape(igroup)**2 - 1._f )
187 : cor = a1 / log( eshape(igroup) + a1 ) &
188 0 : / eshape(igroup)**(ONE/3._f)
189 : endif
190 : endif
191 :
192 : ! Evaluate growth terms only for particle elements that grow.
193 : ! particle number concentration element
194 2101248 : ielem = ienconc(igroup)
195 :
196 : ! condensing gas is <igas>
197 2101248 : igas = igrowgas(ielem)
198 :
199 : ! If the group doesn't grow, but is involved in aerosol
200 : ! freezing, then the gas properties still need to be calculated.
201 2101248 : if( igas .eq. 0 ) igas = inucgas(igroup)
202 :
203 1050624 : if( igas .ne. 0 )then
204 :
205 149188608 : do k = 1, NZ
206 :
207 : ! Latent heat of condensing gas
208 147087360 : if( is_grp_ice(igroup) )then
209 0 : rlh = rlhe(k,igas) + rlhm(k,igas)
210 : else
211 147087360 : rlh = rlhe(k,igas)
212 : endif
213 :
214 : ! Radius-dependent parameters
215 3090935808 : do i = 1, NBIN
216 :
217 2941747200 : br = rlow_wet(k,i,igroup) ! particle bin Boundary Radius
218 :
219 : ! These are Knudsen numbers
220 2941747200 : rknudn = freep(k,igas) / br
221 2941747200 : rknudnt = freept(k,igas) / br
222 :
223 : ! These are "lambdas" used in correction for gas kinetic effects.
224 : rlam = ( 1.33_f*rknudn + 0.71_f ) / ( rknudn + 1._f ) &
225 2941747200 : + ( 4._f*( 1._f - gstick ) ) / ( 3._f*gstick )
226 :
227 : rlamt = ( 1.33_f*rknudnt + 0.71_f ) / ( rknudnt + 1._f ) &
228 2941747200 : + ( 4._f*( 1._f - tstick ) ) / ( 3._f*tstick )
229 :
230 : ! Diffusion coefficient and thermal conductivity modified for
231 : ! free molecular limit and for particle shape.
232 2941747200 : diffus1 = diffus(k,igas)*cor / ( 1._f + rlam*rknudn*cor/phish )
233 2941747200 : thcond1 = thcond(k)*cor / ( 1._f + rlamt*rknudnt*cor/phish )
234 :
235 : ! Save the modified thermal conductivity off so it can be used in pheat.
236 2941747200 : thcondnc(k,i,igroup) = thcond1
237 :
238 : ! Reynolds' number based on particle shape <reyn_shape>
239 2941747200 : if( ishape(igroup) .eq. I_SPHERE )then
240 2941747200 : reyn_shape = re(k,i,igroup)
241 :
242 0 : else if( eshape(igroup) .lt. 1._f )then
243 0 : reyn_shape = re(k,i,igroup) * ( 1._f + 2._f*eshape(igroup) )
244 :
245 : else
246 0 : reyn_shape = re(k,i,igroup) * PI*( 1._f+2._f*eshape(igroup) ) &
247 0 : / ( 2._f*( 1._f + eshape(igroup) ) )
248 : endif
249 :
250 : ! Particle Schmidt number
251 2941747200 : schn = rmu(k) / ( rhoa_cgs(k,igas) * diffus1 )
252 :
253 : ! Prandtl number
254 2941747200 : prnum = rmu(k)*CP/thcond1
255 :
256 : ! Ventilation factors <fv> and <ft> from Pruppacher and Klett
257 2941747200 : x1 = schn **(ONE/3._f) * sqrt( reyn_shape )
258 2941747200 : x2 = prnum**(ONE/3._f) * sqrt( reyn_shape )
259 :
260 2941747200 : if( is_grp_ice(igroup) )then
261 :
262 : ! Ice crystals
263 0 : if( x1 .le. 1._f )then
264 0 : fv = 1._f + 0.14_f*x1**2
265 : else
266 0 : fv = 0.86_f + 0.28_f*x1
267 : endif
268 :
269 0 : if( x2 .le. 1._f )then
270 0 : ft(k,i,igroup) = 1._f + 0.14_f*x2**2
271 : else
272 0 : ft(k,i,igroup) = 0.86_f + 0.28_f*x2
273 : endif
274 : else
275 :
276 : ! Liquid water drops
277 2941747200 : if( x1 .le. 1.4_f )then
278 2931957631 : fv = 1._f + 0.108_f*x1**2
279 : else
280 9789569 : fv = 0.78_f + 0.308_f*x1
281 : endif
282 :
283 2941747200 : if( x2 .le. 1.4_f )then
284 2930486108 : ft(k,i,igroup) = 1._f + 0.108_f*x2**2
285 : else
286 11261092 : ft(k,i,igroup) = 0.78_f + 0.308_f*x2
287 : endif
288 : endif
289 :
290 : ! Growth kernel for particle without radiation or heat conduction at
291 : ! radius lower boundary [g cm^3 / erg / s]
292 0 : gro(k,i,igroup) = 4._f*PI*br &
293 0 : * diffus1*fv*gwtmol(igas) &
294 2941747200 : / ( BK*t(k)*AVG )
295 :
296 : ! Coefficient for conduction term in growth kernel [s/g]
297 0 : gro1(k,i,igroup) = gwtmol(igas)*rlh**2 &
298 0 : / ( RGAS*t(k)**2*ft(k,i,igroup)*thcond1 ) &
299 2941747200 : / ( 4._f*PI*br )
300 :
301 : ! Coefficient for radiation term in growth kernel [g/erg]
302 : ! (note: no radial dependence).
303 3088834560 : if( i .eq. 1 )then
304 147087360 : gro2(k,igroup) = 1._f / rlh
305 : endif
306 :
307 : enddo ! i=1,NBIN
308 : enddo ! k=1,NZ
309 : endif ! igas ne 0
310 : enddo ! igroup=1,NGROUP
311 :
312 : ! Return to caller with time-independent particle growth
313 : ! parameters initialized.
314 : return
315 1050624 : end
|
