Index: palm/trunk/SOURCE/advec_particles.f90
===================================================================
--- palm/trunk/SOURCE/advec_particles.f90	(revision 791)
+++ palm/trunk/SOURCE/advec_particles.f90	(revision 792)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! particle arrays (particles, particles_temp) implemented as pointers in
+! order to speed up sorting (see routine sort_particles)
 !
 ! Former revisions:
@@ -3995,5 +3996,5 @@
     INTEGER ::  i, ilow, j, k, n
 
-    TYPE(particle_type), DIMENSION(1:number_of_particles) ::  particles_temp
+    TYPE(particle_type), DIMENSION(:), POINTER ::  particles_temp
 
 
@@ -4001,6 +4002,26 @@
 
 !
-!-- Initialize the array used for counting and indexing the particles
-    prt_count = 0
+!-- Initialize counters and set pointer of the temporary array into which
+!-- particles are sorted to free memory
+    prt_count  = 0
+    sort_count = sort_count +1
+
+    SELECT CASE ( MOD( sort_count, 2 ) )
+
+       CASE ( 0 )
+
+          particles_temp => part_1
+          part_1 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
+                                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
+                                  0.0, 0, 0, 0, 0 )
+
+       CASE ( 1 )
+
+          particles_temp => part_2
+          part_2 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
+                                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
+                                  0.0, 0, 0, 0, 0 )
+
+    END SELECT
 
 !
@@ -4055,5 +4076,17 @@
     ENDDO
 
-    particles(1:number_of_particles) = particles_temp
+!
+!-- Redirect the particles pointer to the sorted array
+    SELECT CASE ( MOD( sort_count, 2 ) )
+
+       CASE ( 0 )
+
+          particles => part_1
+
+       CASE ( 1 )
+
+          particles => part_2
+
+    END SELECT
 
 !
Index: palm/trunk/SOURCE/init_particles.f90
===================================================================
--- palm/trunk/SOURCE/init_particles.f90	(revision 791)
+++ palm/trunk/SOURCE/init_particles.f90	(revision 792)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! array particles implemented as pointer
 !
 ! Former revisions:
@@ -190,6 +191,19 @@
        ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
                  prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
-                 particle_mask(maximum_number_of_particles),         &
-                 particles(maximum_number_of_particles) )
+                 particle_mask(maximum_number_of_particles),     &
+                 part_1(maximum_number_of_particles),            &
+                 part_2(maximum_number_of_particles) )
+
+       part_1 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
+                               0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
+                               0.0, 0, 0, 0, 0 )
+
+       part_2 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
+                               0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
+                               0.0, 0, 0, 0, 0 )
+
+       sort_count = 0
+
+       particles => part_1
 
        READ ( 90 )  prt_count, prt_start_index
@@ -214,6 +228,11 @@
        ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
                  prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
-                 particle_mask(maximum_number_of_particles),         &
-                 particles(maximum_number_of_particles) )
+                 particle_mask(maximum_number_of_particles),     &
+                 part_1(maximum_number_of_particles),            &
+                 part_2(maximum_number_of_particles)  )
+
+       particles => part_1
+
+       sort_count = 0
 
 !
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 791)
+++ palm/trunk/SOURCE/modules.f90	(revision 792)
@@ -5,5 +5,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! particle arrays (particles, parrticles_temp) implemented as pointers,
+! +particles_1, particles_2, sort_count
 !
 ! Former revisions:
@@ -1166,5 +1167,6 @@
                 offset_ocean_nzt_m1 = 0, particles_per_point = 1,              &
                 particle_file_count = 0, skip_particles_for_tail = 100,        &
-                total_number_of_particles, total_number_of_tails = 0
+                sort_count = 0, total_number_of_particles,                     &
+                total_number_of_tails = 0
 
     INTEGER, PARAMETER ::  max_number_of_particle_groups = 10
