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

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

last commit documented

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