Index: palm/trunk/SOURCE/lpm_set_attributes.f90
===================================================================
--- palm/trunk/SOURCE/lpm_set_attributes.f90	(revision 1321)
+++ palm/trunk/SOURCE/lpm_set_attributes.f90	(revision 1359)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-!
+! New particle structure integrated. 
+! Kind definition added to all floating point numbers.
 !
 ! Former revisions:
@@ -77,5 +78,6 @@
 
     USE particle_attributes,                                                   &
-        ONLY:  number_of_particles, offset_ocean_nzt, particles
+        ONLY:  block_offset, grid_particles, number_of_particles,              &
+               offset_ocean_nzt, particles, prt_count
 
     USE pegrid
@@ -87,7 +89,14 @@
 
     INTEGER(iwp) ::  i        !:
+    INTEGER(iwp) ::  ip       !:
     INTEGER(iwp) ::  j        !:
+    INTEGER(iwp) ::  jp       !:
     INTEGER(iwp) ::  k        !:
+    INTEGER(iwp) ::  kp       !:
     INTEGER(iwp) ::  n        !:
+    INTEGER(iwp) ::  nb       !:
+
+    INTEGER(iwp), DIMENSION(0:7) ::  start_index !:
+    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !:
 
     REAL(wp)    ::  aa        !:
@@ -101,8 +110,6 @@
     REAL(wp)    ::  pt_int_l  !:
     REAL(wp)    ::  pt_int_u  !:
-    REAL(wp)    ::  u_int     !:
     REAL(wp)    ::  u_int_l   !:
     REAL(wp)    ::  u_int_u   !:
-    REAL(wp)    ::  v_int     !:
     REAL(wp)    ::  v_int_l   !:
     REAL(wp)    ::  v_int_u   !:
@@ -113,4 +120,10 @@
     REAL(wp)    ::  y         !:
 
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_int                  !:
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_int                  !:
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  xv                     !:
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  yv                     !:
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  zv                     !:
+
     CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'start' )
 
@@ -122,77 +135,118 @@
 !--    Set particle color depending on the absolute value of the horizontal
 !--    velocity
