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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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