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

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

Bugfix in lpm_collision_kernels

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