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

Last change on this file since 1834 was 1823, checked in by hoffmann, 5 years ago

last commit documented

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