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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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