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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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