Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Feb 21, 2012 12:00:36 PM (11 years ago)
Author:
raasch
Message:

New:
---

Changed:

Optimization of collision kernels. Collision tables can be calculated once at
simulation start for defined radius (and dissipation) classes instead of
re-calculating them at every timestep and for the particle ensemble in
every gridbox.
For this purpose the particle feature color is renamed class.
New parpar parameters radius_classes and dissipation_classes.
(Makefile, advec_particles, check_parameters, data_output_dvrp, header, init_particles, lpm_collision_kernels, modules, package_parin, set_particle_attributes)

Lower limit for droplet radius changed from 1E-7 to 1E-8.
(advec_particles)

Complete re-formatting of collision code (including changes in names of
variables, modules and subroutines).
(advec_particles, lpm_collision_kernels)

Errors:

Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
(advec_particles, lpm_collision_kernels)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/lpm_collision_kernels.f90

    r826 r828  
    1 MODULE lpm_collision_kernels_mod
     1 MODULE lpm_collision_kernels_mod
    22
    33!------------------------------------------------------------------------------!
    44! Current revisions:
    55! -----------------
    6 !
     6! code has been completely reformatted, routine colker renamed
     7! recalculate_kernel,
     8! routine init_kernels added, radius is now communicated to the collision
     9! routines by array radclass
     10!
     11! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
    712!
    813! Former revisions:
    … …  
    2429! Description:
    2530! ------------
    26 ! This routine calculates the effect of (SGS) turbulence on the collision
    27 ! efficiency of droplets.
    28 ! It is based on the original kernel developed by Wang (...)
     31! This module calculates collision efficiencies either due to pure gravitational
     32! effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or
     33! including the effects of (SGS) turbulence (Wang kernel, see Wang and
     34! Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8). The original code has been
     35! provided by L.-P. Wang but is substantially reformatted and speed optimized
     36! here.
     37!
     38! ATTENTION:
     39! Physical quantities (like g, densities, etc.) used in this module still
     40! have to be adjusted to those values used in the main PALM code.
     41! Also, quantities in CGS-units should be converted to SI-units eventually.
    2942!------------------------------------------------------------------------------!
    3043
    … …  
    3346    USE constants
    3447    USE particle_attributes
     48    USE pegrid
     49
    3550
    3651    IMPLICIT NONE
    … …  
    3853    PRIVATE
    3954
    40     PUBLIC  colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi
    41 
    42     INTEGER, SAVE ::  ip, jp, kp, pend, pstart, psum
    43 
    44     REAL, SAVE :: epsilon, eps2, urms, urms2
    45 
    46     REAL, DIMENSION(:),   ALLOCATABLE, SAVE ::  winf
    47     REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  ec, ecf, gck
     55    PUBLIC  ckernel, init_kernels,  rclass_lbound, rclass_ubound, &
     56            recalculate_kernel           
     57
     58    REAL ::  epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2
     59
     60    REAL, DIMENSION(:),   ALLOCATABLE ::  epsclass, radclass, winf
     61    REAL, DIMENSION(:,:), ALLOCATABLE ::  ec, ecf, gck, hkernel, hwratio
     62    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  ckernel
     63
     64    SAVE
    4865
    4966!
    5067!-- Public interfaces
    51     INTERFACE colker
    52        MODULE PROCEDURE colker
    53     END INTERFACE colker
    54 
    55     INTERFACE effic
    56        MODULE PROCEDURE effic
    57     END INTERFACE effic
    58 
    59     INTERFACE fallg
    60        MODULE PROCEDURE fallg
    61     END INTERFACE fallg
    62 
    63     INTERFACE phi
    64        MODULE PROCEDURE phi
    65     END INTERFACE phi
    66 
    67     INTERFACE turbsd
    68        MODULE PROCEDURE turbsd
    69     END INTERFACE turbsd
    70 
    71     INTERFACE turb_enhance_eff
    72        MODULE PROCEDURE turb_enhance_eff
    73     END INTERFACE turb_enhance_eff
    74 
    75     INTERFACE zhi
    76        MODULE PROCEDURE zhi
    77     END INTERFACE zhi
     68    INTERFACE init_kernels
     69       MODULE PROCEDURE init_kernels
     70    END INTERFACE init_kernels
     71
     72    INTERFACE recalculate_kernel
     73       MODULE PROCEDURE recalculate_kernel
     74    END INTERFACE recalculate_kernel
     75
    7876
    7977    CONTAINS
    8078
    81 !------------------------------------------------------------------------------!
    82 ! SUBROUTINE for calculation of collision kernel
    83 !------------------------------------------------------------------------------!
    84     SUBROUTINE colker( i1, j1, k1, kernel )
     79
     80    SUBROUTINE init_kernels
     81!------------------------------------------------------------------------------!
     82! Initialization of the collision efficiency matrix with fixed radius and
     83! dissipation classes, calculated at simulation start only.
     84!------------------------------------------------------------------------------!
     85
     86       IMPLICIT NONE
     87
     88       INTEGER ::  i, j, k
     89
     90
     91!
     92!--    Calculate collision efficiencies for fixed radius- and dissipation
     93!--    classes
     94       IF ( collision_kernel(6:9) == 'fast' )  THEN
     95
     96          ALLOCATE( ckernel(1:radius_classes,1:radius_classes,               &
     97                    0:dissipation_classes), epsclass(1:dissipation_classes), &
     98                    radclass(1:radius_classes) )
     99
     100!
     101!--       Calculate the radius class bounds with logarithmic distances
     102!--       in the interval [1.0E-6, 2.0E-4] m
     103          rclass_lbound = LOG( 1.0E-6 )
     104          rclass_ubound = LOG( 2.0E-4 )
     105          radclass(1)   = 1.0E-6
     106          DO  i = 2, radius_classes
     107             radclass(i) = EXP( rclass_lbound +                                &
     108                                ( rclass_ubound - rclass_lbound ) * ( i-1.0 ) /&
     109                                ( radius_classes - 1.0 ) )
     110!             IF ( myid == 0 )  THEN
     111!                PRINT*, 'i=', i, ' r = ', radclass(i)*1.0E6
     112!             ENDIF
     113          ENDDO
     114!
     115!--       Collision routines expect radius to be in cm
     116          radclass = radclass * 100.0
     117
     118!
     119!--       Set the class bounds for dissipation in interval [0, 1000] cm**2/s**3
     120          DO  i = 1, dissipation_classes
     121             epsclass(i) = 1000.0 * REAL( i ) / dissipation_classes
     122!             IF ( myid == 0 )  THEN
     123!                PRINT*, 'i=', i, ' eps = ', epsclass(i)
     124!             ENDIF
     125          ENDDO
     126!
     127!--       Calculate collision efficiencies of the Wang/ayala kernel
     128          ALLOCATE( ec(1:radius_classes,1:radius_classes),  &
     129                    ecf(1:radius_classes,1:radius_classes), &
     130                    gck(1:radius_classes,1:radius_classes), &
     131                    winf(1:radius_classes) )
     132
     133          DO  k = 1, dissipation_classes
     134
     135             epsilon = epsclass(k)
     136             urms    = 202.0 * ( epsilon / 400.0 )**( 1.0 / 3.0 )
     137
     138             CALL turbsd
     139             CALL turb_enhance_eff
     140             CALL effic
     141
     142             DO  j = 1, radius_classes
     143                DO  i = 1, radius_classes
     144                   ckernel(i,j,k) = ec(i,j) * gck(i,j) * ecf(i,j)
     145                ENDDO
     146             ENDDO
     147
     148          ENDDO
     149
     150!
     151!--       Calculate collision efficiencies of the Hall kernel
     152          ALLOCATE( hkernel(1:radius_classes,1:radius_classes), &
     153                    hwratio(1:radius_classes,1:radius_classes) )
     154
     155          CALL fallg
     156          CALL effic
     157
     158          DO  j = 1, radius_classes
     159             DO  i =  1, radius_classes
     160                hkernel(i,j) = pi * ( radclass(j) + radclass(i) )**2 &
     161                                  * ec(i,j) * ABS( winf(j) - winf(i) )
     162                ckernel(i,j,0) = hkernel(i,j)  ! hall kernel stored on index 0
     163              ENDDO
     164          ENDDO
     165
     166!
     167!--       Test output of efficiencies
     168          IF ( j == -1 )  THEN
     169
     170             PRINT*, '*** Hall kernel'
     171             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
     172             DO  j = 1, radius_classes
     173                WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes )
     174             ENDDO
     175
     176             DO  k = 1, dissipation_classes
     177                DO  i = 1, radius_classes
     178                   DO  j = 1, radius_classes
     179                      IF ( hkernel(i,j) == 0.0 )  THEN
     180                         hwratio(i,j) = 9999999.9
     181                      ELSE
     182                         hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j)
     183                      ENDIF
     184                   ENDDO
     185                ENDDO
     186
     187                PRINT*, '*** epsilon = ', epsclass(k)
     188                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
     189                DO  j = 1, radius_classes
     190!                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( ckernel(i,j,k), i = 1,radius_classes )
     191                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( hwratio(i,j), i = 1,radius_classes )
     192                ENDDO
     193             ENDDO
     194
     195          ENDIF
     196
     197          DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf )
     198
     199          ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
     200
     201       ELSEIF( collision_kernel == 'hall'  .OR.  collision_kernel == 'wang' ) &
     202       THEN
     203!
     204!--       Initial settings for Hall- and Wang-Kernel
     205!--       To be done: move here parts from turbsd, fallg, ecoll, etc.
     206       ENDIF
     207
     208    END SUBROUTINE init_kernels
     209
     210
     211!------------------------------------------------------------------------------!
     212! Calculation of collision kernels during each timestep and for each grid box
     213!------------------------------------------------------------------------------!
     214    SUBROUTINE recalculate_kernel( i1, j1, k1 )
    85215
    86216       USE arrays_3d
    … …  
    94224       IMPLICIT NONE
    95225
    96        INTEGER ::  i, i1, j, j1, k1
    97 
    98        REAL, DIMENSION(prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1,prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1) :: kernel
    99 
    100 !       CALL cpu_log( log_point_s(46), 'colker', 'start' )
    101 
    102        ip = i1
    103        jp = j1
    104        kp = k1
    105 
    106        pstart = prt_start_index(kp,jp,ip)
    107        pend   = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1
    108        psum   = prt_count(kp,jp,ip)
    109 
    110        ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
    111 
    112        IF ( wang_kernel )  THEN
    113 
    114           ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
    115 
    116           epsilon = diss(kp,jp,ip) * 1.0E5  !dissipation rate in cm**2/s**-3
    117           urms     = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 )
    118 
    119           IF ( epsilon <= 0.001 )  THEN
    120 
    121              CALL fallg
    122              CALL effic
    123 
    124              DO  j = pstart, pend
    125                 DO  i =  pstart, pend
    126                    kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
    127                                         particles(i)%radius * 100.0 )**2 &
    128                                     * ec(i,j) * ABS( winf(j) - winf(i) )
    129                 ENDDO
     226       INTEGER ::  i, i1, j, j1, k1, pend, pstart
     227
     228
     229       pstart = prt_start_index(k1,j1,i1)
     230       pend   = prt_start_index(k1,j1,i1) + prt_count(k1,j1,i1) - 1
     231       radius_classes = prt_count(k1,j1,i1)
     232
     233       ALLOCATE( ec(1:radius_classes,1:radius_classes), &
     234                 radclass(1:radius_classes), winf(1:radius_classes) )
     235
     236!
     237!--    Store particle radii on the radclass array. Collision routines
     238!--    expect radii to be in cm.
     239       radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0
     240
     241       epsilon = diss(k1,j1,i1) * 1.0E4   ! dissipation rate in cm**2/s**-3
     242       urms    = 202.0 * ( epsilon / 400.0 )**( 0.33333333333 )
     243
     244       IF ( wang_kernel  .AND.  epsilon > 0.001 )  THEN
     245!
     246!--       Call routines to calculate efficiencies for the Wang kernel
     247          ALLOCATE( gck(1:radius_classes,1:radius_classes), &
     248                    ecf(1:radius_classes,1:radius_classes) )
     249
     250          CALL turbsd
     251          CALL turb_enhance_eff
     252          CALL effic
     253
     254          DO  j = 1, radius_classes
     255             DO  i =  1, radius_classes
     256                ckernel(pstart+i-1,pstart+j-1,1) = ec(i,j) * gck(i,j) * ecf(i,j)
    130257             ENDDO
    131 
    132           ELSE
    133 
    134              CALL turbsd
    135              CALL turb_enhance_eff
    136              CALL effic
    137 
    138              DO  j = pstart, pend, 1
    139                 DO  i =  pstart, pend, 1
    140                    kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
    141                 ENDDO
     258          ENDDO
     259
     260          DEALLOCATE( gck, ecf )
     261
     262       ELSE
     263!
     264!--       Call routines to calculate efficiencies for the Hall kernel
     265          CALL fallg
     266          CALL effic
     267
     268          DO  j = 1, radius_classes
     269             DO  i =  1, radius_classes
     270                ckernel(pstart+i-1,pstart+j-1,1) = pi *                       &
     271                                          ( radclass(j) + radclass(i) )**2    &
     272                                          * ec(i,j) * ABS( winf(j) - winf(i) )
    142273             ENDDO
    143 
    144           ENDIF
    145 
    146           DEALLOCATE(gck, ecf)
    147 
    148        ELSE
    149 
    150 !          CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )
    151           CALL fallg
    152 !          CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )
    153 !          CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )
    154           CALL effic
    155 !          CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )
    156 
    157           DO  j = pstart, pend
    158              DO  i =  pstart, pend
    159                 kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
    160                                      particles(i)%radius * 100.0 )**2 &
    161                                  * ec(i,j) * ABS( winf(j) - winf(i) )
    162              ENDDO
    163274          ENDDO
    164275
    165276       ENDIF
    166277
    167        DEALLOCATE( ec, winf )
    168 
    169 !       CALL cpu_log( log_point_s(46), 'colker', 'stop' )
    170 
    171     END SUBROUTINE colker
    172 
    173 !------------------------------------------------------------------------------!
    174 ! SUBROUTINE for calculation of w, g and gck
     278       ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
     279
     280       DEALLOCATE( ec, radclass, winf )
     281
     282    END SUBROUTINE recalculate_kernel
     283
     284
     285!------------------------------------------------------------------------------!
     286! Calculation of gck
     287! This is from Aayala 2008b, page 37ff.
     288! Necessary input parameters: water density, radii of droplets, air density,
     289! air viscosity, turbulent dissipation rate, taylor microscale reynolds number,
     290! gravitational acceleration  --> to be replaced by PALM parameters
    175291!------------------------------------------------------------------------------!
    176292    SUBROUTINE turbsd
    177 ! from Aayala 2008b, page 37ff, necessary input parameter water density, radii
    178 ! of droplets, air density, air viscosity, turbulent dissipation rate,
    179 ! taylor microscale reynolds number, gravitational acceleration
    180293
    181294       USE constants
    … …  
    186299       IMPLICIT NONE
    187300
    188        REAL :: Relamda, &
    189                Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
    190                bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
    191                vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, &
    192                wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
    193 
    194        REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens
    195 
    196        REAL, DIMENSION(pstart:pend) :: St, tau
     301       INTEGER ::  i, j
    197302
    198303       LOGICAL, SAVE ::  first = .TRUE.
    199304
    200        INTEGER :: i, j
    201 
    202 !
    203 !--   Initial assignment of constants
     305       REAL ::  ao, ao_gr, bbb, be, b1, b2, ccc, c1, c1_gr, c2, d1, d2, eta, &
     306                e1, e2, fao_gr, fr, grfin, lambda, lambda_re, lf, rc, rrp,   &
     307                sst, tauk, tl, t2, tt, t1, vk, vrms1xy, vrms2xy, v1, v1v2xy, &
     308                v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z
     309
     310       REAL, SAVE ::  airdens, airvisc, anu, gravity, waterdens
     311
     312       REAL, DIMENSION(1:radius_classes) ::  st, tau
     313
     314
     315!
     316!--    Initial assignment of constants
    204317       IF ( first )  THEN
    205318
    206319          first = .FALSE.
    207           airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
    208           airdens   = 1.2250   !air density in mg/cm**3
    209           waterdens = 1000.0   !water density in mg/cm**3
    210           gravity   = 980.6650 !in cm/s**2
     320          airvisc   = 0.1818     ! dynamic viscosity in mg/cm*s
     321          airdens   = 1.2250     ! air density in mg/cm**3
     322          waterdens = 1000.0     ! water density in mg/cm**3
     323          gravity   = 980.6650   ! in cm/s**2
    211324          anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
    212325
    213326       ENDIF   
    214327
    215        Relamda = urms**2*sqrt(15.0/epsilon/anu)
    216 
    217        Tl = urms**2/epsilon       !in s
    218        Lf = 0.5 * (urms**3)/epsilon !in cm
    219 
    220        tauk = (anu/epsilon)**0.5       !in s
    221        eta  = (anu**3/epsilon)**0.25  !in cm
    222        vk   = eta/tauk                   !in cm/s
    223 
    224        ao = (11.+7.*Relamda)/(205.+Relamda)
    225 
    226        lambda = urms * sqrt(15.*anu/epsilon)     !in cm
    227 
    228        tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s
    229 
    230        CALL fallg !gives winf in cm/s
    231 
    232        DO i = pstart, pend
    233           tau(i) = winf(i)/gravity    !in s
    234           St(i)  = tau(i)/tauk
     328       lambda    = urms * SQRT( 15.0 * anu / epsilon )     ! in cm
     329       lambda_re = urms**2 * SQRT( 15.0 / epsilon / anu )
     330       tl        = urms**2 / epsilon                       ! in s
     331       lf        = 0.5 * urms**3 / epsilon                 ! in cm
     332       tauk      = SQRT( anu / epsilon )                   ! in s
     333       eta       = ( anu**3 / epsilon )**0.25              ! in cm
     334       vk        = eta / tauk                              ! in cm/s
     335
     336       ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
     337       tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk   ! in s
     338
     339       CALL fallg    ! gives winf in cm/s
     340
     341       DO  i = 1, radius_classes
     342          tau(i) = winf(i) / gravity    ! in s
     343          st(i)  = tau(i) / tauk
    235344       ENDDO
    236345
    237 !***** TO CALCULATE wr ********************************
    238 !from Aayala 2008b, page 38f
    239 
    240        z  = tt/Tl
    241        be = sqrt(2.0)*lambda/Lf
    242 
    243        bbb = sqrt(1.0-2.0*be**2)
    244        d1 = (1.+bbb)/2.0/bbb
    245        e1 = Lf*(1.0+bbb)/2.0   !in cm
    246        d2 = (1.0-bbb)/2.0/bbb
    247        e2 = Lf*(1.0-bbb)/2.0   !in cm
    248 
    249        ccc = sqrt(1.0-2.0*z**2)
    250        b1 = (1.+ccc)/2./ccc
    251        c1 = Tl*(1.+ccc)/2.   !in s
    252        b2 = (1.-ccc)/2./ccc
    253        c2 = Tl*(1.-ccc)/2.   !in s
    254 
    255        DO i = pstart, pend
    256           v1 = winf(i)         !in cm/s
    257           t1 = tau(i)         !in s
    258 
    259           DO j = pstart,i
    260              rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm
    261              v2  = winf(j)         !in cm/s
    262              t2  = tau(j)         !in s
    263 
    264              v1xysq =  b1*d1*PHI(c1,e1,v1,t1) &
    265                      - b1*d2*PHI(c1,e2,v1,t1) &
    266                      - b2*d1*PHI(c2,e1,v1,t1) &
    267                      + b2*d2*PHI(c2,e2,v1,t1)
    268              v1xysq = v1xysq * urms**2/t1    !in cm**2/s**2
    269 
    270              vrms1xy= sqrt(v1xysq)             !in cm/s
    271 
    272              v2xysq =  b1*d1*PHI(c1,e1,v2,t2) &
    273                      - b1*d2*PHI(c1,e2,v2,t2) &
    274                      - b2*d1*PHI(c2,e1,v2,t2) &
    275                      + b2*d2*PHI(c2,e2,v2,t2)
    276              v2xysq = v2xysq * urms**2/t2    !in cm**2/s**2
    277 
    278               vrms2xy= sqrt(v2xysq)             !in cm/s
    279 
    280              IF(winf(i).ge.winf(j)) THEN
     346!
     347!--    Calculate wr (from Aayala 2008b, page 38f)
     348       z   = tt / tl
     349       be  = SQRT( 2.0 ) * lambda / lf
     350       bbb = SQRT( 1.0 - 2.0 * be**2 )
     351       d1  = ( 1.0 + bbb ) / ( 2.0 * bbb )
     352       e1  = lf * ( 1.0 + bbb ) * 0.5   ! in cm
     353       d2  = ( 1.0 - bbb ) * 0.5 / bbb
     354       e2  = lf * ( 1.0 - bbb ) * 0.5   ! in cm
     355       ccc = SQRT( 1.0 - 2.0 * z**2 )
     356       b1  = ( 1.0 + ccc ) * 0.5 / ccc
     357       c1  = tl * ( 1.0 + ccc ) * 0.5   ! in s
     358       b2  = ( 1.0 - ccc ) * 0.5 / ccc
     359       c2  = tl * ( 1.0 - ccc ) * 0.5   ! in s
     360
     361       DO  i = 1, radius_classes
     362
     363          v1 = winf(i)        ! in cm/s
     364          t1 = tau(i)         ! in s
     365
     366          DO  j = 1, i
     367             rrp = radclass(i) + radclass(j)               ! radius in cm
     368             v2  = winf(j)                                 ! in cm/s
     369             t2  = tau(j)                                  ! in s
     370
     371             v1xysq  = b1 * d1 * phi(c1,e1,v1,t1) - b1 * d2 * phi(c1,e2,v1,t1) &
     372                     - b2 * d1 * phi(c2,e1,v1,t1) + b2 * d2 * phi(c2,e2,v1,t1)
     373             v1xysq  = v1xysq * urms**2 / t1                ! in cm**2/s**2
     374             vrms1xy = SQRT( v1xysq )                       ! in cm/s
     375
     376             v2xysq  = b1 * d1 * phi(c1,e1,v2,t2) - b1 * d2 * phi(c1,e2,v2,t2) &
     377                     - b2 * d1 * phi(c2,e1,v2,t2) + b2 * d2 * phi(c2,e2,v2,t2)
     378             v2xysq  = v2xysq * urms**2 / t2                ! in cm**2/s**2
     379             vrms2xy = SQRT( v2xysq )                       ! in cm/s
     380
     381             IF ( winf(i) >= winf(j) )  THEN
    281382                v1 = winf(i)
    282383                t1 = tau(i)
    … …  
    290391             ENDIF
    291392
    292              v1v2xy =  b1*d1*ZHI(c1,e1,v1,t1,v2,t2) &
    293                      - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) &
    294                      - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) &
    295                      + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
    296              fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
    297              v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j)    !in cm**2/s**2
    298 
    299              wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy  !in cm**2/s**2
    300 
    301               IF (wrtur2xy.le.0.0) wrtur2xy=0.0
    302 
    303               wrgrav2=pi/8.*(winf(j)-winf(i))**2
    304 
    305               wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
    306 
    307 
    308 !***** TO CALCULATE gr ********************************
    309 
    310              IF(St(j).gt.St(i)) THEN
    311                 SSt = St(j)
     393             v1v2xy   =  b1 * d1 * zhi(c1,e1,v1,t1,v2,t2) - &
     394                         b1 * d2 * zhi(c1,e2,v1,t1,v2,t2) - &
     395                         b2 * d1 * zhi(c2,e1,v1,t1,v2,t2) + &
     396                         b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
     397             fr       = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
     398             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)  ! in cm**2/s**2
     399             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in cm**2/s**2
     400             IF ( wrtur2xy < 0.0 )  wrtur2xy = 0.0
     401             wrgrav2  = pi / 8.0 * ( winf(j) - winf(i) )**2
     402             wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in cm/s
     403
     404!
     405!--          Calculate gr
     406             IF ( st(j) > st(i) )  THEN
     407                sst = st(j)
    312408             ELSE
    313                 SSt = St(i)
     409                sst = st(i)
    314410             ENDIF
    315411
    316              XX = -0.1988*SSt**4 + 1.5275*SSt**3 &
    317                   -4.2942*SSt**2 + 5.3406*SSt
    318 
    319              IF(XX.le.0.0) XX = 0.0
    320 
    321              YY = 0.1886*exp(20.306/Relamda)
    322 
    323              c1_gr = XX/(gravity/(vk/tauk))**YY
    324 
    325              ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2
    326              fao_gr = 20.115 * (ao_gr/Relamda)**0.5
    327              rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
    328 
    329              grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.)
    330              IF (grFIN.lt.1.0) grFIN = 1.0
    331 
    332              gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN       ! in cm**3/s
     412             xx = -0.1988 * sst**4 + 1.5275 * sst**3 - 4.2942 * sst**2 + &
     413                   5.3406 * sst
     414             IF ( xx < 0.0 )  xx = 0.0
     415             yy = 0.1886 * EXP( 20.306 / lambda_re )
     416
     417             c1_gr  =  xx / ( gravity / vk * tauk )**yy
     418
     419             ao_gr  = ao + ( pi / 8.0) * ( gravity / vk * tauk )**2
     420             fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
     421             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
     422
     423             grfin  = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5 )
     424             IF ( grfin < 1.0 )  grfin = 1.0
     425
     426             gck(i,j) = 2.0 * pi * rrp**2 * wrfin * grfin           ! in cm**3/s
    333427             gck(j,i) = gck(i,j)
    334428
    … …  
    336430       ENDDO
    337431
    338     END SUBROUTINE TurbSD
    339 
    340 !------------------------------------------------------------------------------!
    341 ! PHI as a function
    342 !------------------------------------------------------------------------------!
    343   REAL FUNCTION PHI(a,b,vsett,tau0)
     432    END SUBROUTINE turbsd
     433
     434
     435!------------------------------------------------------------------------------!
     436! phi as a function
     437!------------------------------------------------------------------------------!
     438    REAL FUNCTION phi( a, b, vsett, tau0 )
    344439
    345440       IMPLICIT NONE
    346441
    347        REAL :: a, aa1, b, vsett, tau0
    348 
    349        aa1 = 1./tau0 + 1./a + vsett/b
    350 
    351        PHI = 1./aa1 - vsett/2.0/b/aa1**2  !in s
    352 
    353     END FUNCTION PHI
    354 
    355 !------------------------------------------------------------------------------!
    356 ! ZETA as a function
    357 !------------------------------------------------------------------------------!
    358   REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)
     442       REAL ::  a, aa1, b, tau0, vsett
     443
     444       aa1 = 1.0 / tau0 + 1.0 / a + vsett / b
     445       phi = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
     446
     447    END FUNCTION phi
     448
     449
     450!------------------------------------------------------------------------------!
     451! zeta as a function
     452!------------------------------------------------------------------------------!
     453    REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
    359454
    360455       IMPLICIT NONE
    361456
    362        REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2
    363 
    364         aa1 = vsett2/b - 1./tau2 - 1./a
    365         aa2 = vsett1/b + 1./tau1 + 1./a
    366         aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
    367         aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2
    368         aa5 = vsett2/b + 1./tau2 + 1./a
    369         aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2
    370         ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2          &
    371              + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6      &
    372              + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2)   &
    373                                                            * 1./2./b/aa3             ! in s**2
    374 
    375     END FUNCTION ZHI
    376 
    377 !------------------------------------------------------------------------------!
    378 ! SUBROUTINE for calculation of terminal velocity winf
    379 !------------------------------------------------------------------------------!
    380    SUBROUTINE fallg
     457       REAL ::  a, aa1, aa2, aa3, aa4, aa5, aa6, b, tau1, tau2, vsett1, vsett2
     458
     459       aa1 = vsett2 / b - 1.0 / tau2 - 1.0 / a
     460       aa2 = vsett1 / b + 1.0 / tau1 + 1.0 / a
     461       aa3 = ( vsett1 - vsett2 ) / b + 1.0 / tau1 + 1.0 / tau2
     462       aa4 = ( vsett2 / b )**2 - ( 1.0 / tau2 + 1.0 / a )**2
     463       aa5 = vsett2 / b + 1.0 / tau2 + 1.0 / a
     464       aa6 = 1.0 / tau1 - 1.0 / a + ( 1.0 / tau2 + 1.0 / a) * vsett1 / vsett2
     465       zhi = (1.0 / aa1 - 1.0 / aa2 ) * ( vsett1 - vsett2 ) * 0.5 / b / aa3**2 &
     466           + (4.0 / aa4 - 1.0 / aa5**2 - 1.0 / aa1**2 ) * vsett2 * 0.5 / b /aa6&
     467           + (2.0 * ( b / aa2 - b / aa1 ) - vsett1 / aa2**2 + vsett2 / aa1**2 )&
     468           * 0.5 / b / aa3      ! in s**2
     469
     470    END FUNCTION zhi
     471
     472
     473!------------------------------------------------------------------------------!
     474! Calculation of terminal velocity winf
     475!------------------------------------------------------------------------------!
     476    SUBROUTINE fallg
    381477
    382478       USE constants
    … …  
    385481       USE arrays_3d
    386482
    387       IMPLICIT NONE
    388 
    389       INTEGER :: i, j
    390 
    391       LOGICAL, SAVE ::  first = .TRUE.
    392 
    393       REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py
    394 
    395       REAL :: bond, x, xrey, y
    396 
    397       REAL, DIMENSION(1:7), SAVE  :: b
    398       REAL, DIMENSION(1:6), SAVE  :: c
    399 
    400 !
    401 !--   Initial assignment of constants
    402       IF ( first )  THEN
     483       IMPLICIT NONE
     484
     485       INTEGER ::  i, j
     486
     487       LOGICAL, SAVE ::  first = .TRUE.
     488
     489       REAL, SAVE ::  cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, &
     490                      t0, xlamb
     491
     492       REAL ::  bond, x, xrey, y
     493
     494       REAL, DIMENSION(1:7), SAVE  ::  b
     495       REAL, DIMENSION(1:6), SAVE  ::  c
     496
     497!
     498!--    Initial assignment of constants
     499       IF ( first )  THEN
     500
     501          first = .FALSE.
     502          b = (/  -0.318657E1,  0.992696E0, -0.153193E-2, -0.987059E-3, &
     503                 -0.578878E-3, 0.855176E-4, -0.327815E-5 /)
     504          c = (/  -0.500015E1,  0.523778E1,  -0.204914E1,   0.475294E0, &
     505                 -0.542819E-1, 0.238449E-2 /)
     506
     507          eta   = 1.818E-4          ! in poise = g/(cm s)
     508          xlamb = 6.62E-6           ! in cm
     509          rhow  = 1.0               ! in g/cm**3
     510          rhoa  = 1.225E-3          ! in g/cm**3
     511          grav  = 980.665           ! in cm/s**2
     512          cunh  = 1.257 * xlamb     ! in cm
     513          t0    = 273.15            ! in K
     514          sigma = 76.1 - 0.155 * ( 293.15 - t0 )              ! in N/m = g/s**2
     515          stok  = 2.0  * grav * ( rhow - rhoa ) / ( 9.0 * eta ) ! in 1/(cm s)
     516          stb   = 32.0 * rhoa * ( rhow - rhoa) * grav / (3.0 * eta * eta)
     517                                                                ! in 1/cm**3
     518          phy   = sigma**3 * rhoa**2 / ( eta**4 * grav * ( rhow - rhoa ) )
     519          py    = phy**( 1.0 / 6.0 )
     520
     521       ENDIF
     522
     523       DO  j = 1, radius_classes
     524
     525          IF ( radclass(j) <= 1.0E-3 ) THEN
     526
     527             winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) ! in cm/s
     528
     529          ELSEIF ( radclass(j) > 1.0E-3  .AND.  radclass(j) <= 5.35E-2 )  THEN
     530
     531             x = LOG( stb * radclass(j)**3 )
     532             y = 0.0
     533
     534             DO  i = 1, 7
     535                y = y + b(i) * x**(i-1)
     536             ENDDO
     537
     538             xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
     539             winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )      ! in cm/s
     540
     541          ELSEIF ( radclass(j) > 5.35E-2 )  THEN
     542
     543             IF ( radclass(j) > 0.35 )  THEN
     544                bond = grav * ( rhow - rhoa ) * 0.35**2 / sigma
     545             ELSE
     546                bond = grav * ( rhow - rhoa ) * radclass(j)**2 / sigma
     547             ENDIF
     548
     549             x = LOG( 16.0 * bond * py / 3.0 )
     550             y = 0.0
     551
     552             DO  i = 1, 6
     553                y = y + c(i) * x**(i-1)
     554             ENDDO
     555
     556             xrey = py * EXP( y )
     557
     558             IF ( radclass(j) > 0.35 )  THEN
     559                winf(j) = xrey * eta / ( 2.0 * rhoa * 0.35 )           ! in cm/s
     560             ELSE
     561                winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )    ! in cm/s
     562             ENDIF
     563
     564          ENDIF
     565
     566       ENDDO
     567
     568    END SUBROUTINE fallg
     569
     570
     571!------------------------------------------------------------------------------!
     572! Calculation of collision efficencies for the Hall kernel
     573!------------------------------------------------------------------------------!
     574    SUBROUTINE effic
     575
     576       USE arrays_3d
     577       USE cloud_parameters
     578       USE constants
     579       USE particle_attributes
     580
     581       IMPLICIT NONE
     582
     583       INTEGER ::  i, iq, ir, j, k, kk
     584
     585       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
     586
     587       LOGICAL, SAVE ::  first = .TRUE.
     588
     589       REAL ::  ek, particle_radius, pp, qq, rq
     590
     591       REAL, DIMENSION(1:21), SAVE ::  rat
     592       REAL, DIMENSION(1:15), SAVE ::  r0
     593       REAL, DIMENSION(1:15,1:21), SAVE ::  ecoll
     594
     595!
     596!--    Initial assignment of constants
     597       IF ( first )  THEN
    403598
    404599         first = .FALSE.
    405          b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
    406              -0.578878e-3,0.855176e-4,-0.327815e-5/)
    407          c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
    408              -0.542819e-1, 0.238449e-2/)
    409 
    410          eta   = 1.818e-4     !in poise = g/(cm s)
    411          xlamb = 6.62e-6    !in cm
    412 
    413          rhow = 1.0       !in g/cm**3
    414          rhoa = 1.225e-3    !in g/cm**3
    415 
    416          grav = 980.665     !in cm/s**2
    417          cunh = 1.257 * xlamb !in cm
    418          t0   = 273.15        !in K
    419          sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
    420          stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
    421          stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
    422          phy  = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa))
    423          py   = phy**(1.0/6.0)
    424 
    425       ENDIF
    426 
    427 !particle radius has to be in cm
    428       DO j = pstart, pend
    429 
    430          IF (particles(j)%radius*100.0 .le. 1.e-3) THEN
    431 
    432             winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s
    433 
    434          ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN
    435 
    436             x = log(stb*(particles(j)%radius*100.0)**3)
    437             y = 0.0
    438 
    439             DO i = 1, 7
    440                y = y + b(i) * (x**(i-1))
    441             ENDDO
    442 
    443             xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y)
    444             winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
    445 
    446          ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
    447 
    448             IF (particles(j)%radius*100.0.gt.0.35)  THEN
    449                bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma
    450             ELSE
    451                bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma
    452             ENDIF
    453 
    454             x = log(16.0*bond*py/3.0)
    455             y = 0.0
    456 
    457             DO i = 1, 6
    458                y = y + c(i) * (x**(i-1))
    459             ENDDO
    460 
    461             xrey = py*exp(y)
    462 
    463             IF (particles(j)%radius*100.0 .gt.0.35)  THEN
    464               winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
    465             ELSE
    466                winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
    467             ENDIF
    468 
    469          ENDIF
    470       ENDDO
    471       RETURN
    472    END SUBROUTINE fallg
    473 
    474 !------------------------------------------------------------------------------!
    475 ! SUBROUTINE for calculation of collision efficencies
    476 !------------------------------------------------------------------------------!
    477 
    478    SUBROUTINE effic
    479 
    480       USE arrays_3d
    481       USE constants
    482       USE cloud_parameters
    483       USE particle_attributes
    484 
    485 !collision efficiencies of hall kernel
    486       IMPLICIT NONE
    487 
    488       INTEGER :: i, ir, iq, j, k, kk
    489 
    490       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
    491 
    492       LOGICAL, SAVE ::  first = .TRUE.
    493 
    494       REAL ::  ek, particle_radius, pp, qq, rq
    495 
    496       REAL, DIMENSION(1:21), SAVE :: rat
    497       REAL, DIMENSION(1:15), SAVE :: r0
    498       REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
    499 
    500 !
    501 !--   Initial assignment of constants
    502       IF ( first )  THEN
    503 
    504          first = .FALSE.
    505          r0  = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., &
    506                  300./)
    507          rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,  &
    508                  0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
    509 
    510          ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, &
    511                         0.001,0.001,0.001,0.001,0.001,0.001/)
    512          ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, &
    513                         0.050,0.200,0.500,0.770,0.870,0.970/)
    514          ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, &
    515                         0.430,0.580,0.790,0.930,0.960,1.000/)
    516          ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, &
    517                         0.640,0.750,0.910,0.970,0.980,1.000/)
    518          ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, &
    519                         0.770,0.840,0.950,0.970,1.000,1.000/)
    520          ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, &
    521                         0.840,0.880,0.950,1.000,1.000,1.000/)
    522          ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, &
    523                         0.870,0.900,0.950,1.000,1.000,1.000/)
    524          ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, &
    525                         0.890,0.920,1.000,1.000,1.000,1.000/)
    526          ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, &
    527                         0.900,0.940,1.000,1.000,1.000,1.000/)
    528          ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, &
    529                         0.910,0.950,1.000,1.000,1.000,1.000/)
    530          ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, &
    531                         0.910,0.950,1.000,1.000,1.000,1.000/)
    532          ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, &
    533                         0.910,0.950,1.000,1.000,1.000,1.000/)
    534          ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, &
    535                         0.910,0.950,1.000,1.000,1.000,1.000/)
    536          ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, &
    537                         0.910,0.950,1.000,1.000,1.000,1.000/)
    538          ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, &
    539                         0.920,0.950,1.000,1.000,1.000,1.000/)
    540          ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, &
    541                         0.930,0.970,1.000,1.000,1.000,1.000/)
    542          ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, &
    543                         0.950,1.000,1.000,1.000,1.000,1.000/)
    544          ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, &
    545                         1.000,1.020,1.020,1.020,1.020,1.020/)
    546          ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, &
    547                         1.030,1.040,1.040,1.040,1.040,1.040/)
    548          ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, &
    549                         1.700,2.300,2.300,2.300,2.300,2.300/)
    550          ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, &
    551                         3.000,4.000,4.000,4.000,4.000,4.000/)
    552       ENDIF
    553 
    554 !
    555 !--   Calculate the radius class index of particles with respect to array r
    556       ALLOCATE( ira(pstart:pend) )
    557       DO  j = pstart, pend
    558          particle_radius = particles(j)%radius * 1.0E6
    559          DO  k = 1, 15
    560             IF ( particle_radius < r0(k) )  THEN
    561                ira(j) = k
    562                EXIT
    563             ENDIF
    564          ENDDO
    565          IF ( particle_radius >= r0(15) )  ira(j) = 16
    566       ENDDO
    567 
    568 !
    569 !--   Two-dimensional linear interpolation of the collision efficiency.
    570 !--   Radius has to be in µm
    571       DO  j = pstart, pend
    572          DO  i = pstart, j
    573 
    574             ir = ira(j)
    575 
    576             rq = particles(i)%radius / particles(j)%radius
    577 
    578 !            DO kk = 2, 21
    579 !               IF ( rq <= rat(kk) )  THEN
    580 !                  iq = kk
    581 !                  EXIT
    582 !               ENDIF
    583 !            ENDDO
    584 
    585             iq = INT( rq * 20 ) + 1
    586             iq = MAX(iq , 2)
    587 
    588             IF ( ir < 16 )  THEN
    589 
    590                IF ( ir >= 2 )  THEN
    591                   pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / &
    592                        ( r0(ir) - r0(ir-1) )
    593                   qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    594                   ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)   &
    595                             + pp * ( 1.0-qq ) * ecoll(ir,iq-1)           &
    596                             + qq * ( 1.0-pp ) * ecoll(ir-1,iq)           &
    597                             + pp * qq * ecoll(ir,iq)
    598                ELSE
    599                   qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
    600                   ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
    601                ENDIF
    602 
    603             ELSE
    604                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    605                ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
    606                ec(j,i) = MIN( ek, 1.0 )
    607             ENDIF
    608 
    609             ec(i,j) = ec(j,i)
    610             IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
    611 
    612          ENDDO
    613       ENDDO
    614 
    615       DEALLOCATE( ira )
    616 
    617    END SUBROUTINE effic
    618 
    619 !------------------------------------------------------------------------------!
    620 ! SUBROUTINE for calculation of enhancement factor collision efficencies
    621 !------------------------------------------------------------------------------!
    622    SUBROUTINE turb_enhance_eff
     600         r0  = (/ 6.0, 8.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60., &
     601                  70.0, 100.0, 150.0, 200.0, 300.0 /)
     602         rat = (/ 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, &
     603                  0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, &
     604                  1.00 /)
     605
     606         ecoll(:,1) = (/0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, &
     607                        0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001/)
     608         ecoll(:,2) = (/0.003, 0.003, 0.003, 0.004, 0.005, 0.005, 0.005, &
     609                        0.010, 0.100, 0.050, 0.200, 0.500, 0.770, 0.870, 0.970/)
     610         ecoll(:,3) = (/0.007, 0.007, 0.007, 0.008, 0.009, 0.010, 0.010, &
     611                        0.070, 0.400, 0.430, 0.580, 0.790, 0.930, 0.960, 1.000/)
     612         ecoll(:,4) = (/0.009, 0.009, 0.009, 0.012, 0.015, 0.010, 0.020, &
     613                        0.280, 0.600, 0.640, 0.750, 0.910, 0.970, 0.980, 1.000/)
     614         ecoll(:,5) = (/0.014, 0.014, 0.014, 0.015, 0.016, 0.030, 0.060, &
     615                        0.500, 0.700, 0.770, 0.840, 0.950, 0.970, 1.000, 1.000/)
     616         ecoll(:,6) = (/0.017, 0.017, 0.017, 0.020, 0.022, 0.060, 0.100, &
     617                        0.620, 0.780, 0.840, 0.880, 0.950, 1.000, 1.000, 1.000/)
     618         ecoll(:,7) = (/0.030, 0.030, 0.024, 0.022, 0.032, 0.062, 0.200, &
     619                        0.680, 0.830, 0.870, 0.900, 0.950, 1.000, 1.000, 1.000/)
     620         ecoll(:,8) = (/0.025, 0.025, 0.025, 0.036, 0.043, 0.130, 0.270, &
     621                        0.740, 0.860, 0.890, 0.920, 1.000, 1.000, 1.000, 1.000/)
     622         ecoll(:,9) = (/0.027, 0.027, 0.027, 0.040, 0.052, 0.200, 0.400, &
     623                        0.780, 0.880, 0.900, 0.940, 1.000, 1.000, 1.000, 1.000/)
     624         ecoll(:,10)= (/0.030, 0.030, 0.030, 0.047, 0.064, 0.250, 0.500, &
     625                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
     626         ecoll(:,11)= (/0.040, 0.040, 0.033, 0.037, 0.068, 0.240, 0.550, &
     627                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
     628         ecoll(:,12)= (/0.035, 0.035, 0.035, 0.055, 0.079, 0.290, 0.580, &
     629                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
     630         ecoll(:,13)= (/0.037, 0.037, 0.037, 0.062, 0.082, 0.290, 0.590, &
     631                        0.780, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
     632         ecoll(:,14)= (/0.037, 0.037, 0.037, 0.060, 0.080, 0.290, 0.580, &
     633                        0.770, 0.890, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
     634         ecoll(:,15)= (/0.037, 0.037, 0.037, 0.041, 0.075, 0.250, 0.540, &
     635                        0.760, 0.880, 0.920, 0.950, 1.000, 1.000, 1.000, 1.000/)
     636         ecoll(:,16)= (/0.037, 0.037, 0.037, 0.052, 0.067, 0.250, 0.510, &
     637                        0.770, 0.880, 0.930, 0.970, 1.000, 1.000, 1.000, 1.000/)
     638         ecoll(:,17)= (/0.037, 0.037, 0.037, 0.047, 0.057, 0.250, 0.490, &
     639                        0.770, 0.890, 0.950, 1.000, 1.000, 1.000, 1.000, 1.000/)
     640         ecoll(:,18)= (/0.036, 0.036, 0.036, 0.042, 0.048, 0.230, 0.470, &
     641                        0.780, 0.920, 1.000, 1.020, 1.020, 1.020, 1.020, 1.020/)
     642         ecoll(:,19)= (/0.040, 0.040, 0.035, 0.033, 0.040, 0.112, 0.450, &
     643                        0.790, 1.010, 1.030, 1.040, 1.040, 1.040, 1.040, 1.040/)
     644         ecoll(:,20)= (/0.033, 0.033, 0.033, 0.033, 0.033, 0.119, 0.470, &
     645                        0.950, 1.300, 1.700, 2.300, 2.300, 2.300, 2.300, 2.300/)
     646         ecoll(:,21)= (/0.027, 0.027, 0.027, 0.027, 0.027, 0.125, 0.520, &
     647                        1.400, 2.300, 3.000, 4.000, 4.000, 4.000, 4.000, 4.000/)
     648       ENDIF
     649
     650!
     651!--    Calculate the radius class index of particles with respect to array r
     652       ALLOCATE( ira(1:radius_classes) )
     653       DO  j = 1, radius_classes
     654          particle_radius = radclass(j) * 1.0E4   ! in microm
     655          DO  k = 1, 15
     656             IF ( particle_radius < r0(k) )  THEN
     657                ira(j) = k
     658                EXIT
     659             ENDIF
     660          ENDDO
     661          IF ( particle_radius >= r0(15) )  ira(j) = 16
     662       ENDDO
     663
     664!
     665!--    Two-dimensional linear interpolation of the collision efficiency.
     666!--    Radius has to be in µm
     667       DO  j = 1, radius_classes
     668          DO  i = 1, j
     669
     670             ir = ira(j)
     671             rq = radclass(i) / radclass(j)
     672             iq = INT( rq * 20 ) + 1
     673             iq = MAX( iq , 2)
     674
     675             IF ( ir < 16 )  THEN
     676                IF ( ir >= 2 )  THEN
     677                   pp = ( ( radclass(j) * 1.0E04 ) - r0(ir-1) ) / &
     678                        ( r0(ir) - r0(ir-1) )
     679                   qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     680                   ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)  &
     681                             + pp * ( 1.0-qq ) * ecoll(ir,iq-1)          &
     682                             + qq * ( 1.0-pp ) * ecoll(ir-1,iq)          &
     683                             + pp * qq * ecoll(ir,iq)
     684                ELSE
     685                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     686                   ec(j,i) = (1.0-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
     687                ENDIF
     688             ELSE
     689                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     690                ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
     691                ec(j,i) = MIN( ek, 1.0 )
     692             ENDIF
     693
     694             ec(i,j) = ec(j,i)
     695             IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
     696
     697          ENDDO
     698       ENDDO
     699
     700       DEALLOCATE( ira )
     701
     702    END SUBROUTINE effic
     703
     704
     705!------------------------------------------------------------------------------!
     706! Calculation of enhancement factor for collision efficencies due to turbulence
     707!------------------------------------------------------------------------------!
     708    SUBROUTINE turb_enhance_eff
    623709
    624710       USE constants
    … …  
    627713       USE arrays_3d
    628714
    629       IMPLICIT NONE
    630 
    631       INTEGER :: i, ik, ir, iq, j, k, kk
    632 
    633       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
    634 
    635       REAL    :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3
    636 
    637       LOGICAL, SAVE ::  first = .TRUE.
    638 
    639       REAL, DIMENSION(1:11), SAVE :: rat
    640       REAL, DIMENSION(1:7), SAVE  :: r0
    641       REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400
    642 
    643 !
    644 !--   Initial assignment of constants
    645       IF ( first )  THEN
    646 
    647          first = .FALSE.
    648 
    649          r0  = (/10., 20., 30.,40., 50., 60.,100./)
    650          rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
    651 
    652 ! 100 cm^2/s^3
    653          ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
    654          ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
    655          ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
    656          ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
    657          ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
    658          ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
    659          ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
    660          ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
    661          ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
    662          ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
    663          ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
    664 
    665 ! 400 cm^2/s^3
    666          ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
    667          ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
    668          ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
    669          ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
    670          ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
    671          ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
    672          ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
    673          ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
    674          ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
    675          ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
    676          ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
    677 
    678       ENDIF
    679 
    680 !
    681 !--   Calculate the radius class index of particles with respect to array r
    682       ALLOCATE( ira(pstart:pend) )
    683 
    684       DO  j = pstart, pend
    685          particle_radius = particles(j)%radius * 1.0E6
    686          DO  k = 1, 7
    687             IF ( particle_radius < r0(k) )  THEN
    688                ira(j) = k
    689                EXIT
    690             ENDIF
    691          ENDDO
    692          IF ( particle_radius >= r0(7) )  ira(j) = 8
    693       ENDDO
    694 
    695 ! two-dimensional linear interpolation of the collision efficiency
    696       DO j =  pstart, pend
    697          DO i = pstart, j
    698 
    699             ir = ira(j)
    700 
    701             rq = particles(i)%radius/particles(j)%radius
    702 
    703             DO kk = 2, 11
    704                IF ( rq <= rat(kk) )  THEN
    705                   iq = kk
    706                   EXIT
    707                ENDIF
    708             ENDDO
    709 
    710 ! 0  cm2/s3
    711             y1 = 1.0
    712 ! 100 cm2/s3, 400 cm2/s3
    713             IF (ir.lt.8) THEN
    714                IF (ir.ge.2) THEN
    715                   pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1))
    716                   qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
    717                   y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+  &
    718                             pp*(1.-qq)*ecoll_100(ir,iq-1)+    &
    719                             qq*(1.-pp)*ecoll_100(ir-1,iq)+    &
    720                             pp*qq*ecoll_100(ir,iq)
    721                   y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+  &
    722                             pp*(1.-qq)*ecoll_400(ir,iq-1)+    &
    723                             qq*(1.-pp)*ecoll_400(ir-1,iq)+    &
    724                             pp*qq*ecoll_400(ir,iq)
    725                ELSE
    726                   qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
    727                   y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq)
    728                   y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq)
    729                ENDIF
    730             ELSE
    731             qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
    732             y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
    733             y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
    734          ENDIF
    735 ! linear interpolation
    736 ! dissipation rate in cm2/s3
    737          x1 = 0.0
    738          x2 = 100.0
    739          x3 = 400.0
    740 
    741          IF (epsilon.le.100.) THEN
    742             ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 &
    743                     +  (epsilon-0.)/(100.-0.) * y2
    744          ELSE IF(epsilon.le.600.)THEN
    745             ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 &
    746                     +  (epsilon-100.)/(400.-100.) * y3
    747 
    748          ELSE
    749             ecf(j,i) = (600.-400.)/(100.-400.) * y2 &
    750                     +  (600.-100.)/(400.-100.) * y3
    751          ENDIF
    752 
    753          IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0
    754 
    755          ecf(i,j)=ecf(j,i)
    756       ENDDO
    757    ENDDO
    758 
    759    RETURN
    760    END SUBROUTINE turb_enhance_eff
     715       IMPLICIT NONE
     716
     717       INTEGER :: i, ik, iq, ir, j, k, kk
     718
     719       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
     720
     721       REAL ::  particle_radius, pp, qq, rq, x1, x2, x3, y1, y2, y3
     722
     723       LOGICAL, SAVE ::  first = .TRUE.
     724
     725       REAL, DIMENSION(1:11), SAVE ::  rat
     726       REAL, DIMENSION(1:7), SAVE  ::  r0
     727       REAL, DIMENSION(1:7,1:11), SAVE ::  ecoll_100, ecoll_400
     728
     729!
     730!--    Initial assignment of constants
     731       IF ( first )  THEN
     732
     733          first = .FALSE.
     734
     735          r0  = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 100.0 /)
     736          rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /)
     737!
     738!--       In 100 cm^2/s^3
     739          ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
     740          ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
     741          ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
     742          ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
     743          ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
     744          ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
     745          ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
     746          ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
     747          ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
     748          ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
     749          ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
     750!
     751!--       In 400 cm^2/s^3
     752          ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
     753          ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
     754          ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
     755          ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
     756          ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
     757          ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
     758          ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
     759          ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
     760          ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
     761          ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
     762          ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
     763
     764       ENDIF
     765
     766!
     767!--    Calculate the radius class index of particles with respect to array r0
     768       ALLOCATE( ira(1:radius_classes) )
     769
     770       DO  j = 1, radius_classes
     771          particle_radius = radclass(j) * 1.0E4    ! in microm
     772          DO  k = 1, 7
     773             IF ( particle_radius < r0(k) )  THEN
     774                ira(j) = k
     775                EXIT
     776             ENDIF
     777          ENDDO
     778          IF ( particle_radius >= r0(7) )  ira(j) = 8
     779       ENDDO
     780
     781!
     782!--    Two-dimensional linear interpolation of the collision efficiencies
     783       DO  j =  1, radius_classes
     784          DO  i = 1, j
     785
     786             ir = ira(j)
     787             rq = radclass(i) / radclass(j)
     788
     789             DO  kk = 2, 11
     790                IF ( rq <= rat(kk) )  THEN
     791                   iq = kk
     792                   EXIT
     793                ENDIF
     794             ENDDO
     795
     796             y1 = 1.0      ! 0  cm2/s3
     797!
     798!--          100 cm2/s3, 400 cm2/s3
     799             IF ( ir < 8 )  THEN
     800                IF ( ir >= 2 )  THEN
     801                   pp = ( radclass(j)*1.0E4 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
     802                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     803                   y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) +  &
     804                                pp * ( 1.0-qq ) * ecoll_100(ir,iq-1)   +  &
     805                                qq * ( 10.-pp ) * ecoll_100(ir-1,iq)   +  &
     806                                pp * qq         * ecoll_100(ir,iq)
     807                   y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) +  &
     808                                pp * ( 1.0-qq ) * ecoll_400(ir,iq-1)   +  &
     809                                qq * ( 1.0-pp ) * ecoll_400(ir-1,iq)   +  &
     810                                pp * qq         * ecoll_400(ir,iq)
     811                ELSE
     812                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     813                   y2 = ( 1.0-qq ) * ecoll_100(1,iq-1) + qq * ecoll_100(1,iq)
     814                   y3 = ( 1.0-qq ) * ecoll_400(1,iq-1) + qq * ecoll_400(1,iq)
     815                ENDIF
     816             ELSE
     817                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     818                y2 = ( 1.0-qq ) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
     819                y3 = ( 1.0-qq ) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
     820             ENDIF
     821!
     822!--          Linear interpolation of dissipation rate in cm2/s3
     823             IF ( epsilon <= 100.0 )  THEN
     824                ecf(j,i) = ( epsilon - 100.0 ) / (   0.0 - 100.0 ) * y1 &
     825                         + ( epsilon -   0.0 ) / ( 100.0 -   0.0 ) * y2
     826             ELSEIF ( epsilon <= 600.0 )  THEN
     827                ecf(j,i) = ( epsilon - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
     828                         + ( epsilon - 100.0 ) / ( 400.0 - 100.0 ) * y3
     829             ELSE
     830                ecf(j,i) = (   600.0 - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
     831                         + (   600.0 - 100.0 ) / ( 400.0 - 100.0 ) * y3
     832             ENDIF
     833
     834             IF ( ecf(j,i) < 1.0 )  ecf(j,i) = 1.0
     835
     836             ecf(i,j) = ecf(j,i)
     837
     838          ENDDO
     839       ENDDO
     840
     841    END SUBROUTINE turb_enhance_eff
    761842
    762843 END MODULE lpm_collision_kernels_mod
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |