Index: /palm/trunk/SOURCE/lpm_advec.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_advec.f90	(revision 1313)
+++ /palm/trunk/SOURCE/lpm_advec.f90	(revision 1314)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Vertical logarithmic interpolation of horizontal particle speed for particles 
+! between roughness height and first vertical grid level.
 !
 ! Former revisions:
@@ -65,7 +66,16 @@
              e_mean_int, fs_int, lagr_timescale, random_gauss, vv_int
 
+    REAL ::    height_int, height_p, log_z_z0_int, us_int, z_p, d_z_p_z0
+
     REAL ::  location(1:30,1:3)
     REAL, DIMENSION(1:30) ::  de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
 
+!
+!-- Determine height of Prandtl layer and distance between Prandtl-layer
+!-- height and horizontal mean roughness height, which are required for 
+!-- vertical logarithmic interpolation of horizontal particle speeds 
+!-- (for particles below first vertical grid level).
+    z_p      = zu(nzb+1) - zw(nzb)
+    d_z_p_z0 = 1.0 / ( z_p - z0_av_global )
 
     DO  n = 1, number_of_particles
@@ -75,67 +85,155 @@
 !--    reached
        IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 )  CYCLE
-
-!
-!--    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
+!
+!--    Determine bottom index
        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
+              + offset_ocean_nzt       ! only exact if equidistant
+!
+!--    Interpolation of the u velocity component onto particle position.  
+!--    Particles are interpolation bi-linearly in the horizontal and a 
+!--    linearly in the vertical. An exception is made for particles below
+!--    the first vertical grid level in case of a prandtl layer. In this 
+!--    case the horizontal particle velocity components are determined using
+!--    Monin-Obukhov relations (if branch). 
+!--    First, check if particle is located below first vertical grid level 
+!--    (Prandtl-layer height)
+       IF ( prandtl_layer .AND. particles(n)%z < z_p )  THEN
+!
+!--       Resolved-scale horizontal particle velocity is zero below z0. 
+          IF ( particles(n)%z < z0_av_global )  THEN
+
+             u_int = 0.0
+
+          ELSE
+!
+!--          Determine the sublayer. Further used as index.
+             height_p     = ( particles(n)%z - z0_av_global )           &
+                                  * REAL( number_of_sublayers )         &
+                                  * d_z_p_z0 
+!
+!--          Calculate LOG(z/z0) for exact particle height. Therefore,   
+!--          interpolate linearly between precalculated logarithm. 
+             log_z_z0_int = log_z_z0(INT(height_p))                     &
+                          + ( height_p - INT(height_p) )                &
+                          * ( log_z_z0(INT(height_p)+1)                 &
+                               - log_z_z0(INT(height_p))                &
+                            ) 
+!
+!--          Neutral solution is applied for all situations, e.g. also for 
+!--          unstable and stable situations. Even though this is not exact 
+!--          this saves a lot of CPU time since several calls of intrinsic 
+!--          FORTRAN procedures (LOG, ATAN) are avoided, This is justified 
+!--          as sensitivity studies revealed no significant effect of 
+!--          using the neutral solution also for un/stable situations.
+!--          Calculated left and bottom index on u grid.
+             i      = ( particles(n)%x + 0.5 * dx ) * ddx
+             j      =   particles(n)%y              * ddy
+
+             us_int = 0.5 * ( us(j,i) + us(j,i-1) )  
+
+             u_int  = -usws(j,i) / ( us_int * kappa + 1E-10 )           & 
+                      * log_z_z0_int 
+
+          ENDIF
+!
+!--    Particle above the first grid level. Bi-linear interpolation in the 
+!--    horizontal and linear interpolation in the vertical direction.
        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)  &
+!
+!--       Interpolate u velocity-component, determine left, front, bottom
+!--       index of u-array. Adopt k index from above
+          i = ( particles(n)%x + 0.5 * dx ) * ddx
+          j =   particles(n)%y * ddy
+!
+!--       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
-          u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
+
+          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
+
        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
+!--    Same procedure for interpolation of the v velocity-component.
+       IF ( prandtl_layer .AND. particles(n)%z < z_p )  THEN
+!
+!--       Resolved-scale horizontal particle velocity is zero below z0. 
+          IF ( particles(n)%z < z0_av_global )  THEN
+
+             v_int = 0.0
+
+          ELSE       
+!
+!--          Neutral solution is applied for all situations, e.g. also for 
+!--          unstable and stable situations. Even though this is not exact 
+!--          this saves a lot of CPU time since several calls of intrinsic 
+!--          FORTRAN procedures (LOG, ATAN) are avoided, This is justified 
+!--          as sensitivity studies revealed no significant effect of 
+!--          using the neutral solution also for un/stable situations.
+!--          Calculated left and bottom index on v grid.
+              i      =   particles(n)%x              * ddx
+              j      = ( particles(n)%y + 0.5 * dy ) * ddy
+
+              us_int = 0.5 * ( us(j,i) + us(j-1,i) )    
+
+              v_int  = -vsws(j,i) / ( us_int * kappa + 1E-10 )             &
+                       * log_z_z0_int 
+
+          ENDIF
+!
+!--    Particle above the first grid level. Bi-linear interpolation in the 
+!--    horizontal and linear interpolation in the vertical direction.
        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)  &
