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