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

Last change on this file since 1850 was 1850, checked in by maronga, 6 years ago

added _mod string to several filenames to meet the naming convection for modules

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