-       DO  n = 1, number_of_particles
-!
-!--       Interpolate u velocity-component, determine left, front, bottom
-!--       index of u-array
-          i = ( particles(n)%x + 0.5 * dx ) * ddx
-          j =   particles(n)%y * ddy
-          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
-              + offset_ocean_nzt                     ! only exact if equidistant
-
-!
-!--       Interpolation of the velocity components in the xy-plane
-          x  = particles(n)%x + ( 0.5 - i ) * dx
-          y  = particles(n)%y - j * dy
-          aa = x**2          + y**2
-          bb = ( dx - x )**2 + y**2
-          cc = x**2          + ( dy - y )**2
-          dd = ( dx - x )**2 + ( dy - y )**2
-          gg = aa + bb + cc + dd
-
-          u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
-                    + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
-                    ) / ( 3.0 * gg ) - u_gtrans
-          IF ( k+1 == nzt+1 )  THEN
-             u_int = u_int_l
-          ELSE
-             u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)   &
-                       + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
-                       ) / ( 3.0 * gg ) - u_gtrans
-             u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
-                               ( u_int_u - u_int_l )
-          ENDIF
-
-!
-!--       Same procedure for interpolation of the v velocity-component (adopt
-!--       index k from u velocity-component)
-          i =   particles(n)%x * ddx
-          j = ( particles(n)%y + 0.5 * dy ) * ddy
-
-          x  = particles(n)%x - i * dx
-          y  = particles(n)%y + ( 0.5 - j ) * dy
-          aa = x**2          + y**2
-          bb = ( dx - x )**2 + y**2
-          cc = x**2          + ( dy - y )**2
-          dd = ( dx - x )**2 + ( dy - y )**2
-          gg = aa + bb + cc + dd
-
-          v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
-                 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
-                 ) / ( 3.0 * gg ) - v_gtrans
-          IF ( k+1 == nzt+1 )  THEN
-             v_int = v_int_l
-          ELSE
-             v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
-                       + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
-                       ) / ( 3.0 * gg ) - v_gtrans
-             v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
-                               ( v_int_u - v_int_l )
-          ENDIF
-
-          absuv = SQRT( u_int**2 + v_int**2 )
-
-!
-!--       Limit values by the given interval and normalize to interval [0,1]
-          absuv = MIN( absuv, color_interval(2) )
-          absuv = MAX( absuv, color_interval(1) )
-
-          absuv = ( absuv - color_interval(1) ) / &
-                  ( color_interval(2) - color_interval(1) )
-
-!
-!--       Number of available colors is defined in init_dvrp
-          particles(n)%class = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )
-
+       DO  ip = nxl, nxr
+          DO  jp = nys, nyn
+             DO  kp = nzb+1, nzt
+
+                number_of_particles = prt_count(kp,jp,ip)
+                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+                IF ( number_of_particles <= 0 )  CYCLE
+                start_index = grid_particles(kp,jp,ip)%start_index
+                end_index   = grid_particles(kp,jp,ip)%end_index
+
+                ALLOCATE( u_int(1:number_of_particles), &
+                          v_int(1:number_of_particles), &
+                          xv(1:number_of_particles),    &
+                          yv(1:number_of_particles),    &
+                          zv(1:number_of_particles) )
+
+                xv = particles(1:number_of_particles)%x
+                yv = particles(1:number_of_particles)%y
+                zv = particles(1:number_of_particles)%z
+
+                DO  nb = 0,7
+
+                   i = ip
+                   j = jp + block_offset(nb)%j_off
+                   k = kp + block_offset(nb)%k_off
+
+                   DO  n = start_index(nb), end_index(nb)
+!
+!--                   Interpolation of the velocity components in the xy-plane
+                      x  = xv(n) + ( 0.5_wp - i ) * dx
+                      y  = yv(n) - j * dy
+                      aa = x**2          + y**2
+                      bb = ( dx - x )**2 + y**2
+                      cc = x**2          + ( dy - y )**2
+                      dd = ( dx - x )**2 + ( dy - y )**2
+                      gg = aa + bb + cc + dd
+
+                      u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) *     &
+                                  u(k,j,i+1) + ( gg - cc ) * u(k,j+1,i) +      &
+                                  ( gg - dd ) * u(k,j+1,i+1)                   &
+                                ) / ( 3.0_wp * gg ) - u_gtrans
+
+                      IF ( k+1 == nzt+1 )  THEN
+                         u_int(n) = u_int_l
+                      ELSE
+                         u_int_u = ( ( gg - aa ) * u(k+1,j,i)   + ( gg - bb ) *  &
+                                     u(k+1,j,i+1) + ( gg - cc ) * u(k+1,j+1,i) + &
+                                     ( gg - dd ) * u(k+1,j+1,i+1)                &
+                                   ) / ( 3.0_wp * gg ) - u_gtrans
+                         u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *         &
+                                           ( u_int_u - u_int_l )
+                      ENDIF
+
+                   ENDDO
+
+                   i = ip + block_offset(nb)%i_off
+                   j = jp
+                   k = kp + block_offset(nb)%k_off
+
+                   DO  n = start_index(nb), end_index(nb)
+!
+!--                   Same procedure for interpolation of the v velocity-component
+                      x  = xv(n) - i * dx
+                      y  = yv(n) + ( 0.5_wp - j ) * dy
+                      aa = x**2          + y**2
+                      bb = ( dx - x )**2 + y**2
+                      cc = x**2          + ( dy - y )**2
+                      dd = ( dx - x )**2 + ( dy - y )**2
+                      gg = aa + bb + cc + dd
+
+                      v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) *     &
+                                  v(k,j,i+1) + ( gg - cc ) * v(k,j+1,i) +      &
+                                  ( gg - dd ) * v(k,j+1,i+1)                   &
+                                ) / ( 3.0_wp * gg ) - v_gtrans
+
+                      IF ( k+1 == nzt+1 )  THEN
+                         v_int(n) = v_int_l
+                      ELSE
+                         v_int_u  = ( ( gg - aa ) * v(k+1,j,i) + ( gg - bb ) *    &
+                                      v(k+1,j,i+1) + ( gg - cc ) * v(k+1,j+1,i) + &
+                                      ( gg - dd ) * v(k+1,j+1,i+1)                &
+                                    ) / ( 3.0_wp * gg ) - v_gtrans
+                         v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz *         &
+                                           ( v_int_u - v_int_l )
+                      ENDIF
+
+                   ENDDO
+
+                ENDDO
+
+                DO  n = 1, number_of_particles
+
+                   absuv = SQRT( u_int(n)**2 + v_int(n)**2 )
+
+!
+!--                Limit values by the given interval and normalize to
+!--                interval [0,1]
+                   absuv = MIN( absuv, color_interval(2) )
+                   absuv = MAX( absuv, color_interval(1) )
+
+                   absuv = ( absuv - color_interval(1) ) / &
+                           ( color_interval(2) - color_interval(1) )
+
+!
+!--                Number of available colors is defined in init_dvrp
+                   particles(n)%class = 1 + absuv *                            &
+                                            ( dvrp_colortable_entries_prt - 1 )
+
+                ENDDO
+
+                DEALLOCATE( u_int, v_int, xv, yv, zv )
+
+             ENDDO
+          ENDDO
        ENDDO
 
