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

Last change on this file since 3768 was 3655, checked in by knoop, 3 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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