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

Last change on this file since 1874 was 1874, checked in by maronga, 5 years ago

last commit documented

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