@@ -204,5 +258,5 @@
 !--    (This is also done in flow_statistics, but flow_statistics is called
 !--    after this routine.)
-       sums_l(:,4,0) = 0.0
+       sums_l(:,4,0) = 0.0_wp
        DO  i = nxl, nxr
           DO  j =  nys, nyn
@@ -214,56 +268,81 @@
 
 #if defined( __parallel )
-
        IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
        CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, &
                            MPI_REAL, MPI_SUM, comm2d, ierr )
-
 #else
-
        sums(:,4) = sums_l(:,4,0)
-
 #endif
-
        sums(:,4) = sums(:,4) / ngp_2dh(0)
 
-       DO  n = 1, number_of_particles
-!
-!--       Interpolate temperature to the current particle position
-          i = particles(n)%x * ddx
-          j = particles(n)%y * ddy
-          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
-              + offset_ocean_nzt                     ! only exact if equidistant
-
-          x  = particles(n)%x - i * dx
-          y  = particles(n)%y - j * dy
-          aa = x**2          + y**2
-          bb = ( dx - x )**2 + y**2
-          cc = x**2          + ( dy - y )**2
-          dd = ( dx - x )**2 + ( dy - y )**2
-          gg = aa + bb + cc + dd
-
-          pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) * pt(k,j,i+1)   &
-                     + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
-                     ) / ( 3.0 * gg ) - sums(k,4)
-
-          pt_int_u = ( ( gg-aa ) * pt(k+1,j,i)   + ( gg-bb ) * pt(k+1,j,i+1)   &
-                     + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
-                     ) / ( 3.0 * gg ) - sums(k,4)
-
-          pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
-                              ( pt_int_u - pt_int_l )
-
-!
-!--       Limit values by the given interval and normalize to interval [0,1]
-          pt_int = MIN( pt_int, color_interval(2) )
-          pt_int = MAX( pt_int, color_interval(1) )
-
-          pt_int = ( pt_int - color_interval(1) ) / &
-                   ( color_interval(2) - color_interval(1) )
-
-!
-!--       Number of available colors is defined in init_dvrp
-          particles(n)%class = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
-
+       DO  ip = nxl, nxr
+          DO  jp = nys, nyn
+             DO  kp = nzb+1, nzt
+
+                number_of_particles = prt_count(kp,jp,ip)
+                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+                IF ( number_of_particles <= 0 )  CYCLE
+                start_index = grid_particles(kp,jp,ip)%start_index
+                end_index   = grid_particles(kp,jp,ip)%end_index
+
+                ALLOCATE( xv(1:number_of_particles), &
+                          yv(1:number_of_particles), &
+                          zv(1:number_of_particles) )
+
+                xv = particles(1:number_of_particles)%x
+                yv = particles(1:number_of_particles)%y
+                zv = particles(1:number_of_particles)%z
+
+                DO  nb = 0,7
+
+                   i = ip + block_offset(nb)%i_off
+                   j = jp + block_offset(nb)%j_off
+                   k = kp + block_offset(nb)%k_off
+
+                   DO  n = start_index(nb), end_index(nb)
+!
+!--                   Interpolate temperature to the current particle position
+                      x  = xv(n) - i * dx
+                      y  = yv(n) - j * dy
+                      aa = x**2          + y**2
+                      bb = ( dx - x )**2 + y**2
+                      cc = x**2          + ( dy - y )**2
+                      dd = ( dx - x )**2 + ( dy - y )**2
+                      gg = aa + bb + cc + dd
+
+                      pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) *   &
+                                   pt(k,j,i+1) + ( gg - cc ) * pt(k,j+1,i) +   &
+                                   ( gg - dd ) * pt(k,j+1,i+1)                 &
+                                 ) / ( 3.0_wp * gg ) - sums(k,4)
+
+                      pt_int_u = ( ( gg - aa ) * pt(k+1,j,i) + ( gg - bb ) *     &
+                                   pt(k+1,j,i+1) + ( gg - cc ) * pt(k+1,j+1,i) + &
+                                   ( gg - dd ) * pt(k+1,j+1,i+1)                 &
+                                 ) / ( 3.0_wp * gg ) - sums(k,4)
+
+                      pt_int = pt_int_l + ( zv(n) - zu(k) ) / dz *    &
+                                          ( pt_int_u - pt_int_l )
+
+!
+!--                   Limit values by the given interval and normalize to
+!--                   interval [0,1]
+                      pt_int = MIN( pt_int, color_interval(2) )
+                      pt_int = MAX( pt_int, color_interval(1) )
+
+                      pt_int = ( pt_int - color_interval(1) ) /                &
+                               ( color_interval(2) - color_interval(1) )
+
+!
+!--                   Number of available colors is defined in init_dvrp
+                      particles(n)%class = 1 + pt_int *                        &
+                                           ( dvrp_colortable_entries_prt - 1 )
+
+                   ENDDO
+                ENDDO
+
+                DEALLOCATE( xv, yv, zv )
+
+             ENDDO
+          ENDDO
        ENDDO
 
