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

Last change on this file since 1859 was 1859, checked in by hoffmann, 8 years ago

last commit documented

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