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

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

changes in LPM and bulk cloud microphysics

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