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 the coagulation kernels, ckernel(k,j1,j2,i1,i2)
6 : !! [cm^3 s^-1] and pkernel. Indices correspond to aritrary array of columns <ic, iy>
7 : !! vertical level <k>, aerosol groups <j1,j2> and bins <i1,i2> of colliding particles.
8 : !!
9 : !! ckernel is calculated as a static array for use each timestep
10 : !! ckern0 is also created for a basis to calculate new ckernels each timestep, if desired. (coagwet.f)
11 : !!
12 : !! This routine requires that vertical profiles of temperature <T>,
13 : !! air density <rhoa>, and viscosity <rmu> are defined.
14 : !!
15 : !! @version Oct-1995
16 : !! @author Andy Ackerman
17 1050624 : subroutine setupckern(carma, cstate, rc)
18 :
19 : ! types
20 : use carma_precision_mod
21 : use carma_enums_mod
22 : use carma_constants_mod
23 : use carma_types_mod
24 : use carmastate_mod
25 : use carma_mod
26 :
27 : implicit none
28 :
29 : type(carma_type), intent(in) :: carma !! the carma object
30 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
31 : integer, intent(inout) :: rc !! return code, negative indicates failure
32 :
33 : ! Local declarations
34 : ! 2-D collision efficiency for current group pair under
35 : ! consideration (for extrapolation of input data)
36 2101248 : real(kind=f) :: e_coll2(NBIN,NBIN)
37 : integer, parameter :: NP_DATA = 21 ! number of collector/collected pairs in input data
38 : integer, parameter :: NR_DATA = 12 ! number of radius bins in input data
39 : real(kind=f), parameter :: e_small = 0.0001_f ! smallest collision efficiency
40 : logical, save :: init_data = .FALSE. ! did data_p and data_r get initialized?
41 : real(kind=f), save :: data_p(NP_DATA) ! radius ratios (collected/collector)
42 : real(kind=f), save :: data_r(NR_DATA) ! collector drop radii (um)
43 : real(kind=f), save :: data_e(NP_DATA, NR_DATA) ! geometric collection efficiencies
44 :
45 : integer :: ip
46 : integer :: ig, jg
47 :
48 : ! The probability that two particles that collide through thermal
49 : ! coagulation will stick to each other.
50 : real(kind=f) :: cstick_calc
51 :
52 : integer :: i1, i2, j1, j2, k
53 : integer :: i, j
54 : integer :: igrp
55 : integer :: ibin
56 :
57 : real(kind=f) :: rhoa_cgs
58 : real(kind=f) :: temp1, temp2
59 :
60 : real(kind=f) :: r1
61 : real(kind=f) :: di
62 : real(kind=f) :: gi
63 : real(kind=f) :: rlbi
64 : real(kind=f) :: dti1
65 : real(kind=f) :: dti2
66 : real(kind=f) :: dti
67 :
68 : real(kind=f) :: r2
69 : real(kind=f) :: dj
70 : real(kind=f) :: gj
71 : real(kind=f) :: rlbj
72 : real(kind=f) :: dtj1
73 : real(kind=f) :: dtj2
74 : real(kind=f) :: dtj
75 :
76 : real(kind=f) :: rp
77 : real(kind=f) :: dp
78 : real(kind=f) :: gg
79 : real(kind=f) :: delt
80 : real(kind=f) :: term1
81 : real(kind=f) :: term2
82 : real(kind=f) :: cbr
83 :
84 : real(kind=f) :: r_larg
85 : real(kind=f) :: r_smal
86 : integer :: i_larg
87 : integer :: i_smal
88 : integer :: ig_larg
89 : integer :: ig_smal
90 : real(kind=f) :: d_larg
91 :
92 : real(kind=f) :: re_larg
93 : real(kind=f) :: pe
94 : real(kind=f) :: pe3
95 : real(kind=f) :: ccd
96 :
97 : real(kind=f) :: e_coll
98 : real(kind=f) :: vfc_smal
99 : real(kind=f) :: vfc_larg
100 : real(kind=f) :: sk
101 : real(kind=f) :: e1
102 : real(kind=f) :: e3
103 : real(kind=f) :: e_langmuir
104 : real(kind=f) :: re60
105 :
106 : real(kind=f) :: pr
107 : real(kind=f) :: e_fuchs
108 :
109 : integer :: jp, jj, jr
110 :
111 : real(kind=f) :: pblni
112 : real(kind=f) :: rblni
113 :
114 : real(kind=f) :: term3
115 : real(kind=f) :: term4
116 :
117 : real(kind=f) :: beta
118 : real(kind=f) :: b_coal
119 : real(kind=f) :: a_coal
120 : real(kind=f) :: x_coal
121 : real(kind=f) :: e_coal
122 : real(kind=f) :: vfc_1
123 : real(kind=f) :: vfc_2
124 : real(kind=f) :: cgr
125 :
126 :
127 : ! Add constants for calculating effect of Van Der Waal's forces on coagulation
128 : ! See Chan and Mozurkewich, J. Atmos. Sci., June 2001
129 : real(kind=f), parameter :: vwa1 = 0.0757_f
130 : real(kind=f), parameter :: vwa3 = 0.0015_f
131 : real(kind=f), parameter :: vwb0 = 0.0151_f
132 : real(kind=f), parameter :: vwb1 = -0.186_f
133 : real(kind=f), parameter :: vwb3 = -0.0163_f
134 : real(kind=f), parameter :: ham = 6.4e-13_f ! erg, Hamaker constant
135 : real(kind=f) :: hp, hpln, Enot, Einf
136 2101248 : logical :: use_vw(NGROUP, NGROUP)
137 : integer :: ielem
138 :
139 :
140 : ! Initialization of input data for gravitational collection.
141 : ! The data were compiled by Hall (J. Atmos. Sci. 37, 2486-2507, 1980).
142 :
143 : data data_p/0.00_f,0.05_f,0.10_f,0.15_f,0.20_f,0.25_f,0.30_f,0.35_f,0.40_f,0.45_f, &
144 : 0.50_f,0.55_f,0.60_f,0.65_f,0.70_f,0.75_f,0.80_f,0.85_f,0.90_f,0.95_f,1.00_f/
145 :
146 : data data_r( 1), (data_e(ip, 1),ip=1,NP_DATA) / 10.0_f, &
147 : 0.0001_f, 0.0001_f, 0.0001_f, 0.0001_f, 0.0140_f, 0.0170_f, 0.0190_f, 0.0220_f, &
148 : 0.0270_f, 0.0300_f, 0.0330_f, 0.0350_f, 0.0370_f, 0.0380_f, 0.0380_f, 0.0370_f, &
149 : 0.0360_f, 0.0350_f, 0.0320_f, 0.0290_f, 0.0270_f /
150 : data data_r( 2), (data_e(ip, 2),ip=1,NP_DATA) / 20.0_f, &
151 : 0.0001_f, 0.0001_f, 0.0001_f, 0.0050_f, 0.0160_f, 0.0220_f, 0.0300_f, 0.0430_f, &
152 : 0.0520_f, 0.0640_f, 0.0720_f, 0.0790_f, 0.0820_f, 0.0800_f, 0.0760_f, 0.0670_f, &
153 : 0.0570_f, 0.0480_f, 0.0400_f, 0.0330_f, 0.0270_f /
154 : data data_r( 3), (data_e(ip, 3),ip=1,NP_DATA) / 30.0_f, &
155 : 0.0001_f, 0.0001_f, 0.0020_f, 0.0200_f, 0.0400_f, 0.0850_f, 0.1700_f, 0.2700_f, &
156 : 0.4000_f, 0.5000_f, 0.5500_f, 0.5800_f, 0.5900_f, 0.5800_f, 0.5400_f, 0.5100_f, &
157 : 0.4900_f, 0.4700_f, 0.4500_f, 0.4700_f, 0.5200_f /
158 : data data_r( 4), (data_e(ip, 4),ip=1,NP_DATA) / 40.0_f, &
159 : 0.0001_f, 0.0010_f, 0.0700_f, 0.2800_f, 0.5000_f, 0.6200_f, 0.6800_f, 0.7400_f, &
160 : 0.7800_f, 0.8000_f, 0.8000_f, 0.8000_f, 0.7800_f, 0.7700_f, 0.7600_f, 0.7700_f, &
161 : 0.7700_f, 0.7800_f, 0.7900_f, 0.9500_f, 1.4000_f /
162 : data data_r( 5), (data_e(ip, 5),ip=1,NP_DATA) / 50.0_f, &
163 : 0.0001_f, 0.0050_f, 0.4000_f, 0.6000_f, 0.7000_f, 0.7800_f, 0.8300_f, 0.8600_f, &
164 : 0.8800_f, 0.9000_f, 0.9000_f, 0.9000_f, 0.9000_f, 0.8900_f, 0.8800_f, 0.8800_f, &
165 : 0.8900_f, 0.9200_f, 1.0100_f, 1.3000_f, 2.3000_f /
166 : data data_r( 6), (data_e(ip, 6),ip=1,NP_DATA) / 60.0_f, &
167 : 0.0001_f, 0.0500_f, 0.4300_f, 0.6400_f, 0.7700_f, 0.8400_f, 0.8700_f, 0.8900_f, &
168 : 0.9000_f, 0.9100_f, 0.9100_f, 0.9100_f, 0.9100_f, 0.9100_f, 0.9200_f, 0.9300_f, &
169 : 0.9500_f, 1.0000_f, 1.0300_f, 1.7000_f, 3.0000_f /
170 : data data_r( 7), (data_e(ip, 7),ip=1,NP_DATA) / 70.0_f, &
171 : 0.0001_f, 0.2000_f, 0.5800_f, 0.7500_f, 0.8400_f, 0.8800_f, 0.9000_f, 0.9200_f, &
172 : 0.9400_f, 0.9500_f, 0.9500_f, 0.9500_f, 0.9500_f, 0.9500_f, 0.9500_f, 0.9700_f, &
173 : 1.0000_f, 1.0200_f, 1.0400_f, 2.3000_f, 4.0000_f /
174 : data data_r( 8), (data_e(ip, 8),ip=1,NP_DATA) / 100.0_f, &
175 : 0.0001_f, 0.5000_f, 0.7900_f, 0.9100_f, 0.9500_f, 0.9500_f, 1.0000_f, 1.0000_f, &
176 : 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, &
177 : 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f /
178 : data data_r( 9), (data_e(ip, 9),ip=1,NP_DATA) / 150.0_f, &
179 : 0.0001_f, 0.7700_f, 0.9300_f, 0.9700_f, 0.9700_f, 1.0000_f, 1.0000_f, 1.0000_f, &
180 : 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, &
181 : 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f /
182 : data data_r(10), (data_e(ip,10),ip=1,NP_DATA) / 200.0_f, &
183 : 0.0001_f, 0.8700_f, 0.9600_f, 0.9800_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, &
184 : 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, &
185 : 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f /
186 : data data_r(11), (data_e(ip,11),ip=1,NP_DATA) / 300.0_f, &
187 : 0.0001_f, 0.9700_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, &
188 : 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, &
189 : 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f /
190 : data data_r(12), (data_e(ip,12),ip=1,NP_DATA) / 1000.0_f, &
191 : 0.0001_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, &
192 : 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, &
193 : 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f, 1.0000_f /
194 :
195 :
196 : ! Use constant kernel if <icoagop> = I_COAGOP_CONST
197 1050624 : if( icoagop .eq. I_COAGOP_CONST )then
198 0 : ckernel(:,:,:,:,:) = ck0
199 : else
200 :
201 1050624 : if( icollec .eq. I_COLLEC_DATA )then
202 :
203 : ! Convert <data_r> from um to cm and take logarithm of <data_e>;
204 : ! however, we only want to do this once.
205 : !
206 : ! If we are using Open/MP, we only want one thread to do this
207 : ! operation once. This is a kludge, and this table should probably
208 : ! get set up a different way.
209 : !$OMP CRITICAL(CARMA_HALL)
210 1050624 : if (.not. init_data) then
211 1536 : init_data = .TRUE.
212 :
213 19968 : do i = 1, NR_DATA
214 18432 : data_r(i) = data_r(i)/1.e4_f
215 407040 : do ip = 1, NP_DATA
216 405504 : data_e(ip,i) = log(data_e(ip,i))
217 : enddo
218 : enddo
219 : endif
220 : !$OMP END CRITICAL(CARMA_HALL)
221 : endif
222 :
223 : ! Loop over the grid
224 74594304 : do k = 1, NZ
225 :
226 : ! This is <rhoa> in cartesian coordinates.
227 73543680 : rhoa_cgs = rhoa(k) / zmet(k)
228 :
229 73543680 : temp1 = BK*t(k)
230 73543680 : temp2 = 6._f*PI*rmu(k)
231 :
232 220631040 : do j1 = 1, NGROUP
233 441262080 : do j2 = j1, NGROUP
234 367718400 : use_vw(j1, j2) = is_grp_sulfate(j1) .and. is_grp_sulfate(j2)
235 : end do
236 : end do
237 :
238 : ! Loop over groups!
239 221681664 : do j1 = 1, NGROUP
240 514805760 : do j2 = 1, NGROUP
241 :
242 441262080 : if( icoag(j1,j2) .ne. 0 )then
243 :
244 : ! First particle
245 6177669120 : do i1 = 1, NBIN
246 :
247 5883494400 : r1 = r_wet(k,i1,j1) * rrat(i1,j1)
248 5883494400 : di = temp1*bpm(k,i1,j1)/(temp2*r1)
249 5883494400 : gi = sqrt( 8._f*temp1/(PI*rmass(i1,j1)) )
250 5883494400 : rlbi = 8._f*di/(PI*gi)
251 5883494400 : dti1= (2._f*r1 + rlbi)**3
252 5883494400 : dti2= (4._f*r1*r1 + rlbi*rlbi)**1.5_f
253 5883494400 : dti = 1._f/(6._f*r1*rlbi)
254 5883494400 : dti = dti*(dti1 - dti2) - 2._f*r1
255 :
256 : ! Second particle
257 >12384*10^7 : do i2 = 1, NBIN
258 >11766*10^7 : r2 = r_wet(k,i2,j2) * rrat(i2,j2)
259 >11766*10^7 : dj = temp1*bpm(k,i2,j2)/(temp2*r2)
260 >11766*10^7 : gj = sqrt( 8._f*temp1/(PI*rmass(i2,j2)) )
261 >11766*10^7 : rlbj = 8._f*dj/(PI*gj)
262 >11766*10^7 : dtj1= (2._f*r2 + rlbj)**3
263 >11766*10^7 : dtj2= (4._f*r2*r2 + rlbj*rlbj)**1.5_f
264 >11766*10^7 : dtj = 1._f/(6._f*r2*rlbj)
265 >11766*10^7 : dtj = dtj*(dtj1 - dtj2) - 2._f*r2
266 :
267 : ! Account for the charging effect of small particles (Van Der Waal's forces).
268 : ! Set cstick to E_infinity/Eo, then multiply cbr kernel by Eo
269 : ! See Chan and Mozurkewich, J. Atmos. Sci., June 2001
270 : ! Only applicable to groups with sulfate elements
271 >11766*10^7 : if (use_vw(j1,j2)) then
272 29418312000 : hp = ham / temp1 * (4._f * r1 * r2 / (r1 + r2)**2)
273 29418312000 : hpln = log(1._f + hp)
274 29418312000 : Enot = 1._f + vwa1 * hpln + vwa3 * hpln**3
275 29418312000 : Einf = 1._f + sqrt(hp / 3._f) / (1._f + vwb0*sqrt(hp)) + vwb1 * hpln + vwb3 * hpln**3
276 29418312000 : cstick_calc = Einf / Enot
277 : else
278 88251576000 : cstick_calc = cstick
279 : end if
280 :
281 : ! First calculate thermal coagulation kernel
282 >11766*10^7 : rp = r1 + r2
283 >11766*10^7 : dp = di + dj
284 >11766*10^7 : gg = sqrt(gi*gi + gj*gj)*cstick_calc
285 >11766*10^7 : delt= sqrt(dti*dti + dtj*dtj)
286 >11766*10^7 : term1 = rp/(rp + delt)
287 >11766*10^7 : term2 = 4._f*dp/(gg*rp)
288 :
289 : ! <cbr> is thermal (brownian) coagulation coefficient
290 >11766*10^7 : cbr = 4._f*PI*rp*dp/(term1 + term2)
291 :
292 : ! Determine indices of larger and smaller particles (of the pair)
293 >11766*10^7 : if (r2 .ge. r1) then
294 : r_larg = r2
295 : r_smal = r1
296 : i_larg = i2
297 : i_smal = i1
298 : ig_larg = j2
299 : ig_smal = j1
300 : d_larg = dj
301 : else
302 57364070400 : r_larg = r1
303 57364070400 : r_smal = r2
304 57364070400 : i_larg = i1
305 57364070400 : i_smal = i2
306 57364070400 : ig_larg = j1
307 57364070400 : ig_smal = j2
308 57364070400 : d_larg = di
309 : endif
310 :
311 : ! Calculate enhancement of coagulation due to convective diffusion
312 : ! as described in Pruppacher and Klett (Eqs. 17-12 and 17-14).
313 :
314 : ! Enhancement applies to larger particle.
315 >11766*10^7 : re_larg = re(k,i_larg,ig_larg)
316 :
317 : ! <pe> is Peclet number.
318 >11766*10^7 : pe = re_larg*rmu(k) / (rhoa_cgs*d_larg)
319 >11766*10^7 : pe3 = pe**(1._f/3._f)
320 :
321 : ! <ccd> is convective diffusion coagulation coefficient
322 >11766*10^7 : if (use_ccd(j1,j2)) then
323 : ! Convective diffusion coagulation can be large in specific
324 : ! scavenging processes such as turbulence environment inside
325 : ! volcanic plume or raindrop washing away the aerosols.
326 : ! use_ccd should only set to be .true. if doing processes mentioned above.
327 0 : if( re_larg .lt. 1._f )then
328 0 : ccd = 0.45_f*cbr*pe3
329 : else
330 0 : ccd = 0.45_f*cbr*pe3*re_larg**(ONE/6._f)
331 : endif
332 : else
333 : ! all other conditions, use_ccd should set to .false.
334 : ! and use_ccd should be .false. as default
335 : ccd = 0._f
336 : end if
337 :
338 : ! Next calculate gravitational collection kernel.
339 :
340 : ! First evaluate collection efficiency <e>.
341 >11766*10^7 : if( icollec .eq. I_COLLEC_CONST )then
342 : ! constant value
343 0 : e_coll = grav_e_coll0
344 >11766*10^7 : else if( icollec .eq. I_COLLEC_FUCHS )then
345 : ! Find maximum of Langmuir's formulation and Fuchs' value.
346 : ! First calculate Langmuir's efficiency <e_langmuir>.
347 :
348 : ! <sk> is stokes number.
349 : ! <vfc_{larg,smal}> is the fallspeed in cartesian coordinates.!
350 0 : vfc_smal = vf(k,i_smal,ig_smal) * zmet(k)
351 0 : vfc_larg = vf(k,i_larg,ig_larg) * zmet(k)
352 :
353 0 : sk = vfc_smal * (vfc_larg - vfc_smal) / (r_larg*GRAV)
354 :
355 0 : if( sk .lt. 0.08333334_f )then
356 : e1 = 0._f
357 : else
358 0 : e1 = (sk/(sk + 0.25_f))**2
359 : endif
360 :
361 0 : if( sk .lt. 1.214_f )then
362 : e3 = 0._f
363 : else
364 0 : e3 = 1._f/(1._f+.75_f*log(2._f*sk)/(sk-1.214_f))**2
365 : endif
366 :
367 0 : if( re_larg .lt. 1._f )then
368 : e_langmuir = e3
369 0 : else if( re_larg .gt. 1000._f )then
370 : e_langmuir = e1
371 0 : else if( re_larg .le. 1000._f )then
372 0 : re60 = re_larg/60._f
373 0 : e_langmuir = (e3 + re60*e1)/(1._f + re60)
374 : endif
375 :
376 : ! Next calculate Fuchs' efficiency (valid for r < 10 um).
377 0 : pr = r_smal/r_larg
378 0 : e_fuchs = (pr/(1.414_f*(1._f + pr)))**2
379 :
380 0 : e_coll = max( e_fuchs, e_langmuir )
381 :
382 >11766*10^7 : else if( icollec .eq. I_COLLEC_DATA )then
383 :
384 : ! Interpolate input data (from data statment at beginning of subroutine).
385 >11766*10^7 : pr = r_smal/r_larg
386 :
387 : ! First treat cases outside the data range
388 >11766*10^7 : if( pr .lt. data_p(2) )then
389 :
390 : ! Radius ratio is smaller than lowest nonzero ratio in input data --
391 : ! use constant values (as in Beard and Ochs, 1984) if available,
392 : ! otherwise use very small efficiencty
393 51624551208 : if( i2 .eq. i_larg )then
394 26737525694 : if( i2.eq.1 )then
395 : e_coll = e_small
396 : else
397 26283439309 : e_coll = e_coll2(i1,i2-1)
398 : endif
399 : else
400 24887025514 : if( i1.eq.1 )then
401 : e_coll = e_small
402 : else
403 24580026489 : e_coll = e_coll2(i1-1,i2)
404 : endif
405 : endif
406 :
407 66045336792 : elseif( r_larg .lt. data_r(1) )then
408 : ! Radius of larger particle is smaller than smallest radius in input data --
409 : ! assign very small efficiency.
410 : e_coll = e_small
411 : else
412 :
413 : ! Both droplets are either within grid (interpolate) or larger
414 : ! droplet is larger than maximum on grid (extrapolate) -- in both cases
415 : ! will interpolate on ratio of droplet radii.
416 :
417 : ! Find <jp> such that data_p(jp) <= pr <= data_p(jp+1) and calculate
418 : ! <pblni> = fractional distance of <pr> between points in <data_p>
419 1021937106 : jp = NP_DATA
420 20438742120 : do jj = NP_DATA-1, 2, -1
421 20438742120 : if( pr .le. data_p(jj+1) ) jp = jj
422 : enddo
423 :
424 : ! should not need this if-stmt
425 1021937106 : if( jp .lt. NP_DATA )then
426 1021937106 : pblni = (pr - data_p(jp)) / (data_p(jp+1) - data_p(jp))
427 : else
428 : ! nor this else-stmt
429 0 : if (do_print) write(LUNOPRT, *) 'setupckern::ERROR NP_DATA < jp = ', jp
430 0 : return
431 : endif
432 :
433 1021937106 : if( r_larg .gt. data_r(NR_DATA) )then
434 :
435 : ! Extrapolate on R and interpolate on p
436 : !
437 : ! NOTE: This expression has a bugin it, since jr won't
438 : ! be defined.
439 0 : e_coll = (1._f-pblni)*data_e(jp ,jr) + &
440 0 : ( pblni)*data_e(jp+1,jr)
441 :
442 : else
443 :
444 : ! Find <jr> such that data_r(jr) <= r_larg <= data_r(jr+1) and calculate
445 : ! <rblni> = fractional distance of <r_larg> between points in <data_r>
446 : jr = NR_DATA
447 12263245272 : do jj = NR_DATA-1, 1, -1
448 12263245272 : if( r_larg .le. data_r(jj+1) ) jr = jj
449 : enddo
450 1021937106 : rblni = (r_larg - data_r(jr)) / (data_r(jr+1) - data_r(jr))
451 :
452 : ! Bilinear interpolation of logarithm of data.
453 : e_coll = (1._f-pblni)*(1._f-rblni)*data_e(jp ,jr ) + &
454 : ( pblni)*(1._f-rblni)*data_e(jp+1,jr ) + &
455 : (1._f-pblni)*( rblni)*data_e(jp ,jr+1) + &
456 1021937106 : ( pblni)*( rblni)*data_e(jp+1,jr+1)
457 :
458 : ! (since data_e is logarithm of efficiencies)
459 1021937106 : term1 = (1._f-rblni)*(1._f-pblni)*data_e(jp,jr)
460 :
461 : if( jp .lt. NP_DATA )then
462 1021937106 : term2 = pblni*(1._f-rblni)*data_e(jp+1,jr)
463 : else
464 : term2 = -100._f
465 : endif
466 :
467 1021937106 : if( jr .lt. NR_DATA )then
468 1021937106 : term3 = (1._f-pblni)*rblni*data_e(jp,jr+1)
469 : else
470 : term3 = -100._f
471 : endif
472 :
473 1021937106 : if( jr .lt. NR_DATA .and. jp .lt. NP_DATA )then
474 1021937106 : term4 = pblni*rblni*data_e(jp+1,jr+1)
475 : else
476 : term4 = -100._f
477 : endif
478 :
479 1021937106 : e_coll = exp(term1 + term2 + term3 + term4)
480 : endif
481 : endif
482 :
483 >11766*10^7 : e_coll2(i1,i2) = e_coll
484 : endif
485 :
486 : ! Now calculate coalescence efficiency from Beard and Ochs
487 : ! (J. Geophys. Res. 89, 7165-7169, 1984).
488 >11766*10^7 : beta = log(r_smal*1.e4_f) + 0.44_f*log(r_larg*50._f)
489 >11766*10^7 : b_coal = 0.0946_f*beta - 0.319_f
490 >11766*10^7 : a_coal = sqrt(b_coal**2 + 0.00441_f)
491 : x_coal = (a_coal-b_coal)**(ONE/3._f) &
492 >11766*10^7 : - (a_coal+b_coal)**(ONE/3._f)
493 >11766*10^7 : x_coal = x_coal + 0.459_f
494 :
495 : ! Limit extrapolated values to no less than 50% and no more than 100%
496 >11766*10^7 : x_coal = max(x_coal,.5_f)
497 >11766*10^7 : e_coal = min(x_coal,1._f)
498 :
499 : ! Now use coalescence efficiency and collision efficiency in definition
500 : ! of (geometric) gravitational collection efficiency <cgr>.
501 >11766*10^7 : vfc_1 = vf(k,i1,j1) * zmet(k)
502 >11766*10^7 : vfc_2 = vf(k,i2,j2) * zmet(k)
503 >11766*10^7 : cgr = e_coal * e_coll * PI * rp**2 * abs( vfc_1 - vfc_2 )
504 :
505 : ! Long's (1974) kernel that only depends on size of larger droplet
506 : ! if( r_larg .le. 50.e-4_f )then
507 : ! cgr = 1.1e10_f * vol(i_larg,ig_larg)**2
508 : ! else
509 : ! cgr = 6.33e3_f * vol(i_larg,ig_larg)
510 : ! endif
511 :
512 : ! Now combine all the coagulation and collection kernels into the
513 : ! overall kernel.
514 >12355*10^7 : ckernel(k,i1,i2,j1,j2) = cbr + ccd + cgr
515 :
516 : ! To avoid generation of large, non-physical hydrometeors by
517 : ! coagulation, cut down ckernel for large radii
518 : ! if( ( r1 .gt. 0.18_f .and. r2 .gt. 10.e-4_f ) .or. &
519 : ! ( r2 .gt. 0.18_f .and. r1 .gt. 10.e-4_f ) ) then
520 : ! ckernel(k,i1,i2,j1,j2) = ckernel(k,i1,i2,j1,j2) / 1.e6_f
521 : ! endif
522 :
523 : enddo ! second particle bin
524 : enddo ! first particle bin
525 : endif ! icoag ne 0
526 : enddo ! second particle group
527 : enddo ! first particle group
528 : enddo ! vertical level
529 : endif ! not constant
530 :
531 : ! return to caller with coagulation kernels evaluated.
532 : return
533 1050624 : end
|