@@ -1212,6 +1214,7 @@
     END TYPE particle_type
 
-    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  initial_particles, &
-                                                       particles
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  initial_particles
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE, TARGET ::  part_1, part_2
+    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles
 
     TYPE particle_groups_type
Index: palm/trunk/SOURCE/wang_kernel.f90
===================================================================
--- palm/trunk/SOURCE/wang_kernel.f90	(revision 791)
+++ palm/trunk/SOURCE/wang_kernel.f90	(revision 792)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! speed optimizations + formatting
 !
 ! Former revisions:
@@ -29,25 +29,34 @@
     PRIVATE
 
-    PUBLIC TurbSD, turb_enhance_eff, effic, fallg, PHI, ZHI, colker
-
-    INTEGER, SAVE :: ip, jp, kp, pstart, pend, psum
-    REAL, SAVE :: epsilon, urms
-    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: gck, ec, ecf
-    REAL, DIMENSION(:),   ALLOCATABLE, SAVE :: winf
-    REAL, SAVE :: eps2, urms2
+    PUBLIC  colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi 
+
+    INTEGER, SAVE ::  ip, jp, kp, pend, pstart, psum
+
+    REAL, SAVE :: epsilon, eps2, urms, urms2
+
+    REAL, DIMENSION(:),   ALLOCATABLE, SAVE ::  winf
+    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  ec, ecf, gck
 
 !
 !-- Public interfaces
-    INTERFACE TurbSD
-       MODULE PROCEDURE TurbSD
-    END INTERFACE TurbSD
-
-    INTERFACE PHI
-       MODULE PROCEDURE PHI
-    END INTERFACE PHI
-
-    INTERFACE ZHI
-       MODULE PROCEDURE ZHI
-    END INTERFACE ZHI
+    INTERFACE colker
+       MODULE PROCEDURE colker
+    END INTERFACE colker
+
+    INTERFACE effic
+       MODULE PROCEDURE effic
+    END INTERFACE effic
+
+    INTERFACE fallg
+       MODULE PROCEDURE fallg
+    END INTERFACE fallg
+
+    INTERFACE phi
+       MODULE PROCEDURE phi
+    END INTERFACE phi
+
+    INTERFACE turbsd
+       MODULE PROCEDURE turbsd
+    END INTERFACE turbsd
 
     INTERFACE turb_enhance_eff
@@ -55,15 +64,7 @@
     END INTERFACE turb_enhance_eff
 
-    INTERFACE effic
-       MODULE PROCEDURE effic
-    END INTERFACE effic
-
-    INTERFACE fallg
-       MODULE PROCEDURE fallg
-    END INTERFACE fallg
-
-    INTERFACE colker
-       MODULE PROCEDURE colker
-    END INTERFACE colker
+    INTERFACE zhi
+       MODULE PROCEDURE zhi
+    END INTERFACE zhi
 
     CONTAINS
@@ -72,17 +73,21 @@
 ! SUBROUTINE for calculation of collision kernel
 !------------------------------------------------------------------------------!
-    SUBROUTINE colker(i1,j1,k1,kernel)
+    SUBROUTINE colker( i1, j1, k1, kernel )
 
        USE arrays_3d
        USE cloud_parameters
        USE constants
+       USE cpulog
        USE indices
+       USE interfaces
        USE particle_attributes
 
        IMPLICIT NONE
 
-       INTEGER :: i,i1,j,j1,k1
+       INTEGER ::  i, i1, j, j1, k1
 
        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
+
+!       CALL cpu_log( log_point_s(46), 'colker', 'start' )
 
        ip = i1
@@ -94,21 +99,23 @@
        psum   = prt_count(kp,jp,ip)
 
