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

Last change on this file since 1822 was 1822, checked in by hoffmann, 7 years ago

changes in LPM and bulk cloud microphysics

  • Property svn:keywords set to Id
File size: 38.5 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! PALM kernel has been deleted.
22! Bugfix in the calculation of the turbulent enhancement factor of the
23! collection efficiency.
24!
25! Unused variables removed.
26!
27! Former revisions:
28! -----------------
29! $Id: lpm_collision_kernels.f90 1822 2016-04-07 07:49:42Z hoffmann $
30!
31! 1776 2016-03-02 17:54:58Z hoffmann
32! Bugfix: Collection efficiencies must be calculated for the larger droplet.
33!
34! 1682 2015-10-07 23:56:08Z knoop
35! Code annotations made doxygen readable
36!
37! 1519 2015-01-08 10:20:42Z hoffmann
38! Bugfix: Using the new particle structure, particles are not sorted by size.
39! Hence, computation of collision efficiencies must ensure that the ratio of
40! two colliding droplets is < 1.
41!
42! 1359 2014-04-11 17:15:14Z hoffmann
43! New particle structure integrated.
44! Kind definition added to all floating point numbers.
45!
46! 1346 2014-03-27 13:18:20Z heinze
47! Bugfix: REAL constants provided with KIND-attribute especially in call of
48! intrinsic function like MAX, MIN, SIGN
49!
50! 1322 2014-03-20 16:38:49Z raasch
51! REAL constants defined as wp_kind
52!
53! 1320 2014-03-20 08:40:49Z
54! ONLY-attribute added to USE-statements,
55! kind-parameters added to all INTEGER and REAL declaration statements,
56! kinds are defined in new module kinds,
57! revision history before 2012 removed,
58! comment fields (!:) to be used for variable explanations added to
59! all variable declaration statements
60!
61! 1092 2013-02-02 11:24:22Z raasch
62! unused variables removed
63!
64! 1071 2012-11-29 16:54:55Z franke
65! Bugfix: collision efficiencies for Hall kernel should not be < 1.0E-20
66!
67! 1036 2012-10-22 13:43:42Z raasch
68! code put under GPL (PALM 3.9)
69!
70! 1007 2012-09-19 14:30:36Z franke
71! converted all units to SI units and replaced some parameters by corresponding
72! PALM parameters
73! Bugfix: factor in calculation of enhancement factor for collision efficencies
74! changed from 10. to 1.0
75!
76! 849 2012-03-15 10:35:09Z raasch
77! routine collision_efficiency_rogers added (moved from former advec_particles
78! to here)
79!
80! 835 2012-02-22 11:21:19Z raasch $
81! Bugfix: array diss can be used only in case of Wang kernel
82!
83! 828 2012-02-21 12:00:36Z raasch
84! code has been completely reformatted, routine colker renamed
85! recalculate_kernel,
86! routine init_kernels added, radius is now communicated to the collision
87! routines by array radclass
88!
89! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
90!
91! 825 2012-02-19 03:03:44Z raasch
92! routine renamed from wang_kernel to lpm_collision_kernels,
93! turbulence_effects on collision replaced by wang_kernel
94!
95! 790 2011-11-29 03:11:20Z raasch
96! initial revision
97!
98! Description:
99! ------------
100!> This module calculates collision efficiencies either due to pure gravitational
101!> effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or
102!> including the effects of turbulence (Wang kernel, see Wang and
103!> Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8, and Ayala et al., 2008:
104!> New J. Phys., 10, 075016). The original code has been
105!> provided by L.-P. Wang but is substantially reformatted and speed optimized
106!> here.
107!------------------------------------------------------------------------------!
108 MODULE lpm_collision_kernels_mod
109 
110
111    USE constants,                                                             &
112        ONLY:  pi
113       
114    USE kinds
115
116    USE particle_attributes,                                                   &
117        ONLY:  collision_kernel, dissipation_classes, particles,               &
118               radius_classes
119
120    USE pegrid
121
122
123    IMPLICIT NONE
124
125    PRIVATE
126
127    PUBLIC  ckernel, init_kernels, rclass_lbound, rclass_ubound,               &
128            recalculate_kernel
129
130    REAL(wp) ::  epsilon       !<
131    REAL(wp) ::  rclass_lbound !<
132    REAL(wp) ::  rclass_ubound !<
133    REAL(wp) ::  urms          !<
134
135    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  epsclass  !< dissipation rate class
136    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  radclass  !< radius class
137    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  winf      !<
138   
139    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ec        !<
140    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ecf       !<
141    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gck       !<
142    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hkernel   !<
143    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  hwratio   !<
144   
145    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ckernel !<
146
147    SAVE
148
149!
150!-- Public interfaces
151    INTERFACE init_kernels
152       MODULE PROCEDURE init_kernels
153    END INTERFACE init_kernels
154
155    INTERFACE recalculate_kernel
156       MODULE PROCEDURE recalculate_kernel
157    END INTERFACE recalculate_kernel
158
159
160    CONTAINS
161
162
163!------------------------------------------------------------------------------!
164! Description:
165! ------------
166!> Initialization of the collision efficiency matrix with fixed radius and
167!> dissipation classes, calculated at simulation start only.
168!------------------------------------------------------------------------------!
169 
170    SUBROUTINE init_kernels
171
172       IMPLICIT NONE
173
174       INTEGER(iwp) ::  i !<
175       INTEGER(iwp) ::  j !<
176       INTEGER(iwp) ::  k !<
177
178
179!
180!--    Calculate collision efficiencies for fixed radius- and dissipation
181!--    classes
182       IF ( collision_kernel(6:9) == 'fast' )  THEN
183
184          ALLOCATE( ckernel(1:radius_classes,1:radius_classes,                 &
185                    0:dissipation_classes), epsclass(1:dissipation_classes),   &
186                    radclass(1:radius_classes) )
187
188!
189!--       Calculate the radius class bounds with logarithmic distances
190!--       in the interval [1.0E-6, 2.0E-4] m
191          rclass_lbound = LOG( 1.0E-6_wp )
192          rclass_ubound = LOG( 2.0E-4_wp )
193          radclass(1)   = EXP( rclass_lbound )
194          DO  i = 2, radius_classes
195             radclass(i) = EXP( rclass_lbound +                                &
196                                ( rclass_ubound - rclass_lbound ) *            &
197                                ( i - 1.0_wp ) / ( radius_classes - 1.0_wp ) )
198          ENDDO
199
200!
201!--       Set the class bounds for dissipation in interval [0.0, 0.1] m**2/s**3
202          DO  i = 1, dissipation_classes
203             epsclass(i) = 0.1_wp * REAL( i, KIND=wp ) / dissipation_classes
204          ENDDO
205!
206!--       Calculate collision efficiencies of the Wang/ayala kernel
207          ALLOCATE( ec(1:radius_classes,1:radius_classes),  &
208                    ecf(1:radius_classes,1:radius_classes), &
209                    gck(1:radius_classes,1:radius_classes), &
210                    winf(1:radius_classes) )
211
212          DO  k = 1, dissipation_classes
213
214             epsilon = epsclass(k)
215             urms    = 2.02_wp * ( epsilon / 0.04_wp )**( 1.0_wp / 3.0_wp )
216
217             CALL turbsd
218             CALL turb_enhance_eff
219             CALL effic
220
221             DO  j = 1, radius_classes
222                DO  i = 1, radius_classes
223                   ckernel(i,j,k) = ec(i,j) * gck(i,j) * ecf(i,j)
224                ENDDO
225             ENDDO
226
227          ENDDO
228
229!
230!--       Calculate collision efficiencies of the Hall kernel
231          ALLOCATE( hkernel(1:radius_classes,1:radius_classes), &
232                    hwratio(1:radius_classes,1:radius_classes) )
233
234          CALL fallg
235          CALL effic
236
237          DO  j = 1, radius_classes
238             DO  i =  1, radius_classes
239                hkernel(i,j) = pi * ( radclass(j) + radclass(i) )**2 &
240                                  * ec(i,j) * ABS( winf(j) - winf(i) )
241                ckernel(i,j,0) = hkernel(i,j)  ! hall kernel stored on index 0
242              ENDDO
243          ENDDO
244
245!
246!--       Test output of efficiencies
247          IF ( j == -1 )  THEN
248
249             PRINT*, '*** Hall kernel'
250             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6_wp, &
251                                              i = 1,radius_classes )
252             DO  j = 1, radius_classes
253                WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j),  &
254                                          ( hkernel(i,j), i = 1,radius_classes )
255             ENDDO
256
257             DO  k = 1, dissipation_classes
258                DO  i = 1, radius_classes
259                   DO  j = 1, radius_classes
260                      IF ( hkernel(i,j) == 0.0_wp )  THEN
261                         hwratio(i,j) = 9999999.9_wp
262                      ELSE
263                         hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j)
264                      ENDIF
265                   ENDDO
266                ENDDO
267
268                PRINT*, '*** epsilon = ', epsclass(k)
269                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i) * 1.0E6_wp, &
270                                                 i = 1,radius_classes )
271                DO  j = 1, radius_classes
272                   WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j) * 1.0E6_wp, &
273                                          ( hwratio(i,j), i = 1,radius_classes )
274                ENDDO
275             ENDDO
276
277          ENDIF
278
279          DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf )
280
281       ELSEIF( collision_kernel == 'hall'  .OR.  collision_kernel == 'wang' ) &
282       THEN
283!
284!--       Initial settings for Hall- and Wang-Kernel
285!--       To be done: move here parts from turbsd, fallg, ecoll, etc.
286       ENDIF
287
288    END SUBROUTINE init_kernels
289
290
291!------------------------------------------------------------------------------!
292! Description:
293! ------------
294!> Calculation of collision kernels during each timestep and for each grid box
295!------------------------------------------------------------------------------!
296    SUBROUTINE recalculate_kernel( i1, j1, k1 )
297
298       USE arrays_3d,                                                          &
299           ONLY:  diss
300
301       USE particle_attributes,                                                &
302           ONLY:  prt_count, radius_classes, wang_kernel
303
304       IMPLICIT NONE
305
306       INTEGER(iwp) ::  i      !<
307       INTEGER(iwp) ::  i1     !<
308       INTEGER(iwp) ::  j      !<
309       INTEGER(iwp) ::  j1     !<
310       INTEGER(iwp) ::  k1     !<
311       INTEGER(iwp) ::  pend   !<
312       INTEGER(iwp) ::  pstart !<
313
314
315       pstart = 1
316       pend   = prt_count(k1,j1,i1)
317       radius_classes = prt_count(k1,j1,i1)
318
319       ALLOCATE( ec(1:radius_classes,1:radius_classes), &
320                 radclass(1:radius_classes), winf(1:radius_classes) )
321
322!
323!--    Store particle radii on the radclass array
324       radclass(1:radius_classes) = particles(pstart:pend)%radius
325
326       IF ( wang_kernel )  THEN
327          epsilon = diss(k1,j1,i1)   ! dissipation rate in m**2/s**3
328       ELSE
329          epsilon = 0.0_wp
330       ENDIF
331       urms    = 2.02_wp * ( epsilon / 0.04_wp )**( 0.33333333333_wp )
332
333       IF ( wang_kernel  .AND.  epsilon > 1.0E-7_wp )  THEN
334!
335!--       Call routines to calculate efficiencies for the Wang kernel
336          ALLOCATE( gck(1:radius_classes,1:radius_classes), &
337                    ecf(1:radius_classes,1:radius_classes) )
338
339          CALL turbsd
340          CALL turb_enhance_eff
341          CALL effic
342
343          DO  j = 1, radius_classes
344             DO  i =  1, radius_classes
345                ckernel(pstart+i-1,pstart+j-1,1) = ec(i,j) * gck(i,j) * ecf(i,j)
346             ENDDO
347          ENDDO
348
349          DEALLOCATE( gck, ecf )
350
351       ELSE
352!
353!--       Call routines to calculate efficiencies for the Hall kernel
354          CALL fallg
355          CALL effic
356
357          DO  j = 1, radius_classes
358             DO  i =  1, radius_classes
359                ckernel(pstart+i-1,pstart+j-1,1) = pi *                       &
360                                          ( 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.