source: palm/trunk/SOURCE/lpm_collision_kernels.f90 @ 825

Last change on this file since 825 was 825, checked in by raasch, 10 years ago

New:
---

droplet growth by condensation may include curvature and solution effects.
Steered by new inipar-parameter curvature_solution_effects.
(advec_particles, check_parameters, header, init_cloud_physics, init_particles, modules, parin, read_var_list, write_var_list)

mean/min/max particle radius added as output quantity. (data_output_ptseries, modules)

Changed:


Initialisation of temporary particle array for resorting removed.
(advec_particles)

particle attributes speed_x|y|z_sgs renamed rvar1|2|3.
(advec_particles, data_output_ptseries, modules, init_particles, particle_boundary_conds)

routine wang_kernel and respective module renamed lpm_collision_kernels.
Package (particle) parameters wang_collision_kernel and turbulence_effects_on_collision
replaced by parameter collision_kernel.
(Makefile, advec_particles, check_parameters, diffusion_e, init_3d_model, modules, package_parin, time_integration, new: lpm_collision_kernels, deleted: wang_kernel)

Errors:


  • Property svn:keywords set to Id
File size: 24.9 KB
Line 
1MODULE lpm_collision_kernels_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! routine renamed from wang_kernel to lpm_collision_kernels,
7! turbulence_effects on collision replaced by wang_kernel
8!
9! Former revisions:
10! -----------------
11! $Id: lpm_collision_kernels.f90 825 2012-02-19 03:03:44Z raasch $
12!
13! 799 2011-12-21 17:48:03Z franke
14! speed optimizations and formatting
15! Bugfix: iq=1 is not allowed (routine effic)
16! Bugfix: replaced stop by ec=0.0 in case of very small ec (routine effic)
17!
18! 790 2011-11-29 03:11:20Z raasch
19! initial revision
20!
21! Description:
22! ------------
23! This routine calculates the effect of (SGS) turbulence on the collision
24! efficiency of droplets.
25! It is based on the original kernel developed by Wang (...)
26!------------------------------------------------------------------------------!
27
28    USE arrays_3d
29    USE cloud_parameters
30    USE constants
31    USE particle_attributes
32
33    IMPLICIT NONE
34
35    PRIVATE
36
37    PUBLIC  colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi
38
39    INTEGER, SAVE ::  ip, jp, kp, pend, pstart, psum
40
41    REAL, SAVE :: epsilon, eps2, urms, urms2
42
43    REAL, DIMENSION(:),   ALLOCATABLE, SAVE ::  winf
44    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  ec, ecf, gck
45
46!
47!-- Public interfaces
48    INTERFACE colker
49       MODULE PROCEDURE colker
50    END INTERFACE colker
51
52    INTERFACE effic
53       MODULE PROCEDURE effic
54    END INTERFACE effic
55
56    INTERFACE fallg
57       MODULE PROCEDURE fallg
58    END INTERFACE fallg
59
60    INTERFACE phi
61       MODULE PROCEDURE phi
62    END INTERFACE phi
63
64    INTERFACE turbsd
65       MODULE PROCEDURE turbsd
66    END INTERFACE turbsd
67
68    INTERFACE turb_enhance_eff
69       MODULE PROCEDURE turb_enhance_eff
70    END INTERFACE turb_enhance_eff
71
72    INTERFACE zhi
73       MODULE PROCEDURE zhi
74    END INTERFACE zhi
75
76    CONTAINS
77
78!------------------------------------------------------------------------------!
79! SUBROUTINE for calculation of collision kernel
80!------------------------------------------------------------------------------!
81    SUBROUTINE colker( i1, j1, k1, kernel )
82
83       USE arrays_3d
84       USE cloud_parameters
85       USE constants
86       USE cpulog
87       USE indices
88       USE interfaces
89       USE particle_attributes
90
91       IMPLICIT NONE
92
93       INTEGER ::  i, i1, j, j1, k1
94
95       REAL, DIMENSION(prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1,prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1) :: kernel
96
97!       CALL cpu_log( log_point_s(46), 'colker', 'start' )
98
99       ip = i1
100       jp = j1
101       kp = k1
102
103       pstart = prt_start_index(kp,jp,ip)
104       pend   = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1
105       psum   = prt_count(kp,jp,ip)
106
107       ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
108
109       IF ( wang_kernel )  THEN
110
111          ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
112
113          epsilon = diss(kp,jp,ip) * 1.0E5  !dissipation rate in cm**2/s**-3
114          urms     = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 )
115
116          IF ( epsilon <= 0.001 )  THEN
117
118             CALL fallg
119             CALL effic
120
121             DO  j = pstart, pend
122                DO  i =  pstart, pend
123                   kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
124                                        particles(i)%radius * 100.0 )**2 &
125                                    * ec(i,j) * ABS( winf(j) - winf(i) )
126                ENDDO
127             ENDDO
128
129          ELSE
130
131             CALL turbsd
132             CALL turb_enhance_eff
133             CALL effic
134
135             DO  j = pstart, pend, 1
136                DO  i =  pstart, pend, 1
137                   kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
138                ENDDO
139             ENDDO
140
141          ENDIF
142
143          DEALLOCATE(gck, ecf)
144
145       ELSE
146
147!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )
148          CALL fallg
149!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )
150!          CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )
151          CALL effic
152!          CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )
153
154          DO  j = pstart, pend
155             DO  i =  pstart, pend
156                kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
157                                     particles(i)%radius * 100.0 )**2 &
158                                 * ec(i,j) * ABS( winf(j) - winf(i) )
159             ENDDO
160          ENDDO
161
162       ENDIF
163
164       DEALLOCATE( ec, winf )
165
166!       CALL cpu_log( log_point_s(46), 'colker', 'stop' )
167
168    END SUBROUTINE colker
169
170!------------------------------------------------------------------------------!
171! SUBROUTINE for calculation of w, g and gck
172!------------------------------------------------------------------------------!
173    SUBROUTINE turbsd
174! from Aayala 2008b, page 37ff, necessary input parameter water density, radii
175! of droplets, air density, air viscosity, turbulent dissipation rate,
176! taylor microscale reynolds number, gravitational acceleration
177
178       USE constants
179       USE cloud_parameters
180       USE particle_attributes
181       USE arrays_3d
182
183       IMPLICIT NONE
184
185       REAL :: Relamda, &
186               Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
187               bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
188               vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, &
189               wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
190
191       REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens
192
193       REAL, DIMENSION(pstart:pend) :: St, tau
194
195       LOGICAL, SAVE ::  first = .TRUE.
196
197       INTEGER :: i, j
198
199!
200!--   Initial assignment of constants
201       IF ( first )  THEN
202
203          first = .FALSE.
204          airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
205          airdens   = 1.2250   !air density in mg/cm**3
206          waterdens = 1000.0   !water density in mg/cm**3
207          gravity   = 980.6650 !in cm/s**2
208          anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
209
210       ENDIF   
211
212       Relamda = urms**2*sqrt(15.0/epsilon/anu)
213
214       Tl = urms**2/epsilon       !in s
215       Lf = 0.5 * (urms**3)/epsilon !in cm
216
217       tauk = (anu/epsilon)**0.5       !in s
218       eta  = (anu**3/epsilon)**0.25  !in cm
219       vk   = eta/tauk                   !in cm/s
220
221       ao = (11.+7.*Relamda)/(205.+Relamda)
222
223       lambda = urms * sqrt(15.*anu/epsilon)     !in cm
224
225       tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s
226
227       CALL fallg !gives winf in cm/s
228
229       DO i = pstart, pend
230          tau(i) = winf(i)/gravity    !in s
231          St(i)  = tau(i)/tauk
232       ENDDO
233
234!***** TO CALCULATE wr ********************************
235!from Aayala 2008b, page 38f
236
237       z  = tt/Tl
238       be = sqrt(2.0)*lambda/Lf
239
240       bbb = sqrt(1.0-2.0*be**2)
241       d1 = (1.+bbb)/2.0/bbb
242       e1 = Lf*(1.0+bbb)/2.0   !in cm
243       d2 = (1.0-bbb)/2.0/bbb
244       e2 = Lf*(1.0-bbb)/2.0   !in cm
245
246       ccc = sqrt(1.0-2.0*z**2)
247       b1 = (1.+ccc)/2./ccc
248       c1 = Tl*(1.+ccc)/2.   !in s
249       b2 = (1.-ccc)/2./ccc
250       c2 = Tl*(1.-ccc)/2.   !in s
251
252       DO i = pstart, pend
253          v1 = winf(i)         !in cm/s
254          t1 = tau(i)         !in s
255
256          DO j = pstart,i
257             rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm
258             v2  = winf(j)         !in cm/s
259             t2  = tau(j)         !in s
260
261             v1xysq =  b1*d1*PHI(c1,e1,v1,t1) &
262                     - b1*d2*PHI(c1,e2,v1,t1) &
263                     - b2*d1*PHI(c2,e1,v1,t1) &
264                     + b2*d2*PHI(c2,e2,v1,t1)
265             v1xysq = v1xysq * urms**2/t1    !in cm**2/s**2
266
267             vrms1xy= sqrt(v1xysq)             !in cm/s
268
269             v2xysq =  b1*d1*PHI(c1,e1,v2,t2) &
270                     - b1*d2*PHI(c1,e2,v2,t2) &
271                     - b2*d1*PHI(c2,e1,v2,t2) &
272                     + b2*d2*PHI(c2,e2,v2,t2)
273             v2xysq = v2xysq * urms**2/t2    !in cm**2/s**2
274
275              vrms2xy= sqrt(v2xysq)             !in cm/s
276
277             IF(winf(i).ge.winf(j)) THEN
278                v1 = winf(i)
279                t1 = tau(i)
280                v2 = winf(j)
281                t2 = tau(j)
282             ELSE
283                v1 = winf(j)
284                t1 = tau(j)
285                v2 = winf(i)
286                t2 = tau(i)
287             ENDIF
288
289             v1v2xy =  b1*d1*ZHI(c1,e1,v1,t1,v2,t2) &
290                     - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) &
291                     - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) &
292                     + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
293             fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
294             v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j)    !in cm**2/s**2
295
296             wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy  !in cm**2/s**2
297
298              IF (wrtur2xy.le.0.0) wrtur2xy=0.0
299
300              wrgrav2=pi/8.*(winf(j)-winf(i))**2
301
302              wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
303
304
305!***** TO CALCULATE gr ********************************
306
307             IF(St(j).gt.St(i)) THEN
308                SSt = St(j)
309             ELSE
310                SSt = St(i)
311             ENDIF
312
313             XX = -0.1988*SSt**4 + 1.5275*SSt**3 &
314                  -4.2942*SSt**2 + 5.3406*SSt
315
316             IF(XX.le.0.0) XX = 0.0
317
318             YY = 0.1886*exp(20.306/Relamda)
319
320             c1_gr = XX/(gravity/(vk/tauk))**YY
321
322             ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2
323             fao_gr = 20.115 * (ao_gr/Relamda)**0.5
324             rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
325
326             grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.)
327             IF (grFIN.lt.1.0) grFIN = 1.0
328
329             gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN       ! in cm**3/s
330             gck(j,i) = gck(i,j)
331
332          ENDDO
333       ENDDO
334
335    END SUBROUTINE TurbSD
336
337!------------------------------------------------------------------------------!
338! PHI as a function
339!------------------------------------------------------------------------------!
340  REAL FUNCTION PHI(a,b,vsett,tau0)
341
342       IMPLICIT NONE
343
344       REAL :: a, aa1, b, vsett, tau0
345
346       aa1 = 1./tau0 + 1./a + vsett/b
347
348       PHI = 1./aa1 - vsett/2.0/b/aa1**2  !in s
349
350    END FUNCTION PHI
351
352!------------------------------------------------------------------------------!
353! ZETA as a function
354!------------------------------------------------------------------------------!
355  REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)
356
357       IMPLICIT NONE
358
359       REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2
360
361        aa1 = vsett2/b - 1./tau2 - 1./a
362        aa2 = vsett1/b + 1./tau1 + 1./a
363        aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
364        aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2
365        aa5 = vsett2/b + 1./tau2 + 1./a
366        aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2
367        ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2          &
368             + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6      &
369             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2)   &
370                                                           * 1./2./b/aa3             ! in s**2
371
372    END FUNCTION ZHI
373
374!------------------------------------------------------------------------------!
375! SUBROUTINE for calculation of terminal velocity winf
376!------------------------------------------------------------------------------!
377   SUBROUTINE fallg
378
379       USE constants
380       USE cloud_parameters
381       USE particle_attributes
382       USE arrays_3d
383
384      IMPLICIT NONE
385
386      INTEGER :: i, j
387
388      LOGICAL, SAVE ::  first = .TRUE.
389
390      REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py
391
392      REAL :: bond, x, xrey, y
393
394      REAL, DIMENSION(1:7), SAVE  :: b
395      REAL, DIMENSION(1:6), SAVE  :: c
396
397!
398!--   Initial assignment of constants
399      IF ( first )  THEN
400
401         first = .FALSE.
402         b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
403             -0.578878e-3,0.855176e-4,-0.327815e-5/)
404         c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
405             -0.542819e-1, 0.238449e-2/)
406
407         eta   = 1.818e-4     !in poise = g/(cm s)
408         xlamb = 6.62e-6    !in cm
409
410         rhow = 1.0       !in g/cm**3
411         rhoa = 1.225e-3    !in g/cm**3
412
413         grav = 980.665     !in cm/s**2
414         cunh = 1.257 * xlamb !in cm
415         t0   = 273.15        !in K
416         sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
417         stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
418         stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
419         phy  = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa))
420         py   = phy**(1.0/6.0)
421
422      ENDIF
423
424!particle radius has to be in cm
425      DO j = pstart, pend
426
427         IF (particles(j)%radius*100.0 .le. 1.e-3) THEN
428
429            winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s
430
431         ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN
432
433            x = log(stb*(particles(j)%radius*100.0)**3)
434            y = 0.0
435
436            DO i = 1, 7
437               y = y + b(i) * (x**(i-1))
438            ENDDO
439
440            xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y)
441            winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
442
443         ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
444
445            IF (particles(j)%radius*100.0.gt.0.35)  THEN
446               bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma
447            ELSE
448               bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma
449            ENDIF
450
451            x = log(16.0*bond*py/3.0)
452            y = 0.0
453
454            DO i = 1, 6
455               y = y + c(i) * (x**(i-1))
456            ENDDO
457
458            xrey = py*exp(y)
459
460            IF (particles(j)%radius*100.0 .gt.0.35)  THEN
461              winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
462            ELSE
463               winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
464            ENDIF
465
466         ENDIF
467      ENDDO
468      RETURN
469   END SUBROUTINE fallg
470
471!------------------------------------------------------------------------------!
472! SUBROUTINE for calculation of collision efficencies
473!------------------------------------------------------------------------------!
474
475   SUBROUTINE effic
476
477      USE arrays_3d
478      USE constants
479      USE cloud_parameters
480      USE particle_attributes
481
482!collision efficiencies of hall kernel
483      IMPLICIT NONE
484
485      INTEGER :: i, ir, iq, j, k, kk
486
487      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
488
489      LOGICAL, SAVE ::  first = .TRUE.
490
491      REAL ::  ek, particle_radius, pp, qq, rq
492
493      REAL, DIMENSION(1:21), SAVE :: rat
494      REAL, DIMENSION(1:15), SAVE :: r0
495      REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
496
497!
498!--   Initial assignment of constants
499      IF ( first )  THEN
500
501         first = .FALSE.
502         r0  = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., &
503                 300./)
504         rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,  &
505                 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
506
507         ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, &
508                        0.001,0.001,0.001,0.001,0.001,0.001/)
509         ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, &
510                        0.050,0.200,0.500,0.770,0.870,0.970/)
511         ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, &
512                        0.430,0.580,0.790,0.930,0.960,1.000/)
513         ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, &
514                        0.640,0.750,0.910,0.970,0.980,1.000/)
515         ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, &
516                        0.770,0.840,0.950,0.970,1.000,1.000/)
517         ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, &
518                        0.840,0.880,0.950,1.000,1.000,1.000/)
519         ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, &
520                        0.870,0.900,0.950,1.000,1.000,1.000/)
521         ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, &
522                        0.890,0.920,1.000,1.000,1.000,1.000/)
523         ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, &
524                        0.900,0.940,1.000,1.000,1.000,1.000/)
525         ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, &
526                        0.910,0.950,1.000,1.000,1.000,1.000/)
527         ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, &
528                        0.910,0.950,1.000,1.000,1.000,1.000/)
529         ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, &
530                        0.910,0.950,1.000,1.000,1.000,1.000/)
531         ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, &
532                        0.910,0.950,1.000,1.000,1.000,1.000/)
533         ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, &
534                        0.910,0.950,1.000,1.000,1.000,1.000/)
535         ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, &
536                        0.920,0.950,1.000,1.000,1.000,1.000/)
537         ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, &
538                        0.930,0.970,1.000,1.000,1.000,1.000/)
539         ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, &
540                        0.950,1.000,1.000,1.000,1.000,1.000/)
541         ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, &
542                        1.000,1.020,1.020,1.020,1.020,1.020/)
543         ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, &
544                        1.030,1.040,1.040,1.040,1.040,1.040/)
545         ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, &
546                        1.700,2.300,2.300,2.300,2.300,2.300/)
547         ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, &
548                        3.000,4.000,4.000,4.000,4.000,4.000/)
549      ENDIF
550
551!
552!--   Calculate the radius class index of particles with respect to array r
553      ALLOCATE( ira(pstart:pend) )
554      DO  j = pstart, pend
555         particle_radius = particles(j)%radius * 1.0E6
556         DO  k = 1, 15
557            IF ( particle_radius < r0(k) )  THEN
558               ira(j) = k
559               EXIT
560            ENDIF
561         ENDDO
562         IF ( particle_radius >= r0(15) )  ira(j) = 16
563      ENDDO
564
565!
566!--   Two-dimensional linear interpolation of the collision efficiency.
567!--   Radius has to be in µm
568      DO  j = pstart, pend
569         DO  i = pstart, j
570
571            ir = ira(j)
572
573            rq = particles(i)%radius / particles(j)%radius
574
575!            DO kk = 2, 21
576!               IF ( rq <= rat(kk) )  THEN
577!                  iq = kk
578!                  EXIT
579!               ENDIF
580!            ENDDO
581
582            iq = INT( rq * 20 ) + 1
583            iq = MAX(iq , 2)
584
585            IF ( ir < 16 )  THEN
586
587               IF ( ir >= 2 )  THEN
588                  pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / &
589                       ( r0(ir) - r0(ir-1) )
590                  qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
591                  ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)   &
592                            + pp * ( 1.0-qq ) * ecoll(ir,iq-1)           &
593                            + qq * ( 1.0-pp ) * ecoll(ir-1,iq)           &
594                            + pp * qq * ecoll(ir,iq)
595               ELSE
596                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
597                  ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
598               ENDIF
599
600            ELSE
601               qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
602               ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
603               ec(j,i) = MIN( ek, 1.0 )
604            ENDIF
605
606            ec(i,j) = ec(j,i)
607            IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
608
609         ENDDO
610      ENDDO
611
612      DEALLOCATE( ira )
613
614   END SUBROUTINE effic
615
616!------------------------------------------------------------------------------!
617! SUBROUTINE for calculation of enhancement factor collision efficencies
618!------------------------------------------------------------------------------!
619   SUBROUTINE turb_enhance_eff
620
621       USE constants
622       USE cloud_parameters
623       USE particle_attributes
624       USE arrays_3d
625
626      IMPLICIT NONE
627
628      INTEGER :: i, ik, ir, iq, j, k, kk
629
630      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
631
632      REAL    :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3
633
634      LOGICAL, SAVE ::  first = .TRUE.
635
636      REAL, DIMENSION(1:11), SAVE :: rat
637      REAL, DIMENSION(1:7), SAVE  :: r0
638      REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400
639
640!
641!--   Initial assignment of constants
642      IF ( first )  THEN
643
644         first = .FALSE.
645
646         r0  = (/10., 20., 30.,40., 50., 60.,100./)
647         rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
648
649! 100 cm^2/s^3
650         ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
651         ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
652         ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
653         ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
654         ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
655         ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
656         ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
657         ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
658         ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
659         ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
660         ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
661
662! 400 cm^2/s^3
663         ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
664         ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
665         ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
666         ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
667         ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
668         ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
669         ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
670         ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
671         ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
672         ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
673         ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
674
675      ENDIF
676
677!
678!--   Calculate the radius class index of particles with respect to array r
679      ALLOCATE( ira(pstart:pend) )
680
681      DO  j = pstart, pend
682         particle_radius = particles(j)%radius * 1.0E6
683         DO  k = 1, 7
684            IF ( particle_radius < r0(k) )  THEN
685               ira(j) = k
686               EXIT
687            ENDIF
688         ENDDO
689         IF ( particle_radius >= r0(7) )  ira(j) = 8
690      ENDDO
691
692! two-dimensional linear interpolation of the collision efficiency
693      DO j =  pstart, pend
694         DO i = pstart, j
695
696            ir = ira(j)
697
698            rq = particles(i)%radius/particles(j)%radius
699
700            DO kk = 2, 11
701               IF ( rq <= rat(kk) )  THEN
702                  iq = kk
703                  EXIT
704               ENDIF
705            ENDDO
706
707! 0  cm2/s3
708            y1 = 1.0
709! 100 cm2/s3, 400 cm2/s3
710            IF (ir.lt.8) THEN
711               IF (ir.ge.2) THEN
712                  pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1))
713                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
714                  y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+  &
715                            pp*(1.-qq)*ecoll_100(ir,iq-1)+    &
716                            qq*(1.-pp)*ecoll_100(ir-1,iq)+    &
717                            pp*qq*ecoll_100(ir,iq)
718                  y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+  &
719                            pp*(1.-qq)*ecoll_400(ir,iq-1)+    &
720                            qq*(1.-pp)*ecoll_400(ir-1,iq)+    &
721                            pp*qq*ecoll_400(ir,iq)
722               ELSE
723                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
724                  y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq)
725                  y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq)
726               ENDIF
727            ELSE
728            qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
729            y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
730            y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
731         ENDIF
732! linear interpolation
733! dissipation rate in cm2/s3
734         x1 = 0.0
735         x2 = 100.0
736         x3 = 400.0
737
738         IF (epsilon.le.100.) THEN
739            ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 &
740                    +  (epsilon-0.)/(100.-0.) * y2
741         ELSE IF(epsilon.le.600.)THEN
742            ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 &
743                    +  (epsilon-100.)/(400.-100.) * y3
744
745         ELSE
746            ecf(j,i) = (600.-400.)/(100.-400.) * y2 &
747                    +  (600.-100.)/(400.-100.) * y3
748         ENDIF
749
750         IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0
751
752         ecf(i,j)=ecf(j,i)
753      ENDDO
754   ENDDO
755
756   RETURN
757   END SUBROUTINE turb_enhance_eff
758
759 END MODULE lpm_collision_kernels_mod
Note: See TracBrowser for help on using the repository browser.