source: palm/trunk/SOURCE/lpm_collision_kernels_mod.f90 @ 1858

Last change on this file since 1858 was 1858, checked in by hoffmann, 6 years ago

Interpolation of collision kernels adjusted to more reasonable values

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