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

Last change on this file since 3352 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

  • Property svn:keywords set to Id
File size: 39.2 KB
RevLine 
[1873]1!> @file lpm_collision_kernels.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]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.
[1036]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!
[2718]17! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[790]20! Current revisions:
21! -----------------
[1347]22!
[2001]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: lpm_collision_kernels.f90 3274 2018-09-24 15:42:55Z knoop $
[3274]27! Modularization of all bulk cloud physics code components
28!
29! 2718 2018-01-02 08:49:38Z maronga
[2716]30! Corrected "Former revisions" section
31!
32! 2696 2017-12-14 17:12:51Z kanani
33! Change in file header (GPL part)
[1321]34!
[2716]35! 2101 2017-01-05 16:42:31Z suehring
36!
[2001]37! 2000 2016-08-20 18:09:15Z knoop
38! Forced header and separation lines into 80 columns
39!
[1881]40! 1880 2016-04-20 09:36:50Z hoffmann
41! Bugfix: The index of the larger particle has to be chosen for interpolation.
42!
[1874]43! 1873 2016-04-18 14:50:06Z maronga
44! Module renamed (removed _mod)
45!
[1859]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!
[1851]50! 1850 2016-04-08 13:29:27Z maronga
51! Module renamed
52!
[1823]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!
[1777]60! 1776 2016-03-02 17:54:58Z hoffmann
61! Bugfix: Collection efficiencies must be calculated for the larger droplet.
62!
[1683]63! 1682 2015-10-07 23:56:08Z knoop
64! Code annotations made doxygen readable
65!
[1520]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!
[1360]71! 1359 2014-04-11 17:15:14Z hoffmann
72! New particle structure integrated.
73! Kind definition added to all floating point numbers.
74!
[1347]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!
[1323]79! 1322 2014-03-20 16:38:49Z raasch
80! REAL constants defined as wp_kind
81!
[1321]82! 1320 2014-03-20 08:40:49Z
[1320]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
[1008]89!
[1093]90! 1092 2013-02-02 11:24:22Z raasch
91! unused variables removed
92!
[1072]93! 1071 2012-11-29 16:54:55Z franke
94! Bugfix: collision efficiencies for Hall kernel should not be < 1.0E-20
95!
[1037]96! 1036 2012-10-22 13:43:42Z raasch
97! code put under GPL (PALM 3.9)
98!
[1008]99! 1007 2012-09-19 14:30:36Z franke
[1007]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
[829]104!
[850]105! 849 2012-03-15 10:35:09Z raasch
106! routine collision_efficiency_rogers added (moved from former advec_particles
107! to here)
108!
[836]109! 835 2012-02-22 11:21:19Z raasch $
110! Bugfix: array diss can be used only in case of Wang kernel
111!
[829]112! 828 2012-02-21 12:00:36Z raasch
[828]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
[790]117!
[828]118! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
119!
[826]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!
[791]124! 790 2011-11-29 03:11:20Z raasch
125! initial revision
[790]126!
127! Description:
128! ------------
[1682]129!> This module calculates collision efficiencies either due to pure gravitational
130!> effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or
[1822]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
[1682]134!> provided by L.-P. Wang but is substantially reformatted and speed optimized
135!> here.
[790]136!------------------------------------------------------------------------------!
[1682]137 MODULE lpm_collision_kernels_mod
138 
[790]139
[3274]140    USE basic_constants_and_equations_mod,                                     &
141        ONLY:  g, pi
[1320]142       
143    USE kinds
144
145    USE particle_attributes,                                                   &
[1822]146        ONLY:  collision_kernel, dissipation_classes, particles,               &
147               radius_classes
[1320]148
[828]149    USE pegrid
[790]150
[828]151
[790]152    IMPLICIT NONE
153
154    PRIVATE
155
[1822]156    PUBLIC  ckernel, init_kernels, rclass_lbound, rclass_ubound,               &
157            recalculate_kernel
[790]158
[1682]159    REAL(wp) ::  epsilon       !<
160    REAL(wp) ::  rclass_lbound !<
161    REAL(wp) ::  rclass_ubound !<
162    REAL(wp) ::  urms          !<
[790]163
[1822]164    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  epsclass  !< dissipation rate class
165    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  radclass  !< radius class
166    REAL(wp), DIMENSION(:),   ALLOCATABLE ::  winf      !<
[1320]167   
[1822]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   !<
[1320]173   
[1822]174    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ckernel !<
[792]175
[828]176    SAVE
[792]177
[790]178!
179!-- Public interfaces
[828]180    INTERFACE init_kernels
181       MODULE PROCEDURE init_kernels
182    END INTERFACE init_kernels
[790]183
[828]184    INTERFACE recalculate_kernel
185       MODULE PROCEDURE recalculate_kernel
186    END INTERFACE recalculate_kernel
[790]187
188
[828]189    CONTAINS
[790]190
[792]191
[828]192!------------------------------------------------------------------------------!
[1682]193! Description:
194! ------------
195!> Initialization of the collision efficiency matrix with fixed radius and
196!> dissipation classes, calculated at simulation start only.
[828]197!------------------------------------------------------------------------------!
[1682]198 
199    SUBROUTINE init_kernels
[792]200
[828]201       IMPLICIT NONE
[792]202
[1682]203       INTEGER(iwp) ::  i !<
204       INTEGER(iwp) ::  j !<
205       INTEGER(iwp) ::  k !<
[790]206
[828]207
208!
209!--    Calculate collision efficiencies for fixed radius- and dissipation
210!--    classes
211       IF ( collision_kernel(6:9) == 'fast' )  THEN
212
[1822]213          ALLOCATE( ckernel(1:radius_classes,1:radius_classes,                 &
214                    0:dissipation_classes), epsclass(1:dissipation_classes),   &
[828]215                    radclass(1:radius_classes) )
216
217!
218!--       Calculate the radius class bounds with logarithmic distances
[1858]219!--       in the interval [1.0E-6, 1000.0E-6] m
[1322]220          rclass_lbound = LOG( 1.0E-6_wp )
[1858]221          rclass_ubound = LOG( 1000.0E-6_wp )
[1822]222          radclass(1)   = EXP( rclass_lbound )
[828]223          DO  i = 2, radius_classes
224             radclass(i) = EXP( rclass_lbound +                                &
[1359]225                                ( rclass_ubound - rclass_lbound ) *            &
226                                ( i - 1.0_wp ) / ( radius_classes - 1.0_wp ) )
[828]227          ENDDO
228
229!
[1858]230!--       Set the class bounds for dissipation in interval [0.0, 600.0] cm**2/s**3
[828]231          DO  i = 1, dissipation_classes
[1858]232             epsclass(i) = 0.06_wp * REAL( i, KIND=wp ) / dissipation_classes
[828]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)
[1359]244             urms    = 2.02_wp * ( epsilon / 0.04_wp )**( 1.0_wp / 3.0_wp )
[828]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'
[1359]279             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6_wp, &
[1007]280                                              i = 1,radius_classes )
[828]281             DO  j = 1, radius_classes
[1007]282                WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j),  &
283                                          ( hkernel(i,j), i = 1,radius_classes )
[828]284             ENDDO
285
286             DO  k = 1, dissipation_classes
287                DO  i = 1, radius_classes
288                   DO  j = 1, radius_classes
[1359]289                      IF ( hkernel(i,j) == 0.0_wp )  THEN
290                         hwratio(i,j) = 9999999.9_wp
[828]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)
[1359]298                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i) * 1.0E6_wp, &
[1007]299                                                 i = 1,radius_classes )
[828]300                DO  j = 1, radius_classes
[1359]301                   WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j) * 1.0E6_wp, &
[1007]302                                          ( hwratio(i,j), i = 1,radius_classes )
[828]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
[790]315!------------------------------------------------------------------------------!
[1682]316! Description:
317! ------------
318!> Calculation of collision kernels during each timestep and for each grid box
[790]319!------------------------------------------------------------------------------!
[828]320    SUBROUTINE recalculate_kernel( i1, j1, k1 )
[790]321
[1320]322       USE arrays_3d,                                                          &
323           ONLY:  diss
[790]324
[1320]325       USE particle_attributes,                                                &
[1858]326           ONLY:  number_of_particles, prt_count, radius_classes, wang_kernel
[1320]327
[790]328       IMPLICIT NONE
329
[1682]330       INTEGER(iwp) ::  i      !<
331       INTEGER(iwp) ::  i1     !<
332       INTEGER(iwp) ::  j      !<
333       INTEGER(iwp) ::  j1     !<
334       INTEGER(iwp) ::  k1     !<
[790]335
336
[1858]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
[792]341
[1858]342       ALLOCATE( ec(1:number_of_particles,1:number_of_particles), &
343                 radclass(1:number_of_particles), winf(1:number_of_particles) )
[790]344
[828]345!
[1007]346!--    Store particle radii on the radclass array
[1858]347       radclass(1:number_of_particles) = particles(1:number_of_particles)%radius
[790]348
[835]349       IF ( wang_kernel )  THEN
[1007]350          epsilon = diss(k1,j1,i1)   ! dissipation rate in m**2/s**3
[835]351       ELSE
[1359]352          epsilon = 0.0_wp
[835]353       ENDIF
[1359]354       urms    = 2.02_wp * ( epsilon / 0.04_wp )**( 0.33333333333_wp )
[790]355
[1359]356       IF ( wang_kernel  .AND.  epsilon > 1.0E-7_wp )  THEN
[828]357!
358!--       Call routines to calculate efficiencies for the Wang kernel
[1858]359          ALLOCATE( gck(1:number_of_particles,1:number_of_particles), &
360                    ecf(1:number_of_particles,1:number_of_particles) )
[790]361
[828]362          CALL turbsd
363          CALL turb_enhance_eff
364          CALL effic
[790]365
[1858]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)
[790]369             ENDDO
[828]370          ENDDO
[790]371
[828]372          DEALLOCATE( gck, ecf )
[790]373
374       ELSE
[828]375!
376!--       Call routines to calculate efficiencies for the Hall kernel
[790]377          CALL fallg
378          CALL effic
379
[1858]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) )
[790]384             ENDDO
385          ENDDO
386
387       ENDIF
388
[828]389       DEALLOCATE( ec, radclass, winf )
[790]390
[828]391    END SUBROUTINE recalculate_kernel
[790]392
[828]393
[790]394!------------------------------------------------------------------------------!
[1682]395! Description:
396! ------------
[1822]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
[790]405!------------------------------------------------------------------------------!
[792]406    SUBROUTINE turbsd
[799]407
[1320]408       USE control_parameters,                                                 &
[3274]409           ONLY:  molecular_viscosity
[1320]410   
411       USE particle_attributes,                                                &
412           ONLY:  radius_classes
[790]413
414       IMPLICIT NONE
415
[1682]416       INTEGER(iwp) ::  i     !<
417       INTEGER(iwp) ::  j     !<
[790]418
[1682]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         !<
[790]462
[1822]463       REAL(wp), DIMENSION(1:radius_classes) ::  st  !< Stokes number
464       REAL(wp), DIMENSION(1:radius_classes) ::  tau !< inertial time scale
[790]465
[1822]466       lambda    = urms * SQRT( 15.0_wp * molecular_viscosity / epsilon )
[1322]467       lambda_re = urms**2 * SQRT( 15.0_wp / epsilon / molecular_viscosity )
[1822]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
[1007]472       vk        = eta / tauk
[790]473
[1359]474       ao = ( 11.0_wp + 7.0_wp * lambda_re ) / ( 205.0_wp + lambda_re )
[1822]475       tt = SQRT( 2.0_wp * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk
[799]476
[1822]477!
478!--    Get terminal velocity of droplets
479       CALL fallg
[790]480
[828]481       DO  i = 1, radius_classes
[1822]482          tau(i) = winf(i) / g    ! inertial time scale
483          st(i)  = tau(i) / tauk  ! Stokes number
[790]484       ENDDO
485
[828]486!
[1822]487!--    Calculate average radial relative velocity at contact (wrfin)
[828]488       z   = tt / tl
[1322]489       be  = SQRT( 2.0_wp ) * lambda / lf
[1359]490       bbb = SQRT( 1.0_wp - 2.0_wp * be**2 )
491       d1  = ( 1.0_wp + bbb ) / ( 2.0_wp * bbb )
[1822]492       e1  = lf * ( 1.0_wp + bbb ) * 0.5_wp
[1359]493       d2  = ( 1.0_wp - bbb ) * 0.5_wp / bbb
[1822]494       e2  = lf * ( 1.0_wp - bbb ) * 0.5_wp
[1359]495       ccc = SQRT( 1.0_wp - 2.0_wp * z**2 )
496       b1  = ( 1.0_wp + ccc ) * 0.5_wp / ccc
[1822]497       c1  = tl * ( 1.0_wp + ccc ) * 0.5_wp
[1359]498       b2  = ( 1.0_wp - ccc ) * 0.5_wp / ccc
[1822]499       c2  = tl * ( 1.0_wp - ccc ) * 0.5_wp
[790]500
[828]501       DO  i = 1, radius_classes
[790]502
[1822]503          v1 = winf(i)
504          t1 = tau(i)
[790]505
[828]506          DO  j = 1, i
[1007]507             rrp = radclass(i) + radclass(j)
[1822]508             v2  = winf(j)
509             t2  = tau(j)
[790]510
[1007]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)
[1822]513             v1xysq  = v1xysq * urms**2 / t1
514             vrms1xy = SQRT( v1xysq )
[790]515
[1007]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)
[1822]518             v2xysq  = v2xysq * urms**2 / t2
519             vrms2xy = SQRT( v2xysq )
[790]520
[828]521             IF ( winf(i) >= winf(j) )  THEN
[799]522                v1 = winf(i)
[790]523                t1 = tau(i)
[799]524                v2 = winf(j)
[790]525                t2 = tau(j)
526             ELSE
[799]527                v1 = winf(j)
[790]528                t1 = tau(j)
[799]529                v2 = winf(i)
[790]530                t2 = tau(i)
531             ENDIF
532
[828]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 )
[1822]538             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)
539             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0_wp * v1v2xy
[1359]540             IF ( wrtur2xy < 0.0_wp )  wrtur2xy = 0.0_wp
[1322]541             wrgrav2  = pi / 8.0_wp * ( winf(j) - winf(i) )**2
[1822]542             wrfin    = SQRT( ( 2.0_wp / pi ) * ( wrtur2xy + wrgrav2) )
[790]543
[828]544!
[1822]545!--          Calculate radial distribution function (grfin)
[828]546             IF ( st(j) > st(i) )  THEN
547                sst = st(j)
[790]548             ELSE
[828]549                sst = st(i)
[790]550             ENDIF
551
[1359]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 )
[790]556
[1007]557             c1_gr  =  xx / ( g / vk * tauk )**yy
[790]558
[1322]559             ao_gr  = ao + ( pi / 8.0_wp) * ( g / vk * tauk )**2
[1359]560             fao_gr = 20.115_wp * SQRT( ao_gr / lambda_re )
[1822]561             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta
[790]562
[1359]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
[790]565
[1822]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
[790]570             gck(j,i) = gck(i,j)
571
572          ENDDO
573       ENDDO
574
[828]575    END SUBROUTINE turbsd
[790]576
[1320]577    REAL(wp) FUNCTION phi_w( a, b, vsett, tau0 )
[1822]578!
579!--    Function used in the Ayala et al. (2008) analytical model for turbulent
580!--    effects on the collision kernel
[790]581       IMPLICIT NONE
582
[1682]583       REAL(wp) ::  a     !<
584       REAL(wp) ::  aa1   !<
585       REAL(wp) ::  b     !<
586       REAL(wp) ::  tau0  !<
587       REAL(wp) ::  vsett !<
[790]588
[1359]589       aa1 = 1.0_wp / tau0 + 1.0_wp / a + vsett / b
[1822]590       phi_w = 1.0_wp / aa1  - 0.5_wp * vsett / b / aa1**2
[790]591
[1007]592    END FUNCTION phi_w
[792]593
[1320]594    REAL(wp) FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
[1822]595!
596!--    Function used in the Ayala et al. (2008) analytical model for turbulent
597!--    effects on the collision kernel
[790]598       IMPLICIT NONE
599
[1682]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 !<
[790]612
[1359]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 ) -   &
[1822]623             vsett1 / aa2**2 + vsett2 / aa1**2 ) * 0.5_wp / b / aa3
[799]624
[828]625    END FUNCTION zhi
[790]626
[828]627
[790]628!------------------------------------------------------------------------------!
[1682]629! Description:
630! ------------
[1822]631!> Parameterization of terminal velocity following Rogers et al. (1993, J. Appl.
632!> Meteorol.)
[790]633!------------------------------------------------------------------------------!
[828]634    SUBROUTINE fallg
[790]635
[1320]636       USE particle_attributes,                                                &
637           ONLY:  radius_classes
[790]638
[828]639       IMPLICIT NONE
[790]640
[1822]641       INTEGER(iwp) ::  j                            !<
[790]642
[1822]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
[790]649
[1822]650       REAL(wp)            ::  diameter              !< droplet diameter in mm
[790]651
[799]652
[828]653       DO  j = 1, radius_classes
[790]654
[1822]655          diameter = radclass(j) * 2000.0_wp
[799]656
[1822]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 )
[828]662          ENDIF
[790]663
[828]664       ENDDO
[790]665
[828]666    END SUBROUTINE fallg
[790]667
[828]668
[790]669!------------------------------------------------------------------------------!
[1682]670! Description:
671! ------------
[1822]672!> Interpolation of collision efficiencies (Hall, 1980, J. Atmos. Sci.)
[790]673!------------------------------------------------------------------------------!
[828]674    SUBROUTINE effic
[1320]675 
676       USE particle_attributes,                                                &
677           ONLY:  radius_classes
[790]678
[828]679       IMPLICIT NONE
[790]680
[1682]681       INTEGER(iwp) ::  i  !<
682       INTEGER(iwp) ::  iq !<
683       INTEGER(iwp) ::  ir !<
684       INTEGER(iwp) ::  j  !<
685       INTEGER(iwp) ::  k  !<
[790]686
[1682]687       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ira !<
[790]688
[1682]689       LOGICAL, SAVE ::  first = .TRUE. !<
[790]690
[1682]691       REAL(wp) ::  ek              !<
692       REAL(wp) ::  particle_radius !<
693       REAL(wp) ::  pp              !<
694       REAL(wp) ::  qq              !<
695       REAL(wp) ::  rq              !<
[790]696
[1682]697       REAL(wp), DIMENSION(1:21), SAVE ::  rat        !<
[1320]698       
[1682]699       REAL(wp), DIMENSION(1:15), SAVE ::  r0         !<
[1320]700       
[1682]701       REAL(wp), DIMENSION(1:15,1:21), SAVE ::  ecoll !<
[790]702
[792]703!
[828]704!--    Initial assignment of constants
705       IF ( first )  THEN
[790]706
[792]707         first = .FALSE.
[1822]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,   &
[1359]710                  150.0_wp, 200.0_wp, 300.0_wp /)
[790]711
[1822]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,        &
[1359]715                  0.90_wp, 0.95_wp, 1.00_wp /)
716
[1822]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,    &
[1359]719                          0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp, 0.001_wp /)
[1822]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,    &
[1359]722                          0.200_wp, 0.500_wp, 0.770_wp, 0.870_wp, 0.970_wp /)
[1822]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,    &
[1359]725                          0.580_wp, 0.790_wp, 0.930_wp, 0.960_wp, 1.000_wp /)
[1822]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,    &
[1359]728                          0.750_wp, 0.910_wp, 0.970_wp, 0.980_wp, 1.000_wp /)
[1822]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,    &
[1359]731                          0.840_wp, 0.950_wp, 0.970_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]734                          0.880_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]737                          0.900_wp, 0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]740                          0.920_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]743                          0.940_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]746                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]749                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]752                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]755                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]758                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]761                          0.950_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]764                          0.970_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]767                          1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp, 1.000_wp /)
[1822]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,    &
[1359]770                          1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp, 1.020_wp /)
[1822]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,    &
[1359]773                          1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp, 1.040_wp /)
[1822]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,    &
[1359]776                          2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp, 2.300_wp /)
[1822]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,    &
[1359]779                          4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp, 4.000_wp /)
[828]780       ENDIF
[790]781
[792]782!
[828]783!--    Calculate the radius class index of particles with respect to array r
[1822]784!--    Radius has to be in microns
[828]785       ALLOCATE( ira(1:radius_classes) )
786       DO  j = 1, radius_classes
[1322]787          particle_radius = radclass(j) * 1.0E6_wp
[828]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
[790]796
[792]797!
[828]798!--    Two-dimensional linear interpolation of the collision efficiency.
[1822]799!--    Radius has to be in microns
[828]800       DO  j = 1, radius_classes
801          DO  i = 1, j
[792]802
[1880]803             ir = MAX( ira(i), ira(j) )
[1519]804             rq = MIN( radclass(i) / radclass(j), radclass(j) / radclass(i) )
[828]805             iq = INT( rq * 20 ) + 1
806             iq = MAX( iq , 2)
[792]807
[828]808             IF ( ir < 16 )  THEN
809                IF ( ir >= 2 )  THEN
[1822]810                   pp = ( ( MAX( radclass(j), radclass(i) ) * 1.0E6_wp ) -     &
811                          r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
[1359]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)           &
[828]817                             + pp * qq * ecoll(ir,iq)
818                ELSE
819                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
[1359]820                   ec(j,i) = ( 1.0_wp - qq ) * ecoll(1,iq-1) + qq * ecoll(1,iq)
[828]821                ENDIF
822             ELSE
823                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
[1359]824                ek = ( 1.0_wp - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
[1346]825                ec(j,i) = MIN( ek, 1.0_wp )
[1071]826             ENDIF
[792]827
[1359]828             IF ( ec(j,i) < 1.0E-20_wp )  ec(j,i) = 0.0_wp
[1071]829
[828]830             ec(i,j) = ec(j,i)
[792]831
[828]832          ENDDO
833       ENDDO
[792]834
[828]835       DEALLOCATE( ira )
[792]836
[828]837    END SUBROUTINE effic
[792]838
839
[790]840!------------------------------------------------------------------------------!
[1682]841! Description:
842! ------------
[1822]843!> Interpolation of turbulent enhancement factor for collision efficencies
844!> following Wang and Grabowski (2009, Atmos. Sci. Let.)
[790]845!------------------------------------------------------------------------------!
[828]846    SUBROUTINE turb_enhance_eff
[790]847
[1320]848       USE particle_attributes,                                                &
849           ONLY:  radius_classes
[790]850
[828]851       IMPLICIT NONE
[790]852
[1682]853       INTEGER(iwp) :: i  !<
854       INTEGER(iwp) :: iq !<
855       INTEGER(iwp) :: ir !<
856       INTEGER(iwp) :: j  !<
857       INTEGER(iwp) :: k  !<
858       INTEGER(iwp) :: kk !<
[790]859
[1682]860       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ira !<
[1320]861       
[1682]862       LOGICAL, SAVE ::  first = .TRUE. !<
[790]863
[1682]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              !<
[790]871
[1682]872       REAL(wp), DIMENSION(1:11), SAVE ::  rat           !<
873       REAL(wp), DIMENSION(1:7), SAVE  ::  r0            !<
[1320]874       
[1682]875       REAL(wp), DIMENSION(1:7,1:11), SAVE ::  ecoll_100 !<
876       REAL(wp), DIMENSION(1:7,1:11), SAVE ::  ecoll_400 !<
[799]877
878!
[828]879!--    Initial assignment of constants
880       IF ( first )  THEN
[799]881
[828]882          first = .FALSE.
[799]883
[1359]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 /)
[828]889!
[1822]890!--       Tabulated turbulent enhancement factor at 100 cm**2/s**3
[1359]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 /)
[1822]911          ecoll_100(:,11) = (/ 20.3_wp,   20.3_wp,   14.6_wp,   8.61_wp,  &
[1359]912                                2.60_wp,   2.60_wp,   1.0_wp /)
[828]913!
[1822]914!--       Tabulated turbulent enhancement factor at 400 cm**2/s**3
[1359]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 /)
[799]937
[828]938       ENDIF
[790]939
[828]940!
941!--    Calculate the radius class index of particles with respect to array r0
[1822]942!--    The droplet radius has to be given in microns.
[828]943       ALLOCATE( ira(1:radius_classes) )
[790]944
[828]945       DO  j = 1, radius_classes
[1322]946          particle_radius = radclass(j) * 1.0E6_wp
[828]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
[799]955
956!
[1822]957!--    Two-dimensional linear interpolation of the turbulent enhancement factor.
958!--    The droplet radius has to be given in microns.
[828]959       DO  j =  1, radius_classes
960          DO  i = 1, j
[799]961
[1880]962             ir = MAX( ira(i), ira(j) )
[1519]963             rq = MIN( radclass(i) / radclass(j), radclass(j) / radclass(i) )
[799]964
[828]965             DO  kk = 2, 11
966                IF ( rq <= rat(kk) )  THEN
967                   iq = kk
968                   EXIT
969                ENDIF
970             ENDDO
[790]971
[1822]972             y1 = 1.0_wp  ! turbulent enhancement factor at 0 m**2/s**3
[1007]973
[828]974             IF ( ir < 8 )  THEN
975                IF ( ir >= 2 )  THEN
[1822]976                   pp = ( MAX( radclass(j), radclass(i) ) * 1.0E6_wp -  &
977                          r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
[828]978                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
[1359]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)
[828]987                ELSE
988                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
[1359]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)
[828]991                ENDIF
992             ELSE
993                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
[1359]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)
[828]996             ENDIF
997!
[1822]998!--          Linear interpolation of turbulent enhancement factor
[1359]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
[828]1005             ELSE
[1359]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
[828]1008             ENDIF
[790]1009
[1359]1010             IF ( ecf(j,i) < 1.0_wp )  ecf(j,i) = 1.0_wp
[790]1011
[828]1012             ecf(i,j) = ecf(j,i)
[790]1013
[828]1014          ENDDO
1015       ENDDO
[790]1016
[828]1017    END SUBROUTINE turb_enhance_eff
[790]1018
[825]1019 END MODULE lpm_collision_kernels_mod
Note: See TracBrowser for help on using the repository browser.