@@ -272,19 +351,32 @@
 !--    Set particle color depending on the height above the bottom
 !--    boundary (z=0)
-       DO  n = 1, number_of_particles
-
-          height = particles(n)%z
-!
-!--       Limit values by the given interval and normalize to interval [0,1]
-          height = MIN( height, color_interval(2) )
-          height = MAX( height, color_interval(1) )
-
-          height = ( height - color_interval(1) ) / &
-                   ( color_interval(2) - color_interval(1) )
-
-!
-!--       Number of available colors is defined in init_dvrp
-          particles(n)%class = 1 + height * ( dvrp_colortable_entries_prt - 1 )
-
+       DO  ip = nxl, nxr
+          DO  jp = nys, nyn
+             DO  kp = nzb+1, nzt
+
+                number_of_particles = prt_count(kp,jp,ip)
+                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+                IF ( number_of_particles <= 0 )  CYCLE
+                DO  n = 1, number_of_particles
+
+                   height = particles(n)%z
+!
+!--                Limit values by the given interval and normalize to
+!--                interval [0,1]
+                   height = MIN( height, color_interval(2) )
+                   height = MAX( height, color_interval(1) )
+
+                   height = ( height - color_interval(1) ) / &
+                            ( color_interval(2) - color_interval(1) )
+
+!
+!--                Number of available colors is defined in init_dvrp
+                   particles(n)%class = 1 + height *                           &
+                                            ( dvrp_colortable_entries_prt - 1 )
+
+                ENDDO
+
+             ENDDO
+          ENDDO
        ENDDO
 
@@ -295,44 +387,74 @@
     IF ( particle_dvrpsize == 'absw' )  THEN
 
