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

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

revised renaming of modules

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