-       ALLOCATE(winf(pstart:pend), ec(pstart:pend,pstart:pend))
-
-       IF ( turbulence_effects_on_collision ) THEN
-
-          ALLOCATE(gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend))
-
-          epsilon = diss(kp,jp,ip) * 10000.0  !dissipation rate in cm**2/s**-3
-          urms     = 202.0 * (epsilon/400.0)**(1.0/3.0)
-
-          IF (epsilon <= 0.001) THEN
+       ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
+
+       IF ( turbulence_effects_on_collision )  THEN
+
+          ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
+
+          epsilon = diss(kp,jp,ip) * 1.0E5  !dissipation rate in cm**2/s**-3
+          urms     = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 )
+
+          IF ( epsilon <= 0.001 )  THEN
 
              CALL fallg
              CALL effic
 
-             DO j = pstart, pend, 1
-                DO i =  pstart, pend, 1
-                   kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * ec(i,j) * abs(winf(j)-winf(i))
+             DO  j = pstart, pend
+                DO  i =  pstart, pend
+                   kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
+                                        particles(i)%radius * 100.0 )**2 &
+                                    * ec(i,j) * ABS( winf(j) - winf(i) )
                 ENDDO
              ENDDO
@@ -116,13 +123,14 @@
           ELSE
 
-             CALL TurbSD
+             CALL turbsd
              CALL turb_enhance_eff
              CALL effic
 
-             DO j = pstart, pend, 1
-                DO i =  pstart, pend, 1
+             DO  j = pstart, pend, 1
+                DO  i =  pstart, pend, 1
                    kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
                 ENDDO
              ENDDO
+
           ENDIF
 
@@ -131,10 +139,16 @@
        ELSE
 
+!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )
           CALL fallg
+!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )
+!          CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )
           CALL effic
-
-          DO j = pstart, pend, 1
-             DO i =  pstart, pend, 1
-                kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * ec(i,j) * abs(winf(j)-winf(i))
+!          CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )
+
+          DO  j = pstart, pend
+             DO  i =  pstart, pend
+                kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
+                                     particles(i)%radius * 100.0 )**2 &
+                                 * ec(i,j) * ABS( winf(j) - winf(i) )
              ENDDO
           ENDDO
@@ -142,7 +156,7 @@
        ENDIF
 
-       DEALLOCATE(winf, ec)
-
-       RETURN
+       DEALLOCATE( ec, winf )
+
+!       CALL cpu_log( log_point_s(46), 'colker', 'stop' )
 
     END SUBROUTINE colker
@@ -151,6 +165,8 @@
 ! SUBROUTINE for calculation of w, g and gck
 !------------------------------------------------------------------------------!
-    SUBROUTINE TurbSD
-!from Aayala 2008b, page 37ff, necessary input parameter water density, radii of droplets, air density, air viscosity, turbulent dissipation rate, taylor microscale reynolds number, gravitational acceleration
+    SUBROUTINE turbsd
+! from Aayala 2008b, page 37ff, necessary input parameter water density, radii
+! of droplets, air density, air viscosity, turbulent dissipation rate,
+! taylor microscale reynolds number, gravitational acceleration
        USE constants
        USE cloud_parameters
@@ -301,6 +317,4 @@
        ENDDO
 
-
-       RETURN
     END SUBROUTINE TurbSD
 
@@ -317,5 +331,5 @@
 
        PHI = 1./aa1 - vsett/2.0/b/aa1**2.0  !in s
-       RETURN
+
     END FUNCTION PHI
 
@@ -339,5 +353,4 @@
              + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2.0 + vsett2/aa1**2.0)   &
                                                            * 1./2./b/aa3             ! in s**2
-        RETURN
     END FUNCTION ZHI
 
@@ -430,5 +443,5 @@
          ENDIF
       ENDDO
-      RETURN
+
    END SUBROUTINE fallg
 
@@ -439,90 +452,146 @@
    SUBROUTINE effic
 