-       DO  n = 1, number_of_particles
-!
-!--       Interpolate w-component to the current particle position
-          i = particles(n)%x * ddx
-          j = particles(n)%y * ddy
-          k = particles(n)%z / dz
-
-          x  = particles(n)%x - i * dx
-          y  = particles(n)%y - j * dy
-          aa = x**2          + y**2
-          bb = ( dx - x )**2 + y**2
-          cc = x**2          + ( dy - y )**2
-          dd = ( dx - x )**2 + ( dy - y )**2
-          gg = aa + bb + cc + dd
-
-          w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
-                    + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
-                    ) / ( 3.0 * gg )
-
-          IF ( k+1 == nzt+1 )  THEN
-             w_int = w_int_l
-          ELSE
-             w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + ( gg-bb ) * w(k+1,j,i+1)   &
-                       + ( gg-cc ) * w(k+1,j+1,i) + ( gg-dd ) * w(k+1,j+1,i+1) &
-                       ) / ( 3.0 * gg )
-             w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
-                               ( w_int_u - w_int_l )
-          ENDIF
-
-!
-!--       Limit values by the given interval and normalize to interval [0,1]
-          w_int = ABS( w_int )
-          w_int = MIN( w_int, dvrpsize_interval(2) )
-          w_int = MAX( w_int, dvrpsize_interval(1) )
-
-          w_int = ( w_int - dvrpsize_interval(1) ) / &
-                  ( dvrpsize_interval(2) - dvrpsize_interval(1) )
-
-          particles(n)%dvrp_psize = ( 0.25 + w_int * 0.6 ) * dx
-
+       DO  ip = nxl, nxr
+          DO  jp = nys, nyn
+             DO  kp = nzb+1, nzt
+
+                number_of_particles = prt_count(kp,jp,ip)
+                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+                IF ( number_of_particles <= 0 )  CYCLE
+                start_index = grid_particles(kp,jp,ip)%start_index
+                end_index   = grid_particles(kp,jp,ip)%end_index
+
+                ALLOCATE( xv(1:number_of_particles), &
+                          yv(1:number_of_particles) )
+
+                xv = particles(1:number_of_particles)%x
+                yv = particles(1:number_of_particles)%y
+                zv = particles(1:number_of_particles)%z
+
+                DO  nb = 0,7
+
+                   i = ip + block_offset(nb)%i_off
+                   j = jp + block_offset(nb)%j_off
+                   k = kp-1
+
+                   DO  n = start_index(nb), end_index(nb)
+!
+!--                   Interpolate w-component to the current particle position
+                      x  = xv(n) - i * dx
+                      y  = yv(n) - j * dy
+                      aa = x**2          + y**2
+                      bb = ( dx - x )**2 + y**2
+                      cc = x**2          + ( dy - y )**2
+                      dd = ( dx - x )**2 + ( dy - y )**2
+                      gg = aa + bb + cc + dd
+
+                      w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) *     &
+                                  w(k,j,i+1) + ( gg - cc ) * w(k,j+1,i) +      &
+                                  ( gg - dd ) * w(k,j+1,i+1)                   &
+                                ) / ( 3.0_wp * gg )
+
+                      IF ( k+1 == nzt+1 )  THEN
+                         w_int = w_int_l
+                      ELSE
+                         w_int_u = ( ( gg - aa ) * w(k+1,j,i)   + ( gg - bb ) *  &
+                                     w(k+1,j,i+1) + ( gg - cc ) * w(k+1,j+1,i) + &
+                                     ( gg - dd ) * w(k+1,j+1,i+1)                &
+                                   ) / ( 3.0_wp * gg )
+                         w_int = w_int_l + ( zv(n) - zw(k) ) / dz *     &
+                                           ( w_int_u - w_int_l )
+                      ENDIF
+
+!
+!--                   Limit values by the given interval and normalize to
+!--                   interval [0,1]
+                      w_int = ABS( w_int )
+                      w_int = MIN( w_int, dvrpsize_interval(2) )
+                      w_int = MAX( w_int, dvrpsize_interval(1) )
+
+                      w_int = ( w_int - dvrpsize_interval(1) ) / &
+                              ( dvrpsize_interval(2) - dvrpsize_interval(1) )
+
+                      particles(n)%dvrp_psize = ( 0.25_wp + w_int * 0.6_wp ) * &
+                                                dx
+
+                   ENDDO
+                ENDDO
+
+                DEALLOCATE( xv, yv, zv )
+
+             ENDDO
+          ENDDO
        ENDDO
 