+          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
-          v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
+          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
+
        ENDIF
 
 !
-!--    Same procedure for interpolation of the w velocity-component (adopt
-!--    index i from v velocity-component)
+!--    Same procedure for interpolation of the w velocity-component
        IF ( vertical_particle_advection(particles(n)%group) )  THEN
+          i = particles(n)%x * ddx
           j = particles(n)%y * ddy
           k = particles(n)%z / dz + offset_ocean_nzt_m1
Index: /palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_init.f90	(revision 1313)
+++ /palm/trunk/SOURCE/lpm_init.f90	(revision 1314)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Vertical logarithmic interpolation of horizontal particle speed for particles 
+! between roughness height and first vertical grid level.
 !
 ! Former revisions:
@@ -110,10 +111,11 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, j, n, nn
+    INTEGER ::  i, j, k, n, nn
 #if defined( __parallel )
     INTEGER, DIMENSION(3) ::  blocklengths, displacements, types
 #endif
     LOGICAL ::  uniform_particles_l
-    REAL    ::  pos_x, pos_y, pos_z
+    REAL    ::  height_int, height_p, pos_x, pos_y, pos_z, z_p,            &
+                z0_av_local = 0.0
 
 
@@ -184,4 +186,48 @@
                  de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
                  de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ENDIF
+
+!
+!-- Allocate array required for logarithmic vertical interpolation of 
+!-- horizontal particle velocities between the surface and the first vertical
+!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of 
+!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for 
+!-- several heights. Splitting into 20 sublayers turned out to be sufficient. 
+!-- To obtain exact height levels of particles, linear interpolation is applied
+!-- (see lpm_advec.f90).
+    IF ( prandtl_layer )  THEN
+       
+       ALLOCATE ( log_z_z0(0:number_of_sublayers) ) 
+       z_p         = zu(nzb+1) - zw(nzb)
+
+!
+!--    Calculate horizontal mean value of z0 used for logartihmic
+!--    interpolation. Note: this is not exact for heterogeneous z0. 
+!--    However, sensitivity studies showed that the effect is 
+!--    negligible. 
+       z0_av_local  = SUM( z0(nys:nyn,nxl:nxr) )
+       z0_av_global = 0.0
+
+       CALL MPI_ALLREDUCE(z0_av_local, z0_av_global, 1, MPI_REAL, MPI_SUM, &
+                          comm2d, ierr )
+
+       z0_av_global = z0_av_global  / ( ( ny + 1 ) * ( nx + 1 ) )
+!
+!--    Horizontal wind speed is zero below and at z0
+       log_z_z0(0) = 0.0   
+!
+!--    Calculate vertical depth of the sublayers
+       height_int  = ( z_p - z0_av_global ) / REAL( number_of_sublayers ) 
+!
+!--    Precalculate LOG(z/z0)
+       height_p    = 0.0
+       DO  k = 1, number_of_sublayers
+
+          height_p    = height_p + height_int
+          log_z_z0(k) = LOG( height_p / z0_av_global )
+
+       ENDDO
+
+
     ENDIF
 
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 1313)
+++ /palm/trunk/SOURCE/modules.f90	(revision 1314)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! + log_z_z0, number_of_sublayers, z0_av_global 
 !
 ! Former revisions:
@@ -1428,4 +1428,5 @@
                 maximum_number_of_tailpoints = 100,                            &
                 maximum_number_of_tails = 0,                                   &
+                number_of_sublayers = 20,                                      &
                 number_of_initial_particles = 0, number_of_particles = 0,      &
                 number_of_particle_groups = 1, number_of_tails = 0,            &
@@ -1465,5 +1466,5 @@
                 sgs_wfv_part = 0.3333333, sgs_wfw_part = 0.3333333,            &
                 time_prel = 0.0, time_sort_particles = 0.0,                    &
-                time_write_particle_data = 0.0
+                time_write_particle_data = 0.0, z0_av_global
 
     REAL, DIMENSION(max_number_of_particle_groups) ::  &
@@ -1472,4 +1473,6 @@
                 psn = 9999999.9, psr = 9999999.9, pss = 9999999.9,           &
                 pst = 9999999.9, radius = 9999999.9
+
+    REAL, DIMENSION(:), ALLOCATABLE     ::  log_z_z0 
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::  particle_tail_coordinates