-       USE constants
-       USE cloud_parameters
-       USE particle_attributes
-       USE arrays_3d
-
-!collision efficiencies of hall kernel
+      USE arrays_3d
+      USE constants
+      USE cloud_parameters
+      USE particle_attributes
+
       IMPLICIT NONE
 
       INTEGER :: i, ir, iq, j, k, kk
-      REAL    :: ek, rq, pp, qq
-
-      REAL, DIMENSION(1:21) :: rat
-      REAL, DIMENSION(1:15) :: r0
-      REAL, DIMENSION(1:15,1:21) :: ecoll
-
-      r0 = (/6.,8.,10.,15.,20.,25.,30.,40.,50., &
-              60.,70.,100.,150.,200.,300./)
-      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,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
-
-      ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001/)
-      ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100,0.050,0.200,0.500,0.770,0.870,0.970/)
-      ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400,0.430,0.580,0.790,0.930,0.960,1.000/)
-      ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600,0.640,0.750,0.910,0.970,0.980,1.000/)
-      ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700,0.770,0.840,0.950,0.970,1.000,1.000/)
-      ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780,0.840,0.880,0.950,1.000,1.000,1.000/)
-      ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830,0.870,0.900,0.950,1.000,1.000,1.000/)
-      ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860,0.890,0.920,1.000,1.000,1.000,1.000/)
-      ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880,0.900,0.940,1.000,1.000,1.000,1.000/)
-      ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900,0.910,0.950,1.000,1.000,1.000,1.000/)
-      ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900,0.910,0.950,1.000,1.000,1.000,1.000/)
-      ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900,0.910,0.950,1.000,1.000,1.000,1.000/)
-      ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900,0.910,0.950,1.000,1.000,1.000,1.000/)
-      ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890,0.910,0.950,1.000,1.000,1.000,1.000/)
-      ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880,0.920,0.950,1.000,1.000,1.000,1.000/)
-      ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880,0.930,0.970,1.000,1.000,1.000,1.000/)
-      ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890,0.950,1.000,1.000,1.000,1.000,1.000/)
-      ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920,1.000,1.020,1.020,1.020,1.020,1.020/)
-      ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010,1.030,1.040,1.040,1.040,1.040,1.040/)
-      ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300,1.700,2.300,2.300,2.300,2.300,2.300/)
-      ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300,3.000,4.000,4.000,4.000,4.000,4.000/)
-
-! two-dimensional linear interpolation of the collision efficiency
-! radius has to be in Âµm
-
-      DO j = pstart, pend
-         DO i = pstart, j
-
-            DO k = 2, 15
-               IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN
-                  ir = k
-               ELSEIF ((particles(j)%radius*1.0E06).gt.r0(15)) THEN
-                  ir = 16
-               ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN
-                  ir = 1
-               ENDIF
-            ENDDO
-
-            rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06)
-
-            DO kk = 2, 21
-               IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
-            ENDDO
-
-            IF (ir.lt.16) THEN
-               IF (ir.ge.2) THEN
-                  pp =((particles(j)%radius * 1.0E06) - r0(ir-1)) / (r0(ir)-r0(ir-1))
-                  qq =(rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-                  ec(j,i) = (1.-pp) * (1.-qq) * ecoll(ir-1,iq-1) +  &
-                          pp * (1.-qq) * ecoll(ir,iq-1) +         &
-                          qq * (1.-pp) * ecoll(ir-1,iq) +         &
-                          pp * qq * ecoll(ir,iq)
+
+      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
+
+      LOGICAL, SAVE ::  first = .TRUE.
+
+      REAL ::  ek, particle_radius, pp, qq, rq
+
+      REAL, DIMENSION(1:21), SAVE :: rat
+      REAL, DIMENSION(1:15), SAVE :: r0
+      REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
+
+!
+!--   Initial assignment of constants
+      IF ( first )  THEN
+
+         first = .FALSE.
+         r0  = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., &
+                 300./)
+         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,  &
+                 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
+
+         ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, &
+                        0.001,0.001,0.001,0.001,0.001,0.001/)
+         ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, &
+                        0.050,0.200,0.500,0.770,0.870,0.970/)
+         ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, &
+                        0.430,0.580,0.790,0.930,0.960,1.000/)
+         ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, &
+                        0.640,0.750,0.910,0.970,0.980,1.000/)
+         ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, &
+                        0.770,0.840,0.950,0.970,1.000,1.000/)
+         ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, &
+                        0.840,0.880,0.950,1.000,1.000,1.000/)
+         ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, &
+                        0.870,0.900,0.950,1.000,1.000,1.000/)
+         ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, &
+                        0.890,0.920,1.000,1.000,1.000,1.000/)
+         ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, &
+                        0.900,0.940,1.000,1.000,1.000,1.000/)
+         ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, &
+                        0.910,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, &
+                        0.910,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, &
+                        0.910,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, &
+                        0.910,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, &
+                        0.910,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, &
+                        0.920,0.950,1.000,1.000,1.000,1.000/)
+         ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, &
+                        0.930,0.970,1.000,1.000,1.000,1.000/)
+         ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, &
+                        0.950,1.000,1.000,1.000,1.000,1.000/)
+         ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, &
+                        1.000,1.020,1.020,1.020,1.020,1.020/)
+         ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, &
+                        1.030,1.040,1.040,1.040,1.040,1.040/)
+         ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, &
+                        1.700,2.300,2.300,2.300,2.300,2.300/)
+         ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, &
+                        3.000,4.000,4.000,4.000,4.000,4.000/)
+      ENDIF
+
+!
+!--   Calculate the radius class index of particles with respect to array r
+      ALLOCATE( ira(pstart:pend) )
+      DO  j = pstart, pend
+         particle_radius = particles(j)%radius * 1.0E6
+         DO  k = 1, 15
+            IF ( particle_radius < r0(k) )  THEN
+               ira(j) = k
+               EXIT
+            ENDIF
+         ENDDO
+         IF ( particle_radius >= r0(15) )  ira(j) = 16
+      ENDDO
+
+!
+!--   Two-dimensional linear interpolation of the collision efficiency.
+!--   Radius has to be in Âµm
+      DO  j = pstart, pend
+         DO  i = pstart, j
+
+!            DO k = 2, 15
+!               IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN
+!                  ir = k
+!               ELSEIF ((particles(j)%radius*1.0E06).gt.r0(15)) THEN
+!                  ir = 16
+!               ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN
+!                  ir = 1
+!               ENDIF
+!            ENDDO
+
+            ir = ira(j)
+
+            rq = particles(i)%radius / particles(j)%radius
+
+!            DO kk = 2, 21
+!               IF ( rq .le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
+!            ENDDO
+
+            iq = INT( rq * 20 ) + 1
+
+            IF ( ir < 16 )  THEN
+
+               IF ( ir >= 2 )  THEN
+                  pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / &
+                       ( r0(ir) - r0(ir-1) )
+                  qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                  ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)   &
+                            + pp * ( 1.0-qq ) * ecoll(ir,iq-1)           &
+                            + qq * ( 1.0-pp ) * ecoll(ir-1,iq)           &
+                            + pp * qq * ecoll(ir,iq)
                ELSE
                   qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
                   ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
                ENDIF
+
             ELSE
-               qq = (rq - rat(iq-1))/(rat(iq)-rat(iq-1))
-               ek = (1.-qq)* ecoll(15,iq-1) + qq * ecoll(15,iq)
-               ec(j,i) = min(ek,1.0)
+               qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+               ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
+               ec(j,i) = MIN( ek, 1.0 )
             ENDIF
+
             ec(i,j) = ec(j,i)
-            IF (ec(i,j).lt.1.e-20) stop 99
+            IF ( ec(i,j) < 1.0E-20 )  STOP 99
+
          ENDDO
       ENDDO
-      RETURN
+
+      DEALLOCATE( ira )
+
    END SUBROUTINE effic
 
