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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

  • Property svn:keywords set to Id
File size: 38.9 KB
Line 
1!> @file lpm_collision_kernels.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2016 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! Forced header and separation lines into 80 columns
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_collision_kernels.f90 2000 2016-08-20 18:09:15Z knoop $
27!
28! 1880 2016-04-20 09:36:50Z hoffmann
29! Bugfix: The index of the larger particle has to be chosen for interpolation.
30!
31! 1873 2016-04-18 14:50:06Z maronga
32! Module renamed (removed _mod)
33!
34! 1858 2016-04-13 13:12:11Z hoffmann
35! Interpolation of collision kernels adjusted to more reasonable values.
36! Reformatting of the code.
37!
38! 1850 2016-04-08 13:29:27Z maronga
39! Module renamed
40!
41! 1822 2016-04-07 07:49:42Z hoffmann
42! PALM kernel has been deleted.
43! Bugfix in the calculation of the turbulent enhancement factor of the
44! collection efficiency.
45!
46! Unused variables removed.
47!
48! 1776 2016-03-02 17:54:58Z hoffmann
49! Bugfix: Collection efficiencies must be calculated for the larger droplet.
50!
51! 1682 2015-10-07 23:56:08Z knoop
52! Code annotations made doxygen readable
53!
54! 1519 2015-01-08 10:20:42Z hoffmann
55! Bugfix: Using the new particle structure, particles are not sorted by size.
56! Hence, computation of collision efficiencies must ensure that the ratio of
57! two colliding droplets is < 1.
58!
59! 1359 2014-04-11 17:15:14Z hoffmann
60! New particle structure integrated.
61! Kind definition added to all floating point numbers.
62!
63! 1346 2014-03-27 13:18:20Z heinze
64! Bugfix: REAL constants provided with KIND-attribute especially in call of
65! intrinsic function like MAX, MIN, SIGN
66!
67! 1322 2014-03-20 16:38:49Z raasch
68! REAL constants defined as wp_kind
69!
70! 1320 2014-03-20 08:40:49Z
71! ONLY-attribute added to USE-statements,
72! kind-parameters added to all INTEGER and REAL declaration statements,
73! kinds are defined in new module kinds,
74! revision history before 2012 removed,
75! comment fields (!:) to be used for variable explanations added to
76! all variable declaration statements
77!
78! 1092 2013-02-02 11:24:22Z raasch
79! unused variables removed
80!
81! 1071 2012-11-29 16:54:55Z franke
82! Bugfix: collision efficiencies for Hall kernel should not be < 1.0E-20
83!
84! 1036 2012-10-22 13:43:42Z raasch
85! code put under GPL (PALM 3.9)
86!
87! 1007 2012-09-19 14:30:36Z franke
88! converted all units to SI units and replaced some parameters by corresponding
89! PALM parameters
90! Bugfix: factor in calculation of enhancement factor for collision efficencies
91! changed from 10. to 1.0
92!
93! 849 2012-03-15 10:35:09Z raasch
94! routine collision_efficiency_rogers added (moved from former advec_particles
95! to here)
96!
97! 835 2012-02-22 11:21:19Z raasch $
98! Bugfix: array diss can be used only in case of Wang kernel
99!
100! 828 2012-02-21 12:00:36Z raasch
101! code has been completely reformatted, routine colker renamed
102! recalculate_kernel,
103! routine init_kernels added, radius is now communicated to the collision
104! routines by array radclass
105!
106! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
107!
108! 825 2012-02-19 03:03:44Z raasch
109! routine renamed from wang_kernel to lpm_collision_kernels,
110! turbulence_effects on collision replaced by wang_kernel
111!
112! 790 2011-11-29 03:11:20Z raasch
113! initial revision
114!
115! Description:
116! ------------
117!> This module calculates collision efficiencies either due to pure gravitational
118!> effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or
119!> including the effects of turbulence (Wang kernel, see Wang and
120!> Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8, and Ayala et al., 2008:
121!> New J. Phys., 10, 075016). The original code has been
122!> provided by L.-P. Wang but is substantially reformatted and speed optimized
123!> here.
124!------------------------------------------------------------------------------!
125 MODULE lpm_collision_kernels_mod
126 
127
128    USE constants,                                                             &
129        ONLY:  pi
130       
131    USE kinds
132
133    USE particle_attributes,                                                   &
134        ONLY:  collision_kernel, dissipation_classes, particles,               &
135               radius_classes
136
137    USE pegrid
138
139
140    IMPLICIT NONE
141
142    PRIVATE
143
144    PUBLIC  ckernel, init_kernels, rclass_lbound, rclass_ubound,               &
145            recalculate_kernel
146
147    REAL(wp) ::  epsilon       !<
148    REAL(wp) ::  rclass_lbound !<
149    REAL(wp) ::  rclass_ubound !<
150    REAL(wp) ::  urms          !<
151
152    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  epsclass  !< dissipation rate class
153    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  radclass  !< radius class
154    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  winf      !<
155   
156    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ec        !<
157    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ecf       !<
158    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gck       !<
159    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hkernel   !<
160    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hwratio   !<
161   
162    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ckernel !<
163
164    SAVE
165
166!
167!-- Public interfaces
168    INTERFACE init_kernels
169       MODULE PROCEDURE init_kernels
170    END INTERFACE init_kernels
171
172    INTERFACE recalculate_kernel
173       MODULE PROCEDURE recalculate_kernel
174    END INTERFACE recalculate_kernel
175
176
177    CONTAINS
178
179
180!------------------------------------------------------------------------------!
181! Description:
182! ------------
183!> Initialization of the collision efficiency matrix with fixed radius and
184!> dissipation classes, calculated at simulation start only.
185!------------------------------------------------------------------------------!
186 
187    SUBROUTINE init_kernels
188
189       IMPLICIT NONE
190
191       INTEGER(iwp) ::  i !<
192       INTEGER(iwp) ::  j !<
193       INTEGER(iwp) ::  k !<
194
195
196!
197!--    Calculate collision efficiencies for fixed radius- and dissipation
198!--    classes
199       IF ( collision_kernel(6:9) == 'fast' )  THEN
200
201          ALLOCATE( ckernel(1:radius_classes,1:radius_classes,                 &
202                    0:dissipation_classes), epsclass(1:dissipation_classes),   &
203                    radclass(1:radius_classes) )
204
205!
206!--       Calculate the radius class bounds with logarithmic distances
207!--       in the interval [1.0E-6, 1000.0E-6] m
208          rclass_lbound = LOG( 1.0E-6_wp )
209          rclass_ubound = LOG( 1000.0E-6_wp )
210          radclass(1)   = EXP( rclass_lbound )
211          DO  i = 2, radius_classes
212             radclass(i) = EXP( rclass_lbound +                                &
213                                ( rclass_ubound - rclass_lbound ) *            &
214                                ( i - 1.0_wp ) / ( radius_classes - 1.0_wp ) )
215          ENDDO
216
217!
218!--       Set the class bounds for dissipation in interval [0.0, 600.0] cm**2/s**3
219          DO  i = 1, dissipation_classes
220             epsclass(i) = 0.06_wp * REAL( i, KIND=wp ) / dissipation_classes
221          ENDDO
222!
223!--       Calculate collision efficiencies of the Wang/ayala kernel
224          ALLOCATE( ec(1:radius_classes,1:radius_classes),  &
225                    ecf(1:radius_classes,1:radius_classes), &
226                    gck(1:radius_classes,1:radius_classes), &
227                    winf(1:radius_classes) )
228
229          DO  k = 1, dissipation_classes
230
231             epsilon = epsclass(k)
232             urms    = 2.02_wp * ( epsilon / 0.04_wp )**( 1.0_wp / 3.0_wp )
233
234             CALL turbsd
235             CALL turb_enhance_eff
236             CALL effic
237
238             DO  j = 1, radius_classes
239                DO  i = 1, radius_classes
240                   ckernel(i,j,k) = ec(i,j) * gck(i,j) * ecf(i,j)
241                ENDDO
242             ENDDO
243
244          ENDDO
245
246!
247!--       Calculate collision efficiencies of the Hall kernel
248          ALLOCATE( hkernel(1:radius_classes,1:radius_classes), &
249                    hwratio(1:radius_classes,1:radius_classes) )
250
251          CALL fallg
252          CALL effic
253
254          DO  j = 1, radius_classes
255             DO  i =  1, radius_classes
256                hkernel(i,j) = pi * ( radclass(j) + radclass(i) )**2 &
257                                  * ec(i,j) * ABS( winf(j) - winf(i) )
258                ckernel(i,j,0) = hkernel(i,j)  ! hall kernel stored on index 0
259              ENDDO
260          ENDDO
261
262!
263!--       Test output of efficiencies
264          IF ( j == -1 )  THEN
265
266             PRINT*, '*** Hall kernel'
267             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6_wp, &
268                                              i = 1,radius_classes )
269             DO  j = 1, radius_classes
270                WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j),  &
271                                          ( hkernel(i,j), i = 1,radius_classes )
272             ENDDO
273
274             DO  k = 1, dissipation_classes
275                DO  i = 1, radius_classes
276                   DO  j = 1, radius_classes
277                      IF ( hkernel(i,j) == 0.0_wp )  THEN
278                         hwratio(i,j) = 9999999.9_wp
279                      ELSE
280                         hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j)
281                      ENDIF
282                   ENDDO
283                ENDDO
284
285                PRINT*, '*** epsilon = ', epsclass(k)
286                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i) * 1.0E6_wp, &
287                                                 i = 1,radius_classes )
288                DO  j = 1, radius_classes
289                   WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j) * 1.0E6_wp, &
290                                          ( hwratio(i,j), i = 1,radius_classes )
291                ENDDO
292             ENDDO
293
294          ENDIF
295
296          DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf )
297
298       ENDIF
299
300    END SUBROUTINE init_kernels
301
302
303!------------------------------------------------------------------------------!
304! Description:
305! ------------
306!> Calculation of collision kernels during each timestep and for each grid box
307!------------------------------------------------------------------------------!
308    SUBROUTINE recalculate_kernel( i1, j1, k1 )
309
310       USE arrays_3d,                                                          &
311           ONLY:  diss
312
313       USE particle_attributes,                                                &
314           ONLY:  number_of_particles, prt_count, radius_classes, wang_kernel
315
316       IMPLICIT NONE
317
318       INTEGER(iwp) ::  i      !<
319       INTEGER(iwp) ::  i1     !<
320       INTEGER(iwp) ::  j      !<
321       INTEGER(iwp) ::  j1     !<
322       INTEGER(iwp) ::  k1     !<
323
324
325       number_of_particles = prt_count(k1,j1,i1)
326       radius_classes      = number_of_particles   ! necessary to use the same
327                                                   ! subroutines as for
328                                                   ! precalculated kernels
329
330       ALLOCATE( ec(1:number_of_particles,1:number_of_particles), &
331                 radclass(1:number_of_particles), winf(1:number_of_particles) )
332
333!
334!--    Store particle radii on the radclass array
335       radclass(1:number_of_particles) = particles(1:number_of_particles)%radius
336
337       IF ( wang_kernel )  THEN
338          epsilon = diss(k1,j1,i1)   ! dissipation rate in m**2/s**3
339       ELSE
340          epsilon = 0.0_wp
341       ENDIF
342       urms    = 2.02_wp * ( epsilon / 0.04_wp )**( 0.33333333333_wp )
343
344       IF ( wang_kernel  .AND.  epsilon > 1.0E-7_wp )  THEN
345!
346!--       Call routines to calculate efficiencies for the Wang kernel
347          ALLOCATE( gck(1:number_of_particles,1:number_of_particles), &
348                    ecf(1:number_of_particles,1:number_of_particles) )
349
350          CALL turbsd
351          CALL turb_enhance_eff
352          CALL effic
353
354          DO  j = 1, number_of_particles
355             DO  i =  1, number_of_particles
356                ckernel(1+i-1,1+j-1,1) = ec(i,j) * gck(i,j) * ecf(i,j)
357             ENDDO
358          ENDDO
359
360          DEALLOCATE( gck, ecf )
361
362       ELSE
363!
364!--       Call routines to calculate efficiencies for the Hall kernel
365          CALL fallg
366          CALL effic
367
368          DO  j = 1, number_of_particles
369             DO  i =  1, number_of_particles
370                ckernel(i,j,1) = pi * ( radclass(j) + radclass(i) )**2         &
371                                    * ec(i,j) * ABS( winf(j) - winf(i) )
372             ENDDO
373          ENDDO
374
375       ENDIF
376
377       DEALLOCATE( ec, radclass, winf )
378
379    END SUBROUTINE recalculate_kernel
380
381
382!------------------------------------------------------------------------------!
383! Description:
384! ------------
385!> Calculation of effects of turbulence on the geometric collision kernel
386!> (by including the droplets' average radial relative velocities and their
387!> radial distribution function) following the analytic model by Aayala et al.
388!> (2008, New J. Phys.). For details check the second part 2 of the publication,
389!> page 37ff.
390!>
391!> Input parameters, which need to be replaced by PALM parameters:
392!>    water density, air density
393!------------------------------------------------------------------------------!
394    SUBROUTINE turbsd
395
396       USE control_parameters,                                                 &
397           ONLY:  g, molecular_viscosity
398   
399       USE particle_attributes,                                                &
400           ONLY:  radius_classes
401
402       IMPLICIT NONE
403
404       INTEGER(iwp) ::  i     !<
405       INTEGER(iwp) ::  j     !<
406
407       REAL(wp) ::  ao        !<
408       REAL(wp) ::  ao_gr     !<
409       REAL(wp) ::  bbb       !<
410       REAL(wp) ::  be        !<
411       REAL(wp) ::  b1        !<
412       REAL(wp) ::  b2        !<
413       REAL(wp) ::  ccc       !<
414       REAL(wp) ::  c1        !<
415       REAL(wp) ::  c1_gr     !<
416       REAL(wp) ::  c2        !<
417       REAL(wp) ::  d1        !<
418       REAL(wp) ::  d2        !<
419       REAL(wp) ::  eta       !<
420       REAL(wp) ::  e1        !<
421       REAL(wp) ::  e2        !<
422       REAL(wp) ::  fao_gr    !<
423       REAL(wp) ::  fr        !<
424       REAL(wp) ::  grfin     !<
425       REAL(wp) ::  lambda    !<
426       REAL(wp) ::  lambda_re !<
427       REAL(wp) ::  lf        !<
428       REAL(wp) ::  rc        !<
429       REAL(wp) ::  rrp       !<
430       REAL(wp) ::  sst       !<
431       REAL(wp) ::  tauk      !<
432       REAL(wp) ::  tl        !<
433       REAL(wp) ::  t2        !<
434       REAL(wp) ::  tt        !<
435       REAL(wp) ::  t1        !<
436       REAL(wp) ::  vk        !<
437       REAL(wp) ::  vrms1xy   !<
438       REAL(wp) ::  vrms2xy   !<
439       REAL(wp) ::  v1        !<
440       REAL(wp) ::  v1v2xy    !<
441       REAL(wp) ::  v1xysq    !<
442       REAL(wp) ::  v2        !<
443       REAL(wp) ::  v2xysq    !<
444       REAL(wp) ::  wrfin     !<
445       REAL(wp) ::  wrgrav2   !<
446       REAL(wp) ::  wrtur2xy  !<
447       REAL(wp) ::  xx        !<
448       REAL(wp) ::  yy        !<
449       REAL(wp) ::  z         !<
450
451       REAL(wp), DIMENSION(1:radius_classes) ::  st  !< Stokes number
452       REAL(wp), DIMENSION(1:radius_classes) ::  tau !< inertial time scale
453
454       lambda    = urms * SQRT( 15.0_wp * molecular_viscosity / epsilon )
455       lambda_re = urms**2 * SQRT( 15.0_wp / epsilon / molecular_viscosity )
456       tl        = urms**2 / epsilon
457       lf        = 0.5_wp * urms**3 / epsilon
458       tauk      = SQRT( molecular_viscosity / epsilon )
459       eta       = ( molecular_viscosity**3 / epsilon )**0.25_wp
460       vk        = eta / tauk
461
462       ao = ( 11.0_wp + 7.0_wp * lambda_re ) / ( 205.0_wp + lambda_re )
463       tt = SQRT( 2.0_wp * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk
464
465!
466!--    Get terminal velocity of droplets
467       CALL fallg
468
469       DO  i = 1, radius_classes
470          tau(i) = winf(i) / g    ! inertial time scale
471          st(i)  = tau(i) / tauk  ! Stokes number
472       ENDDO
473
474!
475!--    Calculate average radial relative velocity at contact (wrfin)
476       z   = tt / tl
477       be  = SQRT( 2.0_wp ) * lambda / lf
478       bbb = SQRT( 1.0_wp - 2.0_wp * be**2 )
479       d1  = ( 1.0_wp + bbb ) / ( 2.0_wp * bbb )
480       e1  = lf * ( 1.0_wp + bbb ) * 0.5_wp
481       d2  = ( 1.0_wp - bbb ) * 0.5_wp / bbb
482       e2  = lf * ( 1.0_wp - bbb ) * 0.5_wp
483       ccc = SQRT( 1.0_wp - 2.0_wp * z**2 )
484       b1  = ( 1.0_wp + ccc ) * 0.5_wp / ccc
485       c1  = tl * ( 1.0_wp + ccc ) * 0.5_wp
486       b2  = ( 1.0_wp - ccc ) * 0.5_wp / ccc
487       c2  = tl * ( 1.0_wp - ccc ) * 0.5_wp
488
489       DO  i = 1, radius_classes
490
491          v1 = winf(i)
492          t1 = tau(i)
493
494          DO  j = 1, i
495             rrp = radclass(i) + radclass(j)
496             v2  = winf(j)
497             t2  = tau(j)
498
499             v1xysq  = b1 * d1 * phi_w(c1,e1,v1,t1) - b1 * d2 * phi_w(c1,e2,v1,t1) &
500                     - b2 * d1 * phi_w(c2,e1,v1,t1) + b2 * d2 * phi_w(c2,e2,v1,t1)
501             v1xysq  = v1xysq * urms**2 / t1
502             vrms1xy = SQRT( v1xysq )
503
504             v2xysq  = b1 * d1 * phi_w(c1,e1,v2,t2) - b1 * d2 * phi_w(c1,e2,v2,t2) &
505                     - b2 * d1 * phi_w(c2,e1,v2,t2) + b2 * d2 * phi_w(c2,e2,v2,t2)
506             v2xysq  = v2xysq * urms**2 / t2
507             vrms2xy = SQRT( v2xysq )
508
509             IF ( winf(i) >= winf(j) )  THEN
510                v1 = winf(i)
511                t1 = tau(i)
512                v2 = winf(j)
513                t2 = tau(j)
514             ELSE
515                v1 = winf(j)
516                t1 = tau(j)
517                v2 = winf(i)
518                t2 = tau(i)
519             ENDIF
520
521             v1v2xy   =  b1 * d1 * zhi(c1,e1,v1,t1,v2,t2) - &
522                         b1 * d2 * zhi(c1,e2,v1,t1,v2,t2) - &
523                         b2 * d1 * zhi(c2,e1,v1,t1,v2,t2) + &
524                         b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
525             fr       = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
526             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)
527             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0_wp * v1v2xy
528             IF ( wrtur2xy < 0.0_wp )  wrtur2xy = 0.0_wp
529             wrgrav2  = pi / 8.0_wp * ( winf(j) - winf(i) )**2
530             wrfin    = SQRT( ( 2.0_wp / pi ) * ( wrtur2xy + wrgrav2) )
531
532!
533!--          Calculate radial distribution function (grfin)
534             IF ( st(j) > st(i) )  THEN
535                sst = st(j)
536             ELSE
537                sst = st(i)
538             ENDIF
539
540             xx = -0.1988_wp * sst**4 + 1.5275_wp * sst**3 - 4.2942_wp *       &
541                   sst**2 + 5.3406_wp * sst
542             IF ( xx < 0.0_wp )  xx = 0.0_wp
543             yy = 0.1886_wp * EXP( 20.306_wp / lambda_re )
544
545             c1_gr  =  xx / ( g / vk * tauk )**yy
546
547             ao_gr  = ao + ( pi / 8.0_wp) * ( g / vk * tauk )**2
548             fao_gr = 20.115_wp * SQRT( ao_gr / lambda_re )
549             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta
550
551             grfin  = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5_wp )
552             IF ( grfin < 1.0_wp )  grfin = 1.0_wp
553
554!
555!--          Calculate general collection kernel (without the consideration of
556!--          collection efficiencies)
557             gck(i,j) = 2.0_wp * pi * rrp**2 * wrfin * grfin
558             gck(j,i) = gck(i,j)
559
560          ENDDO
561       ENDDO
562
563    END SUBROUTINE turbsd
564
565    REAL(wp) FUNCTION phi_w( a, b, vsett, tau0 )
566!
567!--    Function used in the Ayala et al. (2008) analytical model for turbulent
568!--    effects on the collision kernel
569       IMPLICIT NONE
570
571       REAL(wp) ::  a     !<
572       REAL(wp) ::  aa1   !<
573       REAL(wp) ::  b     !<
574       REAL(wp) ::  tau0  !<
575       REAL(wp) ::  vsett !<
576
577       aa1 = 1.0_wp / tau0 + 1.0_wp / a + vsett / b
578       phi_w = 1.0_wp / aa1  - 0.5_wp * vsett / b / aa1**2
579
580    END FUNCTION phi_w
581
582    REAL(wp) FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
583!
584!--    Function used in the Ayala et al. (2008) analytical model for turbulent
585!--    effects on the collision kernel
586       IMPLICIT NONE
587
588       REAL(wp) ::  a      !<
589       REAL(wp) ::  aa1    !<
590       REAL(wp) ::  aa2    !<
591       REAL(wp) ::  aa3    !<
592       REAL(wp) ::  aa4    !<
593       REAL(wp) ::  aa5    !<
594       REAL(wp) ::  aa6    !<
595       REAL(wp) ::  b      !<
596       REAL(wp) ::  tau1   !<
597       REAL(wp) ::  tau2   !<
598       REAL(wp) ::  vsett1 !<
599       REAL(wp) ::  vsett2 !<
600
601       aa1 = vsett2 / b - 1.0_wp / tau2 - 1.0_wp / a
602       aa2 = vsett1 / b + 1.0_wp / tau1 + 1.0_wp / a
603       aa3 = ( vsett1 - vsett2 ) / b + 1.0_wp / tau1 + 1.0_wp / tau2
604       aa4 = ( vsett2 / b )**2 - ( 1.0_wp / tau2 + 1.0_wp / a )**2
605       aa5 = vsett2 / b + 1.0_wp / tau2 + 1.0_wp / a
606       aa6 = 1.0_wp / tau1 - 1.0_wp / a + ( 1.0_wp / tau2 + 1.0_wp / a) *      &
607             vsett1 / vsett2
608       zhi = (1.0_wp / aa1 - 1.0_wp / aa2 ) * ( vsett1 - vsett2 ) * 0.5_wp /   &
609             b / aa3**2 + ( 4.0_wp / aa4 - 1.0_wp / aa5**2 - 1.0_wp / aa1**2 ) &
610             * vsett2 * 0.5_wp / b /aa6 + ( 2.0_wp * ( b / aa2 - b / aa1 ) -   &
611             vsett1 / aa2**2 + vsett2 / aa1**2 ) * 0.5_wp / b / aa3
612
613    END FUNCTION zhi
614
615
616!------------------------------------------------------------------------------!
617! Description:
618! ------------
619!> Parameterization of terminal velocity following Rogers et al. (1993, J. Appl.
620!> Meteorol.)
621!------------------------------------------------------------------------------!
622    SUBROUTINE fallg
623
624       USE particle_attributes,                                                &
625           ONLY:  radius_classes
626
627       IMPLICIT NONE
628
629       INTEGER(iwp) ::  j                            !<
630
631       REAL(wp), PARAMETER ::  k_cap_rog = 4.0_wp    !< parameter
632       REAL(wp), PARAMETER ::  k_low_rog = 12.0_wp   !< parameter
633       REAL(wp), PARAMETER ::  a_rog     = 9.65_wp   !< parameter
634       REAL(wp), PARAMETER ::  b_rog     = 10.43_wp  !< parameter
635       REAL(wp), PARAMETER ::  c_rog     = 0.6_wp    !< parameter
636       REAL(wp), PARAMETER ::  d0_rog    = 0.745_wp  !< seperation diameter
637
638       REAL(wp)            ::  diameter              !< droplet diameter in mm
639
640
641       DO  j = 1, radius_classes
642
643          diameter = radclass(j) * 2000.0_wp
644
645          IF ( diameter <= d0_rog )  THEN
646             winf(j) = k_cap_rog * diameter * ( 1.0_wp -                       &
647                                                EXP( -k_low_rog * diameter ) )
648          ELSE
649             winf(j) = a_rog - b_rog * EXP( -c_rog * diameter )
650          ENDIF
651
652       ENDDO
653
654    END SUBROUTINE fallg
655
656
657!------------------------------------------------------------------------------!
658! Description:
659! ------------
660!> Interpolation of collision efficiencies (Hall, 1980, J. Atmos. Sci.)
661!------------------------------------------------------------------------------!
662    SUBROUTINE effic
663 
664       USE particle_attributes,                                                &
665           ONLY:  radius_classes
666
667       IMPLICIT NONE
668
669       INTEGER(iwp) ::  i  !<
670       INTEGER(iwp) ::  iq !<
671       INTEGER(iwp) ::  ir !<
672       INTEGER(iwp) ::  j  !<
673       INTEGER(iwp) ::  k  !<
674
675       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ira !<
676
677       LOGICAL, SAVE ::  first = .TRUE. !<
678
679       REAL(wp) ::  ek              !<
680       REAL(wp) ::  particle_radius !<
681       REAL(wp) ::  pp              !<
682       REAL(wp) ::  qq              !<
683       REAL(wp) ::  rq              !<
684
685       REAL(wp), DIMENSION(1:21), SAVE ::  rat        !<
686       
687       REAL(wp), DIMENSION(1:15), SAVE ::  r0         !<
688       
689       REAL(wp), DIMENSION(1:15,1:21), SAVE ::  ecoll !<
690
691!
692!--    Initial assignment of constants
693       IF ( first )  THEN
694
695         first = .FALSE.
696         r0  = (/   6.0_wp,   8.0_wp,  10.0_wp, 15.0_wp,  20.0_wp,  25.0_wp,   &
697                   30.0_wp,  40.0_wp,  50.0_wp, 60.0_wp,  70.0_wp, 100.0_wp,   &
698                  150.0_wp, 200.0_wp, 300.0_wp /)
699
700         rat = (/ 0.00_wp, 0.05_wp, 0.10_wp, 0.15_wp, 0.20_wp, 0.25_wp,        &
701                  0.30_wp, 0.35_wp, 0.40_wp, 0.45_wp, 0.50_wp, 0.55_wp,        &
702                  0.60_wp, 0.65_wp, 0.70_wp, 0.75_wp, 0.80_wp, 0.85_wp,        &
703                  0.90_wp, 0.95_wp, 1.00_wp /)
704
705         ecoll(:,1)  = (/ 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp,    &
706                          0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp,    &
707                          0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp /)
708         ecoll(:,2)  = (/ 0.003_wp, 0.003_wp, 0.003_wp, 0.004_wp, 0.005_wp,    &
709                          0.005_wp, 0.005_wp, 0.010_wp, 0.100_wp, 0.050_wp,    &
710                          0.200_wp, 0.500_wp, 0.770_wp, 0.870_wp, 0.970_wp /)
711         ecoll(:,3)  = (/ 0.007_wp, 0.007_wp, 0.007_wp, 0.008_wp, 0.009_wp,    &
712                          0.010_wp, 0.010_wp, 0.070_wp, 0.400_wp, 0.430_wp,    &
713                          0.580_wp, 0.790_wp, 0.930_wp, 0.960_wp, 1.000_wp /)
714         ecoll(:,4)  = (/ 0.009_wp, 0.009_wp, 0.009_wp, 0.012_wp, 0.015_wp,    &
715                          0.010_wp, 0.020_wp, 0.280_wp, 0.600_wp, 0.640_wp,    &
716                          0.750_wp, 0.910_wp, 0.970_wp, 0.980_wp, 1.000_wp /)
717         ecoll(:,5)  = (/ 0.014_wp, 0.014_wp, 0.014_wp, 0.015_wp, 0.016_wp,    &
718                          0.030_wp, 0.060_wp, 0.500_wp, 0.700_wp, 0.770_wp,    &
719                          0.840_wp, 0.950_wp, 0.970_wp, 1.000_wp, 1.000_wp /)
720         ecoll(:,6)  = (/ 0.017_wp, 0.017_wp, 0.017_wp, 0.020_wp, 0.022_wp,    &
721                          0.060_wp, 0.100_wp, 0.620_wp, 0.780_wp, 0.840_wp,    &
722                          0.880_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
723         ecoll(:,7)  = (/ 0.030_wp, 0.030_wp, 0.024_wp, 0.022_wp, 0.032_wp,    &
724                          0.062_wp, 0.200_wp, 0.680_wp, 0.830_wp, 0.870_wp,    &
725                          0.900_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
726         ecoll(:,8)  = (/ 0.025_wp, 0.025_wp, 0.025_wp, 0.036_wp, 0.043_wp,    &
727                          0.130_wp, 0.270_wp, 0.740_wp, 0.860_wp, 0.890_wp,    &
728                          0.920_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
729         ecoll(:,9)  = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.040_wp, 0.052_wp,    &
730                          0.200_wp, 0.400_wp, 0.780_wp, 0.880_wp, 0.900_wp,    &
731                          0.940_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
732         ecoll(:,10) = (/ 0.030_wp, 0.030_wp, 0.030_wp, 0.047_wp, 0.064_wp,    &
733                          0.250_wp, 0.500_wp, 0.800_wp, 0.900_wp, 0.910_wp,    &
734                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
735         ecoll(:,11) = (/ 0.040_wp, 0.040_wp, 0.033_wp, 0.037_wp, 0.068_wp,    &
736                          0.240_wp, 0.550_wp, 0.800_wp, 0.900_wp, 0.910_wp,    &
737                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
738         ecoll(:,12) = (/ 0.035_wp, 0.035_wp, 0.035_wp, 0.055_wp, 0.079_wp,    &
739                          0.290_wp, 0.580_wp, 0.800_wp, 0.900_wp, 0.910_wp,    &
740                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
741         ecoll(:,13) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.062_wp, 0.082_wp,    &
742                          0.290_wp, 0.590_wp, 0.780_wp, 0.900_wp, 0.910_wp,    &
743                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
744         ecoll(:,14) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.060_wp, 0.080_wp,    &
745                          0.290_wp, 0.580_wp, 0.770_wp, 0.890_wp, 0.910_wp,    &
746                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
747         ecoll(:,15) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.041_wp, 0.075_wp,    &
748                          0.250_wp, 0.540_wp, 0.760_wp, 0.880_wp, 0.920_wp,    &
749                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
750         ecoll(:,16) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.052_wp, 0.067_wp,    &
751                          0.250_wp, 0.510_wp, 0.770_wp, 0.880_wp, 0.930_wp,    &
752                          0.970_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
753         ecoll(:,17) = (/ 0.037_wp, 0.037_wp, 0.037_wp, 0.047_wp, 0.057_wp,    &
754                          0.250_wp, 0.490_wp, 0.770_wp, 0.890_wp, 0.950_wp,    &
755                          1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
756         ecoll(:,18) = (/ 0.036_wp, 0.036_wp, 0.036_wp, 0.042_wp, 0.048_wp,    &
757                          0.230_wp, 0.470_wp, 0.780_wp, 0.920_wp, 1.000_wp,    &
758                          1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp /)
759         ecoll(:,19) = (/ 0.040_wp, 0.040_wp, 0.035_wp, 0.033_wp, 0.040_wp,    &
760                          0.112_wp, 0.450_wp, 0.790_wp, 1.010_wp, 1.030_wp,    &
761                          1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp /)
762         ecoll(:,20) = (/ 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp, 0.033_wp,    &
763                          0.119_wp, 0.470_wp, 0.950_wp, 1.300_wp, 1.700_wp,    &
764                          2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp /)
765         ecoll(:,21) = (/ 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp, 0.027_wp,    &
766                          0.125_wp, 0.520_wp, 1.400_wp, 2.300_wp, 3.000_wp,    &
767                          4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp /)
768       ENDIF
769
770!
771!--    Calculate the radius class index of particles with respect to array r
772!--    Radius has to be in microns
773       ALLOCATE( ira(1:radius_classes) )
774       DO  j = 1, radius_classes
775          particle_radius = radclass(j) * 1.0E6_wp
776          DO  k = 1, 15
777             IF ( particle_radius < r0(k) )  THEN
778                ira(j) = k
779                EXIT
780             ENDIF
781          ENDDO
782          IF ( particle_radius >= r0(15) )  ira(j) = 16
783       ENDDO
784
785!
786!--    Two-dimensional linear interpolation of the collision efficiency.
787!--    Radius has to be in microns
788       DO  j = 1, radius_classes
789          DO  i = 1, j
790
791             ir = MAX( ira(i), ira(j) )
792             rq = MIN( radclass(i) / radclass(j), radclass(j) / radclass(i) )
793             iq = INT( rq * 20 ) + 1
794             iq = MAX( iq , 2)
795
796             IF ( ir < 16 )  THEN
797                IF ( ir >= 2 )  THEN
798                   pp = ( ( MAX( radclass(j), radclass(i) ) * 1.0E6_wp ) -     &
799                          r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
800                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
801                   ec(j,i) = ( 1.0_wp - pp ) * ( 1.0_wp - qq )                 &
802                             * ecoll(ir-1,iq-1)                                &
803                             + pp * ( 1.0_wp - qq ) * ecoll(ir,iq-1)           &
804                             + qq * ( 1.0_wp - pp ) * ecoll(ir-1,iq)           &
805                             + pp * qq * ecoll(ir,iq)
806                ELSE
807                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
808                   ec(j,i) = ( 1.0_wp - qq ) * ecoll(1,iq-1) + qq * ecoll(1,iq)
809                ENDIF
810             ELSE
811                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
812                ek = ( 1.0_wp - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
813                ec(j,i) = MIN( ek, 1.0_wp )
814             ENDIF
815
816             IF ( ec(j,i) < 1.0E-20_wp )  ec(j,i) = 0.0_wp
817
818             ec(i,j) = ec(j,i)
819
820          ENDDO
821       ENDDO
822
823       DEALLOCATE( ira )
824
825    END SUBROUTINE effic
826
827
828!------------------------------------------------------------------------------!
829! Description:
830! ------------
831!> Interpolation of turbulent enhancement factor for collision efficencies
832!> following Wang and Grabowski (2009, Atmos. Sci. Let.)
833!------------------------------------------------------------------------------!
834    SUBROUTINE turb_enhance_eff
835
836       USE particle_attributes,                                                &
837           ONLY:  radius_classes
838
839       IMPLICIT NONE
840
841       INTEGER(iwp) :: i  !<
842       INTEGER(iwp) :: iq !<
843       INTEGER(iwp) :: ir !<
844       INTEGER(iwp) :: j  !<
845       INTEGER(iwp) :: k  !<
846       INTEGER(iwp) :: kk !<
847
848       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ira !<
849       
850       LOGICAL, SAVE ::  first = .TRUE. !<
851
852       REAL(wp) ::  particle_radius !<
853       REAL(wp) ::  pp              !<
854       REAL(wp) ::  qq              !<
855       REAL(wp) ::  rq              !<
856       REAL(wp) ::  y1              !<
857       REAL(wp) ::  y2              !<
858       REAL(wp) ::  y3              !<
859
860       REAL(wp), DIMENSION(1:11), SAVE ::  rat           !<
861       REAL(wp), DIMENSION(1:7), SAVE  ::  r0            !<
862       
863       REAL(wp), DIMENSION(1:7,1:11), SAVE ::  ecoll_100 !<
864       REAL(wp), DIMENSION(1:7,1:11), SAVE ::  ecoll_400 !<
865
866!
867!--    Initial assignment of constants
868       IF ( first )  THEN
869
870          first = .FALSE.
871
872          r0  = (/  10.0_wp, 20.0_wp, 30.0_wp, 40.0_wp, 50.0_wp, 60.0_wp,  &
873                   100.0_wp /)
874
875          rat = (/ 0.0_wp, 0.1_wp, 0.2_wp, 0.3_wp, 0.4_wp, 0.5_wp, 0.6_wp, &
876                   0.7_wp, 0.8_wp, 0.9_wp, 1.0_wp /)
877!
878!--       Tabulated turbulent enhancement factor at 100 cm**2/s**3
879          ecoll_100(:,1)  = (/  1.74_wp,   1.74_wp,   1.773_wp, 1.49_wp,  &
880                                1.207_wp,  1.207_wp,  1.0_wp /)
881          ecoll_100(:,2)  = (/  1.46_wp,   1.46_wp,   1.421_wp, 1.245_wp, &
882                                1.069_wp,  1.069_wp,  1.0_wp /)
883          ecoll_100(:,3)  = (/  1.32_wp,   1.32_wp,   1.245_wp, 1.123_wp, &
884                                1.000_wp,  1.000_wp,  1.0_wp /)
885          ecoll_100(:,4)  = (/  1.250_wp,  1.250_wp,  1.148_wp, 1.087_wp, &
886                                1.025_wp,  1.025_wp,  1.0_wp /)
887          ecoll_100(:,5)  = (/  1.186_wp,  1.186_wp,  1.066_wp, 1.060_wp, &
888                                1.056_wp,  1.056_wp,  1.0_wp /)
889          ecoll_100(:,6)  = (/  1.045_wp,  1.045_wp,  1.000_wp, 1.014_wp, &
890                                1.028_wp,  1.028_wp,  1.0_wp /)
891          ecoll_100(:,7)  = (/  1.070_wp,  1.070_wp,  1.030_wp, 1.038_wp, &
892                                1.046_wp,  1.046_wp,  1.0_wp /)
893          ecoll_100(:,8)  = (/  1.000_wp,  1.000_wp,  1.054_wp, 1.042_wp, &
894                                1.029_wp,  1.029_wp,  1.0_wp /)
895          ecoll_100(:,9)  = (/  1.223_wp,  1.223_wp,  1.117_wp, 1.069_wp, &
896                                1.021_wp,  1.021_wp,  1.0_wp /)
897          ecoll_100(:,10) = (/  1.570_wp,  1.570_wp,  1.244_wp, 1.166_wp, &
898                                1.088_wp,  1.088_wp,  1.0_wp /)
899          ecoll_100(:,11) = (/ 20.3_wp,   20.3_wp,   14.6_wp,   8.61_wp,  &
900                                2.60_wp,   2.60_wp,   1.0_wp /)
901!
902!--       Tabulated turbulent enhancement factor at 400 cm**2/s**3
903          ecoll_400(:,1)  = (/  4.976_wp,  4.976_wp,  3.593_wp,  2.519_wp, &
904                                1.445_wp,  1.445_wp,  1.0_wp /)
905          ecoll_400(:,2)  = (/  2.984_wp,  2.984_wp,  2.181_wp,  1.691_wp, &
906                                1.201_wp,  1.201_wp,  1.0_wp /)
907          ecoll_400(:,3)  = (/  1.988_wp,  1.988_wp,  1.475_wp,  1.313_wp, &
908                                1.150_wp,  1.150_wp,  1.0_wp /)
909          ecoll_400(:,4)  = (/  1.490_wp,  1.490_wp,  1.187_wp,  1.156_wp, &
910                                1.126_wp,  1.126_wp,  1.0_wp /)
911          ecoll_400(:,5)  = (/  1.249_wp,  1.249_wp,  1.088_wp,  1.090_wp, &
912                                1.092_wp,  1.092_wp,  1.0_wp /)
913          ecoll_400(:,6)  = (/  1.139_wp,  1.139_wp,  1.130_wp,  1.091_wp, &
914                                1.051_wp,  1.051_wp,  1.0_wp /)
915          ecoll_400(:,7)  = (/  1.220_wp,  1.220_wp,  1.190_wp,  1.138_wp, &
916                                1.086_wp,  1.086_wp,  1.0_wp /)
917          ecoll_400(:,8)  = (/  1.325_wp,  1.325_wp,  1.267_wp,  1.165_wp, &
918                                1.063_wp,  1.063_wp,  1.0_wp /)
919          ecoll_400(:,9)  = (/  1.716_wp,  1.716_wp,  1.345_wp,  1.223_wp, &
920                                1.100_wp,  1.100_wp,  1.0_wp /)
921          ecoll_400(:,10) = (/  3.788_wp,  3.788_wp,  1.501_wp,  1.311_wp, &
922                                1.120_wp,  1.120_wp,  1.0_wp /)
923          ecoll_400(:,11) = (/ 36.52_wp,  36.52_wp,  19.16_wp,  22.80_wp,  &
924                               26.0_wp,   26.0_wp,    1.0_wp /)
925
926       ENDIF
927
928!
929!--    Calculate the radius class index of particles with respect to array r0
930!--    The droplet radius has to be given in microns.
931       ALLOCATE( ira(1:radius_classes) )
932
933       DO  j = 1, radius_classes
934          particle_radius = radclass(j) * 1.0E6_wp
935          DO  k = 1, 7
936             IF ( particle_radius < r0(k) )  THEN
937                ira(j) = k
938                EXIT
939             ENDIF
940          ENDDO
941          IF ( particle_radius >= r0(7) )  ira(j) = 8
942       ENDDO
943
944!
945!--    Two-dimensional linear interpolation of the turbulent enhancement factor.
946!--    The droplet radius has to be given in microns.
947       DO  j =  1, radius_classes
948          DO  i = 1, j
949
950             ir = MAX( ira(i), ira(j) )
951             rq = MIN( radclass(i) / radclass(j), radclass(j) / radclass(i) )
952
953             DO  kk = 2, 11
954                IF ( rq <= rat(kk) )  THEN
955                   iq = kk
956                   EXIT
957                ENDIF
958             ENDDO
959
960             y1 = 1.0_wp  ! turbulent enhancement factor at 0 m**2/s**3
961
962             IF ( ir < 8 )  THEN
963                IF ( ir >= 2 )  THEN
964                   pp = ( MAX( radclass(j), radclass(i) ) * 1.0E6_wp -  &
965                          r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
966                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
967                   y2 = ( 1.0_wp - pp ) * ( 1.0_wp - qq ) * ecoll_100(ir-1,iq-1) + &
968                                pp * ( 1.0_wp - qq ) * ecoll_100(ir,iq-1)        + &
969                                qq * ( 1.0_wp - pp ) * ecoll_100(ir-1,iq)        + &
970                                pp * qq              * ecoll_100(ir,iq)
971                   y3 = ( 1.0-pp ) * ( 1.0_wp - qq ) * ecoll_400(ir-1,iq-1)      + &
972                                pp * ( 1.0_wp - qq ) * ecoll_400(ir,iq-1)        + &
973                                qq * ( 1.0_wp - pp ) * ecoll_400(ir-1,iq)        + &
974                                pp * qq              * ecoll_400(ir,iq)
975                ELSE
976                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
977                   y2 = ( 1.0_wp - qq ) * ecoll_100(1,iq-1) + qq * ecoll_100(1,iq)
978                   y3 = ( 1.0_wp - qq ) * ecoll_400(1,iq-1) + qq * ecoll_400(1,iq)
979                ENDIF
980             ELSE
981                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
982                y2 = ( 1.0_wp - qq ) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
983                y3 = ( 1.0_wp - qq ) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
984             ENDIF
985!
986!--          Linear interpolation of turbulent enhancement factor
987             IF ( epsilon <= 0.01_wp )  THEN
988                ecf(j,i) = ( epsilon - 0.01_wp ) / ( 0.0_wp  - 0.01_wp ) * y1 &
989                         + ( epsilon - 0.0_wp  ) / ( 0.01_wp - 0.0_wp  ) * y2
990             ELSEIF ( epsilon <= 0.06_wp )  THEN
991                ecf(j,i) = ( epsilon - 0.04_wp ) / ( 0.01_wp - 0.04_wp ) * y2 &
992                         + ( epsilon - 0.01_wp ) / ( 0.04_wp - 0.01_wp ) * y3
993             ELSE
994                ecf(j,i) = ( 0.06_wp - 0.04_wp ) / ( 0.01_wp - 0.04_wp ) * y2 &
995                         + ( 0.06_wp - 0.01_wp ) / ( 0.04_wp - 0.01_wp ) * y3
996             ENDIF
997
998             IF ( ecf(j,i) < 1.0_wp )  ecf(j,i) = 1.0_wp
999
1000             ecf(i,j) = ecf(j,i)
1001
1002          ENDDO
1003       ENDDO
1004
1005    END SUBROUTINE turb_enhance_eff
1006
1007 END MODULE lpm_collision_kernels_mod
Note: See TracBrowser for help on using the repository browser.