Index: /palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/average_3d_data.f90	(revision 977)
+++ /palm/trunk/SOURCE/average_3d_data.f90	(revision 978)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! +z0h_av
 !
 ! Former revisions:
@@ -284,4 +285,11 @@
              ENDDO
 
+          CASE ( 'z0h*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   z0h_av(j,i) = z0h_av(j,i) / REAL( average_count_3d )
+                ENDDO
+             ENDDO
+
           CASE DEFAULT
 !
Index: /palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/boundary_conds.f90	(revision 977)
+++ /palm/trunk/SOURCE/boundary_conds.f90	(revision 978)
@@ -4,5 +4,8 @@
 ! Current revisions:
 ! -----------------
-!
+! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
+! Outflow boundary conditions for the velocity components can be set to Neumann
+! conditions or to radiation conditions with a horizontal averaged phase
+! velocity. 
 ! 
 !  
@@ -190,5 +193,4 @@
           q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
 
-
        ENDIF
 
@@ -198,8 +200,15 @@
 !--    Since in prognostic_equations (cache optimized version) these levels are
 !--    handled as a prognostic level, boundary values have to be restored here.
+!--    For the SGS-TKE, Neumann boundary conditions are used at the inflow.
        IF ( inflow_s )  THEN
           v_p(:,nys,:) = v_p(:,nys-1,:)
+          IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
+       ELSEIF ( inflow_n )  THEN
+          IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
        ELSEIF ( inflow_l ) THEN
           u_p(:,:,nxl) = u_p(:,:,nxl-1)
+          IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
+       ELSEIF ( inflow_r )  THEN
+          IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
        ENDIF
 
@@ -227,22 +236,35 @@
 
 !
-!-- Radiation boundary condition for the velocities at the respective outflow
+!-- Neumann or Radiation boundary condition for the velocities at the 
+!-- respective outflow
     IF ( outflow_s )  THEN
 
-       c_max = dy / dt_3d
-
-       DO i = nxlg, nxrg
+       IF ( bc_ns_dirneu )  THEN
+          u(:,-1,:) = u(:,0,:)
+          v(:,0,:)  = v(:,1,:)
+          w(:,-1,:) = w(:,0,:)          
+       ELSEIF ( bc_ns_dirrad )  THEN
+
+          c_max = dy / dt_3d
+
+          c_u_m_l = 0.0 
+          c_v_m_l = 0.0
+          c_w_m_l = 0.0
+
+          c_u_m = 0.0 
+          c_v_m = 0.0
+          c_w_m = 0.0
+
+!
+!--       Calculate the phase speeds for u,v, and w, first local and then 
+!--       average parallel along the outflow boundary.  
           DO k = nzb+1, nzt+1
-
-!
-!--          Calculate the phase speeds for u,v, and w. In case of using a
-!--          Runge-Kutta scheme, do this for the first substep only and then
-!--          reuse this values for the further substeps.
-             IF ( intermediate_timestep_count == 1 )  THEN
+             DO i = nxl, nxr
 
                 denom = u_m_s(k,0,i) - u_m_s(k,1,i)
 
                 IF ( denom /= 0.0 )  THEN
-                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
+                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) )             &
+                              / ( denom * tsc(2) )
                    IF ( c_u(k,i) < 0.0 )  THEN
                       c_u(k,i) = 0.0
@@ -257,5 +279,6 @@
 
                 IF ( denom /= 0.0 )  THEN
-                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / denom
+                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) )             &
+                              / ( denom * tsc(2) )
                    IF ( c_v(k,i) < 0.0 )  THEN
                       c_v(k,i) = 0.0
@@ -270,5 +293,6 @@
 
                 IF ( denom /= 0.0 )  THEN
-                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
+                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) )             &
+                              / ( denom * tsc(2) )
                    IF ( c_w(k,i) < 0.0 )  THEN
                       c_w(k,i) = 0.0
@@ -280,47 +304,77 @@
                 ENDIF
 
-!
-!--             Save old timelevels for the next timestep
-                u_m_s(k,:,i) = u(k,0:1,i)
-                v_m_s(k,:,i) = v(k,1:2,i)
-                w_m_s(k,:,i) = w(k,0:1,i)
-
-             ENDIF
-
-!
-!--          Calculate the new velocities
-             u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u(k,i) * &
+                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
+                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
+                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
+
+             ENDDO
+          ENDDO
+
+#if defined( __parallel )   
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
+          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dx, ierr )   
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
+          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dx, ierr )  
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
+          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dx, ierr )  
+#else
+          c_u_m = c_u_m_l
+          c_v_m = c_v_m_l
+          c_w_m = c_w_m_l
+#endif
+
+          c_u_m = c_u_m / (nx+1)
+          c_v_m = c_v_m / (nx+1)
+          c_w_m = c_w_m / (nx+1)
+
+!
+!--       Save old timelevels for the next timestep
+          IF ( intermediate_timestep_count == 1 )  THEN
+             u_m_s(:,:,:) = u(:,0:1,:)
+             v_m_s(:,:,:) = v(:,1:2,:)
+             w_m_s(:,:,:) = w(:,0:1,:)
+          ENDIF
+
+!
+!--       Calculate the new velocities
+          DO k = nzb+1, nzt+1 
+             DO i = nxlg, nxrg
+                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *          &
                                        ( u(k,-1,i) - u(k,0,i) ) * ddy
 
-             v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v(k,i) * &
+                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *          &
                                        ( v(k,0,i) - v(k,1,i) ) * ddy
 
-             w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w(k,i) * &
+                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *          &
                                        ( w(k,-1,i) - w(k,0,i) ) * ddy
-
-          ENDDO
-       ENDDO
-
-!
-!--    Bottom boundary at the outflow
-       IF ( ibc_uv_b == 0 )  THEN
-          u_p(nzb,-1,:) = 0.0 
-          v_p(nzb,0,:)  = 0.0  
-       ELSE                    
-          u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
-          v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
-       ENDIF
-       w_p(nzb,-1,:) = 0.0
-
-!
-!--    Top boundary at the outflow
-       IF ( ibc_uv_t == 0 )  THEN
-          u_p(nzt+1,-1,:) = u_init(nzt+1)
-          v_p(nzt+1,0,:)  = v_init(nzt+1)
-       ELSE
-          u_p(nzt+1,-1,:) = u(nzt,-1,:)
-          v_p(nzt+1,0,:)  = v(nzt,0,:)
-       ENDIF
-       w_p(nzt:nzt+1,-1,:) = 0.0
+             ENDDO
+          ENDDO
+
+!
+!--       Bottom boundary at the outflow
+          IF ( ibc_uv_b == 0 )  THEN
+             u_p(nzb,-1,:) = 0.0 
+             v_p(nzb,0,:)  = 0.0  
+          ELSE                    
+             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
+             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
+          ENDIF
+          w_p(nzb,-1,:) = 0.0
+
+!
+!--       Top boundary at the outflow
+          IF ( ibc_uv_t == 0 )  THEN
+             u_p(nzt+1,-1,:) = u_init(nzt+1)
+             v_p(nzt+1,0,:)  = v_init(nzt+1)
+          ELSE
+             u_p(nzt+1,-1,:) = u(nzt,-1,:)
+             v_p(nzt+1,0,:)  = v(nzt,0,:)
+          ENDIF
+          w_p(nzt:nzt+1,-1,:) = 0.0
+
+       ENDIF
 
     ENDIF
@@ -328,19 +382,31 @@
     IF ( outflow_n )  THEN
 
-       c_max = dy / dt_3d
-
-       DO i = nxlg, nxrg
+       IF ( bc_ns_neudir )  THEN
+          u(:,ny+1,:) = u(:,ny,:)
+          v(:,ny+1,:) = v(:,ny,:)
+          w(:,ny+1,:) = w(:,ny,:)          
+       ELSEIF ( bc_ns_dirrad )  THEN
+
+          c_max = dy / dt_3d
+
+          c_u_m_l = 0.0 
+          c_v_m_l = 0.0
+          c_w_m_l = 0.0
+
+          c_u_m = 0.0 
+          c_v_m = 0.0
+          c_w_m = 0.0
+
+!
+!--       Calculate the phase speeds for u,v, and w, first local and then 
+!--       average parallel along the outflow boundary.  
           DO k = nzb+1, nzt+1
-
-!
-!--          Calculate the phase speeds for u,v, and w. In case of using a
-!--          Runge-Kutta scheme, do this for the first substep only and then
-!--          reuse this values for the further substeps.
-             IF ( intermediate_timestep_count == 1 )  THEN
+             DO i = nxl, nxr
 
                 denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
 
                 IF ( denom /= 0.0 )  THEN
-                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
+                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) )           &
+                              / ( denom * tsc(2) )
                    IF ( c_u(k,i) < 0.0 )  THEN
                       c_u(k,i) = 0.0
@@ -355,5 +421,6 @@
 
                 IF ( denom /= 0.0 )  THEN
-                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
+                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) )           &
+                              / ( denom * tsc(2) )
                    IF ( c_v(k,i) < 0.0 )  THEN
                       c_v(k,i) = 0.0
@@ -368,5 +435,6 @@
 
                 IF ( denom /= 0.0 )  THEN
-                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
+                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) )           &
+                              / ( denom * tsc(2) )
                    IF ( c_w(k,i) < 0.0 )  THEN
                       c_w(k,i) = 0.0
@@ -378,47 +446,77 @@
                 ENDIF
 
-!
-!--             Swap timelevels for the next timestep
-                u_m_n(k,:,i) = u(k,ny-1:ny,i)
-                v_m_n(k,:,i) = v(k,ny-1:ny,i)
-                w_m_n(k,:,i) = w(k,ny-1:ny,i)
-
-             ENDIF
-
-!
-!--          Calculate the new velocities
-             u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u(k,i) * &
-                                           ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
-
-             v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v(k,i) * &
-                                           ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
-
-             w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w(k,i) * &
-                                           ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
-
-          ENDDO
-       ENDDO
-
-!
-!--    Bottom boundary at the outflow
-       IF ( ibc_uv_b == 0 )  THEN
-          u_p(nzb,ny+1,:) = 0.0 
-          v_p(nzb,ny+1,:) = 0.0   
-       ELSE                    
-          u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
-          v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
-       ENDIF
-       w_p(nzb,ny+1,:) = 0.0
-
-!
-!--    Top boundary at the outflow
-       IF ( ibc_uv_t == 0 )  THEN
-          u_p(nzt+1,ny+1,:) = u_init(nzt+1)
-          v_p(nzt+1,ny+1,:) = v_init(nzt+1)
-       ELSE
-          u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
-          v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
-       ENDIF
-       w_p(nzt:nzt+1,ny+1,:) = 0.0
+                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
+                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
+                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
+
+             ENDDO
+          ENDDO
+
+#if defined( __parallel )   
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
+          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dx, ierr )   
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
+          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dx, ierr )  
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
+          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dx, ierr )  
+#else
+          c_u_m = c_u_m_l
+          c_v_m = c_v_m_l
+          c_w_m = c_w_m_l
+#endif
+
+          c_u_m = c_u_m / (nx+1)
+          c_v_m = c_v_m / (nx+1)
+          c_w_m = c_w_m / (nx+1)
+
+!
+!--       Save old timelevels for the next timestep
+          IF ( intermediate_timestep_count == 1 )  THEN
+                u_m_n(:,:,:) = u(:,ny-1:ny,:)
+                v_m_n(:,:,:) = v(:,ny-1:ny,:)
+                w_m_n(:,:,:) = w(:,ny-1:ny,:)
+          ENDIF
+
+!
+!--       Calculate the new velocities
+          DO k = nzb+1, nzt+1 
+             DO i = nxlg, nxrg
+                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *      &
+                                       ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
+
+                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *     &
+                                       ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
+
+                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *      &
+                                       ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
+             ENDDO
+          ENDDO
+
+!
+!--       Bottom boundary at the outflow
+          IF ( ibc_uv_b == 0 )  THEN
+             u_p(nzb,ny+1,:) = 0.0 
+             v_p(nzb,ny+1,:) = 0.0   
+          ELSE                    
+             u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
+             v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
+          ENDIF
+          w_p(nzb,ny+1,:) = 0.0
+
+!
+!--       Top boundary at the outflow
+          IF ( ibc_uv_t == 0 )  THEN
+             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
+             v_p(nzt+1,ny+1,:) = v_init(nzt+1)
+          ELSE
+             u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
+             v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
+          ENDIF
+          w_p(nzt:nzt+1,ny+1,:) = 0.0
+
+       ENDIF
 
     ENDIF
@@ -426,19 +524,31 @@
     IF ( outflow_l )  THEN
 
-       c_max = dx / dt_3d
-
-       DO j = nysg, nyng
+       IF ( bc_lr_neudir )  THEN
+          u(:,:,-1) = u(:,:,0)
+          v(:,:,0)  = v(:,:,1)
+          w(:,:,-1) = w(:,:,0)          
+       ELSEIF ( bc_ns_dirrad )  THEN
+
+          c_max = dx / dt_3d
+
+          c_u_m_l = 0.0 
+          c_v_m_l = 0.0
+          c_w_m_l = 0.0
+
+          c_u_m = 0.0 
+          c_v_m = 0.0
+          c_w_m = 0.0
+
+!
+!--       Calculate the phase speeds for u,v, and w, first local and then 
+!--       average parallel along the outflow boundary.  
           DO k = nzb+1, nzt+1
-
-!
-!--          Calculate the phase speeds for u,v, and w. In case of using a
-!--          Runge-Kutta scheme, do this for the first substep only and then
-!--          reuse this values for the further substeps.
-             IF ( intermediate_timestep_count == 1 )  THEN
+             DO j = nys, nyn
 
                 denom = u_m_l(k,j,1) - u_m_l(k,j,2)
 
                 IF ( denom /= 0.0 )  THEN
-                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / denom
+                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) )             &
+                              / ( denom * tsc(2) )
                    IF ( c_u(k,j) < 0.0 )  THEN
                       c_u(k,j) = 0.0
@@ -453,5 +563,6 @@
 
                 IF ( denom /= 0.0 )  THEN
-                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
+                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) )             &
+                              / ( denom * tsc(2) )
                    IF ( c_v(k,j) < 0.0 )  THEN
                       c_v(k,j) = 0.0
@@ -466,5 +577,6 @@
 
                 IF ( denom /= 0.0 )  THEN
-                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
+                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) )             &
+                              / ( denom * tsc(2) )
                    IF ( c_w(k,j) < 0.0 )  THEN
                       c_w(k,j) = 0.0
@@ -476,47 +588,77 @@
                 ENDIF
 
-!
-!--             Swap timelevels for the next timestep
-                u_m_l(k,j,:) = u(k,j,1:2)
-                v_m_l(k,j,:) = v(k,j,0:1)
-                w_m_l(k,j,:) = w(k,j,0:1)
-
-             ENDIF
-
-!
-!--          Calculate the new velocities
-             u_p(k,j,0)  = u(k,j,0)  - dt_3d * tsc(2) * c_u(k,j) * &
+                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
+                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
+                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
+
+             ENDDO
+          ENDDO
+
+#if defined( __parallel )   
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, ierr )   
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, ierr )  
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, ierr )  
+#else
+          c_u_m = c_u_m_l
+          c_v_m = c_v_m_l
+          c_w_m = c_w_m_l
+#endif
+
+          c_u_m = c_u_m / (ny+1)
+          c_v_m = c_v_m / (ny+1)
+          c_w_m = c_w_m / (ny+1)
+
+!
+!--       Save old timelevels for the next timestep
+          IF ( intermediate_timestep_count == 1 )  THEN
+                u_m_l(:,:,:) = u(:,:,1:2)
+                v_m_l(:,:,:) = v(:,:,0:1)
+                w_m_l(:,:,:) = w(:,:,0:1)
+          ENDIF
+
+!
+!--       Calculate the new velocities
+          DO k = nzb+1, nzt+1 
+             DO i = nxlg, nxrg
+                u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
                                        ( u(k,j,0) - u(k,j,1) ) * ddx
 
-             v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v(k,j) * &
+                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *          &
                                        ( v(k,j,-1) - v(k,j,0) ) * ddx
 
-             w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w(k,j) * &
+                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *          &
                                        ( w(k,j,-1) - w(k,j,0) ) * ddx
-
-          ENDDO
-       ENDDO
-
-!
-!--    Bottom boundary at the outflow
-       IF ( ibc_uv_b == 0 )  THEN
-          u_p(nzb,:,0)  = 0.0 
-          v_p(nzb,:,-1) = 0.0
-       ELSE                    
-          u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
-          v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
-       ENDIF
-       w_p(nzb,:,-1) = 0.0
-
-!
-!--    Top boundary at the outflow
-       IF ( ibc_uv_t == 0 )  THEN
-          u_p(nzt+1,:,-1) = u_init(nzt+1)
-          v_p(nzt+1,:,-1) = v_init(nzt+1)
-       ELSE
-          u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
-          v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
-       ENDIF
-       w_p(nzt:nzt+1,:,-1) = 0.0
+             ENDDO
+          ENDDO
+
+!
+!--       Bottom boundary at the outflow
+          IF ( ibc_uv_b == 0 )  THEN
+             u_p(nzb,:,0)  = 0.0 
+             v_p(nzb,:,-1) = 0.0
+          ELSE                    
+             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
+             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
+          ENDIF
+          w_p(nzb,:,-1) = 0.0
+
+!
+!--       Top boundary at the outflow
+          IF ( ibc_uv_t == 0 )  THEN
+             u_p(nzt+1,:,-1) = u_init(nzt+1)
+             v_p(nzt+1,:,-1) = v_init(nzt+1)
+          ELSE
+             u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
+             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
+          ENDIF
+          w_p(nzt:nzt+1,:,-1) = 0.0
+
+       ENDIF
 
     ENDIF
@@ -524,19 +666,31 @@
     IF ( outflow_r )  THEN
 
-       c_max = dx / dt_3d
-
-       DO j = nysg, nyng
+       IF ( bc_lr_dirneu )  THEN
+          u(:,:,nx+1) = u(:,:,nx)
+          v(:,:,nx+1) = v(:,:,nx)
+          w(:,:,nx+1) = w(:,:,nx)          
+       ELSEIF ( bc_ns_dirrad )  THEN
+
+          c_max = dx / dt_3d
+
+          c_u_m_l = 0.0 
+          c_v_m_l = 0.0
+          c_w_m_l = 0.0
+
+          c_u_m = 0.0 
+          c_v_m = 0.0
+          c_w_m = 0.0
+
+!
+!--       Calculate the phase speeds for u,v, and w, first local and then 
+!--       average parallel along the outflow boundary.  
           DO k = nzb+1, nzt+1
-
-!
-!--          Calculate the phase speeds for u,v, and w. In case of using a
-!--          Runge-Kutta scheme, do this for the first substep only and then
-!--          reuse this values for the further substeps.
-             IF ( intermediate_timestep_count == 1 )  THEN
+             DO j = nys, nyn
 
                 denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
 
                 IF ( denom /= 0.0 )  THEN
-                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
+                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) )           &
+                              / ( denom * tsc(2) )
                    IF ( c_u(k,j) < 0.0 )  THEN
                       c_u(k,j) = 0.0
@@ -551,5 +705,6 @@
 
                 IF ( denom /= 0.0 )  THEN
-                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
+                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) )           &
+                              / ( denom * tsc(2) )
                    IF ( c_v(k,j) < 0.0 )  THEN
                       c_v(k,j) = 0.0
@@ -564,5 +719,6 @@
 
                 IF ( denom /= 0.0 )  THEN
-                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
+                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) )           &
+                              / ( denom * tsc(2) )
                    IF ( c_w(k,j) < 0.0 )  THEN
                       c_w(k,j) = 0.0
@@ -574,50 +730,79 @@
                 ENDIF
 
-!
-!--             Swap timelevels for the next timestep
-                u_m_r(k,j,:) = u(k,j,nx-1:nx)
-                v_m_r(k,j,:) = v(k,j,nx-1:nx)
-                w_m_r(k,j,:) = w(k,j,nx-1:nx)
-
-             ENDIF
-
-!
-!--          Calculate the new velocities
-             u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u(k,j) * &
-                                           ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
-
-             v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v(k,j) * &
-                                           ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
-
-             w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w(k,j) * &
-                                           ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
-
-          ENDDO
-       ENDDO
-
-!
-!--    Bottom boundary at the outflow
-       IF ( ibc_uv_b == 0 )  THEN
-          u_p(nzb,:,nx+1) = 0.0
-          v_p(nzb,:,nx+1) = 0.0 
-       ELSE                    
-          u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
-          v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
-       ENDIF
-       w_p(nzb,:,nx+1) = 0.0
-
-!
-!--    Top boundary at the outflow
-       IF ( ibc_uv_t == 0 )  THEN
-          u_p(nzt+1,:,nx+1) = u_init(nzt+1)
-          v_p(nzt+1,:,nx+1) = v_init(nzt+1)
-       ELSE
-          u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
-          v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
-       ENDIF
-       w(nzt:nzt+1,:,nx+1) = 0.0
+                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
+                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
+                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
+
+             ENDDO
+          ENDDO
+
+#if defined( __parallel )   
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, ierr )   
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, ierr )  
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, ierr )  
+#else
+          c_u_m = c_u_m_l
+          c_v_m = c_v_m_l
+          c_w_m = c_w_m_l
+#endif
+
+          c_u_m = c_u_m / (ny+1)
+          c_v_m = c_v_m / (ny+1)
+          c_w_m = c_w_m / (ny+1)
+
+!
+!--       Save old timelevels for the next timestep
+          IF ( intermediate_timestep_count == 1 )  THEN
+                u_m_r(:,:,:) = u(:,:,nx-1:nx)
+                v_m_r(:,:,:) = v(:,:,nx-1:nx)
+                w_m_r(:,:,:) = w(:,:,nx-1:nx)
+          ENDIF
+
+!
+!--       Calculate the new velocities
+          DO k = nzb+1, nzt+1 
+             DO i = nxlg, nxrg
+                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
+                                       ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
+
+                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *      &
+                                       ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
+
+                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *      &
+                                       ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
+             ENDDO
+          ENDDO
+
+!
+!--       Bottom boundary at the outflow
+          IF ( ibc_uv_b == 0 )  THEN
+             u_p(nzb,:,nx+1) = 0.0
+             v_p(nzb,:,nx+1) = 0.0 
+          ELSE                    
+             u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
+             v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
+          ENDIF
+          w_p(nzb,:,nx+1) = 0.0
+
+!
+!--       Top boundary at the outflow
+          IF ( ibc_uv_t == 0 )  THEN
+             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
+             v_p(nzt+1,:,nx+1) = v_init(nzt+1)
+          ELSE
+             u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
+             v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
+          ENDIF
+          w(nzt:nzt+1,:,nx+1) = 0.0
+
+       ENDIF
 
     ENDIF
 
- 
  END SUBROUTINE boundary_conds
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 977)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 978)
@@ -4,5 +4,8 @@
 ! Current revisions:
 ! -----------------
-!
+! setting of bc_lr/ns_dirneu/neudir
+! outflow damping layer removed
+! check for z0h*
+! check for pt_damping_width
 !
 ! Former revisions:
@@ -1379,5 +1382,6 @@
 !-- Lateral boundary conditions
     IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/radiation'  .AND. &
-         bc_lr /= 'radiation/dirichlet' )  THEN
+         bc_lr /= 'radiation/dirichlet' .AND. bc_lr /= 'dirichlet/neumann' &
+         .AND. bc_lr /= 'neumann/dirichlet' )  THEN
        message_string = 'unknown boundary condition: bc_lr = "' // &
                         TRIM( bc_lr ) // '"'
@@ -1385,5 +1389,6 @@
     ENDIF
     IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/radiation'  .AND. &
-         bc_ns /= 'radiation/dirichlet' )  THEN
+         bc_ns /= 'radiation/dirichlet' .AND. bc_ns /= 'dirichlet/neumann' &
+         .AND. bc_ns /= 'neumann/dirichlet' )  THEN
        message_string = 'unknown boundary condition: bc_ns = "' // &
                         TRIM( bc_ns ) // '"'
@@ -1396,7 +1401,11 @@
     IF ( bc_lr == 'dirichlet/radiation' )  bc_lr_dirrad = .TRUE.
     IF ( bc_lr == 'radiation/dirichlet' )  bc_lr_raddir = .TRUE.
+    IF ( bc_lr == 'dirichlet/neumann' )    bc_lr_dirneu = .TRUE.
+    IF ( bc_lr == 'neumann/dirichlet' )    bc_lr_neudir = .TRUE.
     IF ( bc_ns /= 'cyclic' )               bc_ns_cyc    = .FALSE.
     IF ( bc_ns == 'dirichlet/radiation' )  bc_ns_dirrad = .TRUE.
     IF ( bc_ns == 'radiation/dirichlet' )  bc_ns_raddir = .TRUE.
+    IF ( bc_ns == 'dirichlet/neumann' )    bc_ns_dirneu = .TRUE.
+    IF ( bc_ns == 'neumann/dirichlet' )    bc_ns_neudir = .TRUE.
 
 !
@@ -2659,5 +2668,5 @@
              unit = 'psu'
 
-          CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'qsws*', 'shf*', 'z0*' )
+          CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'qsws*', 'shf*', 'z0*', 'z0h*' )
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
                 message_string = 'illegal value for data_output: "' // &
@@ -2700,4 +2709,5 @@
              IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
              IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
+             IF ( TRIM( var ) == 'z0h*'    )  unit = 'm'
 
 
@@ -2964,16 +2974,9 @@
 
 !
-!-- In case of non-cyclic lateral boundaries, set the default maximum value
-!-- for the horizontal diffusivity used within the outflow damping layer,
-!-- and check/set the width of the damping layer
+!-- In case of non-cyclic lateral boundaries and a damping layer for the
+!-- potential temperature, check the width of the damping layer
     IF ( bc_lr /= 'cyclic' ) THEN
-       IF ( km_damp_max == -1.0 )  THEN
-          km_damp_max = 0.5 * dx
-       ENDIF
-       IF ( outflow_damping_width == -1.0 )  THEN
-          outflow_damping_width = MIN( 20, nx/2 )
-       ENDIF
-       IF ( outflow_damping_width <= 0  .OR.  outflow_damping_width > nx )  THEN
-          message_string = 'outflow_damping width out of range'
+       IF ( pt_damping_width < 0.0 .OR. pt_damping_width > REAL( nx * dx ) )  THEN
+          message_string = 'pt_damping_width out of range'
           CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
        ENDIF
@@ -2981,12 +2984,6 @@
 
     IF ( bc_ns /= 'cyclic' )  THEN
-       IF ( km_damp_max == -1.0 )  THEN
-          km_damp_max = 0.5 * dy
-       ENDIF
-       IF ( outflow_damping_width == -1.0 )  THEN
-          outflow_damping_width = MIN( 20, ny/2 )
-       ENDIF
-       IF ( outflow_damping_width <= 0  .OR.  outflow_damping_width > ny )  THEN
-          message_string = 'outflow_damping width out of range'
+       IF ( pt_damping_width < 0.0 .OR. pt_damping_width > REAL( ny * dy ) )  THEN
+          message_string = 'pt_damping_width out of range'
           CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
        ENDIF
@@ -3108,15 +3105,15 @@
     ENDIF
 
-    IF ( bc_lr == 'radiation/dirichlet' )  THEN
+    IF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' )  THEN
        dist_nxr    = nx - inflow_disturbance_begin
        dist_nxl(1) = nx - inflow_disturbance_end
-    ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
+    ELSEIF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' )  THEN
        dist_nxl    = inflow_disturbance_begin
        dist_nxr(1) = inflow_disturbance_end
     ENDIF
-    IF ( bc_ns == 'dirichlet/radiation' )  THEN
+    IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' )  THEN
        dist_nyn    = ny - inflow_disturbance_begin
        dist_nys(1) = ny - inflow_disturbance_end
-    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
+    ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' )  THEN
        dist_nys    = inflow_disturbance_begin
        dist_nyn(1) = inflow_disturbance_end
@@ -3126,5 +3123,5 @@
 !-- A turbulent inflow requires Dirichlet conditions at the respective inflow
 !-- boundary (so far, a turbulent inflow is realized from the left side only)
-    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
+    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' .AND.  bc_lr /= 'dirichlet/neumann' )  THEN
        message_string = 'turbulent_inflow = .T. requires a Dirichlet ' // &
                         'condition at the inflow boundary'
Index: /palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_2d.f90	(revision 977)
+++ /palm/trunk/SOURCE/data_output_2d.f90	(revision 978)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! +z0h
 !
 ! Former revisions:
@@ -594,4 +594,22 @@
                       DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) =  z0_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'z0h*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  z0h(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  z0h_av(j,i)
                       ENDDO
                    ENDDO
Index: /palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_u.f90	(revision 977)
+++ /palm/trunk/SOURCE/diffusion_u.f90	(revision 978)
@@ -4,4 +4,6 @@
 ! Current revisions:
 ! -----------------
+! outflow damping layer removed
+! kmym_x/_y and kmyp_x/_y change to kmym and kmyp
 !
 ! Former revisions:
@@ -63,6 +65,5 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, &
-                            v, w )
+    SUBROUTINE diffusion_u( ddzu, ddzw, km, tend, u, usws, uswst, v, w )
 
        USE control_parameters
@@ -73,6 +74,6 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nysg:nyng)
+       REAL    ::  kmym, kmyp, kmzm, kmzp
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
        REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
@@ -95,19 +96,8 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmyp_x = 0.25 * &
-                         ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-                kmym_x = 0.25 * &
-                         ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
-                kmyp_y = kmyp_x
-                kmym_y = kmym_x
-!
-!--             Increase diffusion at the outflow boundary in case of
-!--             non-cyclic lateral boundaries. Damping is only needed for
-!--             velocity components parallel to the outflow boundary in
-!--             the direction normal to the outflow boundary.
-                IF ( .NOT. bc_ns_cyc )  THEN
-                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
-                   kmym_y = MAX( kmym_y, km_damp_y(j) )
-                ENDIF
+                kmyp = 0.25 * &
+                       ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
+                kmym = 0.25 * &
+                       ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
                 tend(k,j,i) = tend(k,j,i)                                    &
@@ -116,8 +106,8 @@
                       &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
                       &         ) * ddx2                                     &
-                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
-                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
-                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
-                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
+                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy       &
+                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx       &
+                      &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
+                      &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
                       &   ) * ddy
              ENDDO
@@ -128,19 +118,8 @@
 
                 DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
-                   kmyp_x = 0.25 * &
-                            ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-                   kmym_x = 0.25 * &
-                            ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
-                   kmyp_y = kmyp_x
-                   kmym_y = kmym_x
-!
-!--                Increase diffusion at the outflow boundary in case of
-!--                non-cyclic lateral boundaries. Damping is only needed for
-!--                velocity components parallel to the outflow boundary in
-!--                the direction normal to the outflow boundary.
-                   IF ( .NOT. bc_ns_cyc )  THEN
-                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
-                      kmym_y = MAX( kmym_y, km_damp_y(j) )
-                   ENDIF
+                   kmyp = 0.25 * &
+                          ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
+                   kmym = 0.25 * &
+                          ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
                    tend(k,j,i) = tend(k,j,i)                                   &
@@ -150,10 +129,10 @@
                                          ) * ddx2                              &
                                  + (   fyp(j,i) * (                            &
-                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
-                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
+                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
+                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
                                                   )                            &
                                      - fym(j,i) * (                            &
-                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
-                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
+                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
+                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
                                                   )                            &
                                      + wall_u(j,i) * usvs(k,j,i)               &
@@ -239,5 +218,5 @@
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
+    SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, tend, u, usws, &
                                uswst, v, w )
 
@@ -249,6 +228,6 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nysg:nyng)
+       REAL    ::  kmym, kmyp, kmzm, kmzp
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
        REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(nzb:nzt+1)      ::  usvs
@@ -261,18 +240,6 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-          kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
-          kmyp_y = kmyp_x
-          kmym_y = kmym_x
-
-!
-!--       Increase diffusion at the outflow boundary in case of non-cyclic
-!--       lateral boundaries. Damping is only needed for velocity components
-!--       parallel to the outflow boundary in the direction normal to the
-!--       outflow boundary.
-          IF ( .NOT. bc_ns_cyc )  THEN
-             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
-             kmym_y = MAX( kmym_y, km_damp_y(j) )
-          ENDIF
+          kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
+          kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
           tend(k,j,i) = tend(k,j,i)                                          &
@@ -281,8 +248,8 @@
                       &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
                       &         ) * ddx2                                     &
-                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
-                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
-                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
-                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
+                      & + ( kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy       &
+                      &   + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx       &
+                      &   - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
+                      &   - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
                       &   ) * ddy
        ENDDO
@@ -298,17 +265,6 @@
 
           DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
-             kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
-             kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
-             kmyp_y = kmyp_x
-             kmym_y = kmym_x
-!
-!--          Increase diffusion at the outflow boundary in case of
-!--          non-cyclic lateral boundaries. Damping is only needed for
-!--          velocity components parallel to the outflow boundary in
-!--          the direction normal to the outflow boundary.
-             IF ( .NOT. bc_ns_cyc )  THEN
-                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
-                kmym_y = MAX( kmym_y, km_damp_y(j) )
-             ENDIF
+             kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
+             kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
 
              tend(k,j,i) = tend(k,j,i)                                         &
@@ -318,10 +274,10 @@
                                          ) * ddx2                              &
                                  + (   fyp(j,i) * (                            &
-                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
-                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
+                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy   &
+                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx   &
                                                   )                            &
                                      - fym(j,i) * (                            &
-                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
-                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
+                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
+                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
                                                   )                            &
                                      + wall_u(j,i) * usvs(k)                   &
Index: /palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_v.f90	(revision 977)
+++ /palm/trunk/SOURCE/diffusion_v.f90	(revision 978)
@@ -4,4 +4,6 @@
 ! Current revisions:
 ! -----------------
+! outflow damping layer removed
+! kmxm_x/_y and kmxp_x/_y change to kmxm and kmxp
 !
 ! Former revisions:
@@ -61,6 +63,5 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, &
-                            vswst, w )
+    SUBROUTINE diffusion_v( ddzu, ddzw, km, tend, u, v, vsws, vswst, w )
 
        USE control_parameters
@@ -71,6 +72,6 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg)
+       REAL    ::  kmxm, kmxp, kmzm, kmzp
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
        REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst
@@ -93,25 +94,14 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmxp_x = 0.25 * &
-                         ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
-                kmxm_x = 0.25 * &
-                         ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
-                kmxp_y = kmxp_x
-                kmxm_y = kmxm_x
-!
-!--             Increase diffusion at the outflow boundary in case of
-!--             non-cyclic lateral boundaries. Damping is only needed for
-!--             velocity components parallel to the outflow boundary in
-!--             the direction normal to the outflow boundary.
-                IF ( .NOT. bc_lr_cyc )  THEN
-                   kmxp_x = MAX( kmxp_x, km_damp_x(i) )
-                   kmxm_x = MAX( kmxm_x, km_damp_x(i) )
-                ENDIF
+                kmxp = 0.25 * &
+                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
+                kmxm = 0.25 * &
+                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
 
                 tend(k,j,i) = tend(k,j,i)                                    &
-                      & + ( kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx     &
-                      &   + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy     &
-                      &   - kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
-                      &   - kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
+                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx       &
+                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy       &
+                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
+                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
                       &   ) * ddx                                            &
                       & + 2.0 * (                                            &
@@ -126,20 +116,9 @@
 
                 DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
-                   kmxp_x = 0.25 * &
-                            ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
-                   kmxm_x = 0.25 * &
-                            ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
-                   kmxp_y = kmxp_x
-                   kmxm_y = kmxm_x
-!
-!--                Increase diffusion at the outflow boundary in case of
-!--                non-cyclic lateral boundaries. Damping is only needed for
-!--                velocity components parallel to the outflow boundary in
-!--                the direction normal to the outflow boundary.
-                   IF ( .NOT. bc_lr_cyc )  THEN
-                      kmxp_x = MAX( kmxp_x, km_damp_x(i) )
-                      kmxm_x = MAX( kmxm_x, km_damp_x(i) )
-                   ENDIF
-
+                   kmxp = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
+                   kmxm = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
+                   
                    tend(k,j,i) = tend(k,j,i)                                   &
                                  + 2.0 * (                                     &
@@ -148,10 +127,10 @@
                                          ) * ddy2                              &
                                  + (   fxp(j,i) * (                            &
-                                  kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx &
-                                + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
+                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
+                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
                                                   )                            &
                                      - fxm(j,i) * (                            &
-                                  kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
-                                + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
+                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
+                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
                                                   )                            &
                                      + wall_v(j,i) * vsus(k,j,i)               &
@@ -237,5 +216,5 @@
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
+    SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, tend, u, v, &
                                vsws, vswst, w )
 
@@ -247,6 +226,6 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg)
+       REAL    ::  kmxm, kmxp, kmzm, kmzp
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
        REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(nzb:nzt+1)      ::  vsus
@@ -259,23 +238,12 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
-          kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
-          kmxp_y = kmxp_x
-          kmxm_y = kmxm_x
-!
-!--       Increase diffusion at the outflow boundary in case of non-cyclic
-!--       lateral boundaries. Damping is only needed for velocity components
-!--       parallel to the outflow boundary in the direction normal to the
-!--       outflow boundary.
-          IF ( .NOT. bc_lr_cyc )  THEN
-             kmxp_x = MAX( kmxp_x, km_damp_x(i) )
-             kmxm_x = MAX( kmxm_x, km_damp_x(i) )
-          ENDIF
+          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
+          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
 
           tend(k,j,i) = tend(k,j,i)                                          &
-                      & + ( kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx     &
-                      &   + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy     &
-                      &   - kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
-                      &   - kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
+                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx       &
+                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy       &
+                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx           &
+                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy           &
                       &   ) * ddx                                            &
                       & + 2.0 * (                                            &
@@ -295,19 +263,8 @@
 
           DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
-             kmxp_x = 0.25 * &
-                      ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
-             kmxm_x = 0.25 * &
-                      ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
-             kmxp_y = kmxp_x
-             kmxm_y = kmxm_x
-!
-!--          Increase diffusion at the outflow boundary in case of
-!--          non-cyclic lateral boundaries. Damping is only needed for
-!--          velocity components parallel to the outflow boundary in
-!--          the direction normal to the outflow boundary.
-             IF ( .NOT. bc_lr_cyc )  THEN
-                kmxp_x = MAX( kmxp_x, km_damp_x(i) )
-                kmxm_x = MAX( kmxm_x, km_damp_x(i) )
-             ENDIF
+             kmxp = 0.25 * &
+                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
+             kmxm = 0.25 * &
+                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
 
              tend(k,j,i) = tend(k,j,i)                                         &
@@ -317,10 +274,10 @@
                                          ) * ddy2                              &
                                  + (   fxp(j,i) * (                            &
-                                  kmxp_x * ( v(k,j,i+1) - v(k,j,i)     ) * ddx &
-                                + kmxp_y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
+                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
+                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
                                                   )                            &
                                      - fxm(j,i) * (                            &
-                                  kmxm_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
-                                + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
+                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
+                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
                                                   )                            &
                                      + wall_v(j,i) * vsus(k)                   &
Index: /palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_w.f90	(revision 977)
+++ /palm/trunk/SOURCE/diffusion_w.f90	(revision 978)
@@ -4,4 +4,7 @@
 ! Current revisions:
 ! -----------------
+! outflow damping layer removed
+! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
+! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
 !
 ! Former revisions:
@@ -54,6 +57,5 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
-                            w )
+    SUBROUTINE diffusion_w( ddzu, ddzw, km, tend, u, v, w )
 
        USE control_parameters
@@ -64,8 +66,6 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
-                   kmyp_z
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg),        & 
-                   km_damp_y(nysg:nyng)
+       REAL    ::  kmxm, kmxp, kmym, kmyp
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
        REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
@@ -88,40 +88,23 @@
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
-                kmxp_x = 0.25 * &
-                         ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-                kmxm_x = 0.25 * &
-                         ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-                kmxp_z = kmxp_x
-                kmxm_z = kmxm_x
-                kmyp_y = 0.25 * &
-                         ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-                kmym_y = 0.25 * &
-                         ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-                kmyp_z = kmyp_y
-                kmym_z = kmym_y
-!
-!--             Increase diffusion at the outflow boundary in case of
-!--             non-cyclic lateral boundaries. Damping is only needed for
-!--             velocity components parallel to the outflow boundary in
-!--             the direction normal to the outflow boundary.
-                IF ( .NOT. bc_lr_cyc )  THEN
-                   kmxp_x = MAX( kmxp_x, km_damp_x(i) )
-                   kmxm_x = MAX( kmxm_x, km_damp_x(i) )
-                ENDIF
-                IF ( .NOT. bc_ns_cyc )  THEN
-                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
-                   kmym_y = MAX( kmym_y, km_damp_y(j) )
-                ENDIF
+                kmxp = 0.25 * &
+                       ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+                kmxm = 0.25 * &
+                       ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+                kmyp = 0.25 * &
+                       ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+                kmym = 0.25 * &
+                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
 
                 tend(k,j,i) = tend(k,j,i)                                      &
-                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
-                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
-                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
-                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
+                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
+                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
+                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
+                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
                       &   ) * ddx                                              &
-                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
-                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
-                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
-                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
+                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
+                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
+                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
+                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
                       &   ) * ddy                                              &
                       & + 2.0 * (                                              &
@@ -138,48 +121,31 @@
 !
 !--                Interpolate eddy diffusivities on staggered gridpoints
-                   kmxp_x = 0.25 * &
-                            ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-                   kmxm_x = 0.25 * &
-                            ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-                   kmxp_z = kmxp_x
-                   kmxm_z = kmxm_x
-                   kmyp_y = 0.25 * &
-                            ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-                   kmym_y = 0.25 * &
-                            ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-                   kmyp_z = kmyp_y
-                   kmym_z = kmym_y
-!
-!--                Increase diffusion at the outflow boundary in case of
-!--                non-cyclic lateral boundaries. Damping is only needed for
-!--                velocity components parallel to the outflow boundary in
-!--                the direction normal to the outflow boundary.
-                   IF ( .NOT. bc_lr_cyc )  THEN
-                      kmxp_x = MAX( kmxp_x, km_damp_x(i) )
-                      kmxm_x = MAX( kmxm_x, km_damp_x(i) )
-                   ENDIF
-                   IF ( .NOT. bc_ns_cyc )  THEN
-                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
-                      kmym_y = MAX( kmym_y, km_damp_y(j) )
-                   ENDIF
+                   kmxp = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+                   kmxm = 0.25 * &
+                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+                   kmyp = 0.25 * &
+                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+                   kmym = 0.25 * &
+                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
 
                    tend(k,j,i) = tend(k,j,i)                                   &
                                  + (   fwxp(j,i) * (                           &
-                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
-                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
+                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
+                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
                                                    )                           &
                                      - fwxm(j,i) * (                           &
-                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
-                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
+                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
+                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
                                                    )                           &
                                      + wall_w_x(j,i) * wsus(k,j,i)             &
                                    ) * ddx                                     &
                                  + (   fwyp(j,i) * (                           &
-                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
-                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
+                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
+                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
                                                    )                           &
                                      - fwym(j,i) * (                           &
-                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
-                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
+                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
+                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
                                                    )                           &
                                      + wall_w_y(j,i) * wsvs(k,j,i)             &
@@ -201,6 +167,5 @@
 ! Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
-                               tend, u, v, w )
+    SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, tend, u, v, w )
 
        USE control_parameters
@@ -211,8 +176,6 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
-                   kmyp_z
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg),        & 
-                   km_damp_y(nysg:nyng)
+       REAL    ::  kmxm, kmxp, kmym, kmyp
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
        REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
@@ -223,36 +186,19 @@
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
-          kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-          kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-          kmxp_z = kmxp_x
-          kmxm_z = kmxm_x
-          kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-          kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-          kmyp_z = kmyp_y
-          kmym_z = kmym_y
-!
-!--       Increase diffusion at the outflow boundary in case of non-cyclic
-!--       lateral boundaries. Damping is only needed for velocity components
-!--       parallel to the outflow boundary in the direction normal to the
-!--       outflow boundary.
-          IF ( .NOT. bc_lr_cyc )  THEN
-             kmxp_x = MAX( kmxp_x, km_damp_x(i) )
-             kmxm_x = MAX( kmxm_x, km_damp_x(i) )
-          ENDIF
-          IF ( .NOT. bc_ns_cyc )  THEN
-             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
-             kmym_y = MAX( kmym_y, km_damp_y(j) )
-          ENDIF
+          kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+          kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+          kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+          kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
 
           tend(k,j,i) = tend(k,j,i)                                            &
-                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
-                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
-                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
-                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
+                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
+                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
+                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
+                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
                       &   ) * ddx                                              &
-                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
-                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
-                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
-                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
+                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
+                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
+                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
+                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
                       &   ) * ddy                                              &
                       & + 2.0 * (                                              &
@@ -285,44 +231,27 @@
 !
 !--          Interpolate eddy diffusivities on staggered gridpoints
-             kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
-             kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
-             kmxp_z = kmxp_x
-             kmxm_z = kmxm_x
-             kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
-             kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
-             kmyp_z = kmyp_y
-             kmym_z = kmym_y
-!
-!--          Increase diffusion at the outflow boundary in case of
-!--          non-cyclic lateral boundaries. Damping is only needed for
-!--          velocity components parallel to the outflow boundary in
-!--          the direction normal to the outflow boundary.
-             IF ( .NOT. bc_lr_cyc )  THEN
-                kmxp_x = MAX( kmxp_x, km_damp_x(i) )
-                kmxm_x = MAX( kmxm_x, km_damp_x(i) )
-             ENDIF
-             IF ( .NOT. bc_ns_cyc )  THEN
-                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
-                kmym_y = MAX( kmym_y, km_damp_y(j) )
-             ENDIF
+             kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
+             kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
+             kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
+             kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
 
              tend(k,j,i) = tend(k,j,i)                                         &
                                  + (   fwxp(j,i) * (                           &
-                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
-                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
+                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
+                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
                                                    )                           &
                                      - fwxm(j,i) * (                           &
-                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
-                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
+                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
+                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
                                                    )                           &
                                      + wall_w_x(j,i) * wsus(k)                 &
                                    ) * ddx                                     &
                                  + (   fwyp(j,i) * (                           &
-                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
-                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
+                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
+                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
                                                    )                           &
                                      - fwym(j,i) * (                           &
-                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
-                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
+                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
+                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
                                                    )                           &
                                      + wall_w_y(j,i) * wsvs(k)                 &
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 977)
+++ /palm/trunk/SOURCE/header.f90	(revision 978)
@@ -4,5 +4,7 @@
 ! Current revisions:
 ! -----------------
-!
+! -km_damp_max, outflow_damping_width
+! +pt_damping_factor, pt_damping_width
+! +z0h
 !
 ! Former revisions:
@@ -704,5 +706,6 @@
 
     IF ( prandtl_layer )  THEN
-       WRITE ( io, 305 )  (zu(1)-zu(0)), roughness_length, kappa, &
+       WRITE ( io, 305 )  (zu(1)-zu(0)), roughness_length, &
+                          z0h_factor*roughness_length, kappa, &
                           rif_min, rif_max
        IF ( .NOT. constant_heatflux )  WRITE ( io, 308 )
@@ -721,5 +724,5 @@
     WRITE ( io, 317 )  bc_lr, bc_ns
     IF ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )  THEN
-       WRITE ( io, 318 )  outflow_damping_width, km_damp_max
+       WRITE ( io, 318 )  pt_damping_width, pt_damping_factor       
        IF ( turbulent_inflow )  THEN
           WRITE ( io, 319 )  recycling_width, recycling_plane, &
@@ -1779,5 +1782,6 @@
 305 FORMAT (//'    Prandtl-Layer between bottom surface and first ', &
                'computational u,v-level:'// &
-             '       zp = ',F6.2,' m   z0 = ',F6.4,' m   kappa = ',F4.2/ &
+             '       zp = ',F6.2,' m   z0 = ',F6.4,' m   z0h = ',F7.5,&
+             ' m   kappa = ',F4.2/ &
              '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
 306 FORMAT ('       Predefined constant heatflux:   ',F9.6,' K m/s')
@@ -1797,6 +1801,6 @@
             '       left/right:  ',A/    &
             '       north/south: ',A)
-318 FORMAT (/'       outflow damping layer width: ',I3,' gridpoints with km_', &
-                    'max =',F5.1,' m**2/s')
+318 FORMAT (/'       pt damping layer width = ',F7.2,' m, pt ', &
+                    'damping factor = ',F6.4)
 319 FORMAT ('       turbulence recycling at inflow switched on'/ &
             '       width of recycling domain: ',F7.1,' m   grid index: ',I4/ &
Index: /palm/trunk/SOURCE/init_1d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_1d_model.f90	(revision 977)
+++ /palm/trunk/SOURCE/init_1d_model.f90	(revision 978)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! roughness length for scalar quantities z0h1d added
 !
 ! Former revisions:
@@ -149,4 +149,5 @@
     vsws1d = 0.0; vsws1d_m = 0.0
     z01d = roughness_length
+    z0h1d = z0h_factor * z01d 
     IF ( humidity .OR. passive_scalar )  qs1d = 0.0
 
@@ -445,6 +446,6 @@
 !--                Stable stratification
                    ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /      &
-                          ( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
-                                          ( zu(nzb+1) - z01d ) / zu(nzb+1) &
+                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * &
+                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
                           )
                 ELSE
@@ -452,5 +453,5 @@
 !--                Unstable stratification
                    a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) )
-                   b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z01d )
+                   b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0h1d )
 !
 !--                In the borderline case the formula for stable stratification 
@@ -459,6 +460,6 @@
                    IF ( a == 0.0  .OR.  b == 0.0 )  THEN
                       ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /      &
-                             ( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
-                                             ( zu(nzb+1) - z01d ) / zu(nzb+1) &
+                             ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * &
+                                             ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
                              )
                    ELSE
@@ -578,6 +579,6 @@
 !--                Stable stratification
                    qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /        &
-                          ( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
-                                          ( zu(nzb+1) - z01d ) / zu(nzb+1) &
+                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * &
+                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
                           )
                 ELSE
@@ -585,5 +586,5 @@
 !--                Unstable stratification
                    a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) )
-                   b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z01d )
+                   b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0h1d )
 !
 !--                In the borderline case the formula for stable stratification 
@@ -592,6 +593,6 @@
                    IF ( a == 1.0  .OR.  b == 1.0 )  THEN
                       qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /        &
-                             ( LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * &
-                                             ( zu(nzb+1) - z01d ) / zu(nzb+1) &
+                             ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * &
+                                             ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &
                              )
                    ELSE
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 977)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 978)
@@ -7,5 +7,10 @@
 ! Current revisions:
 ! ------------------
-! 
+! outflow damping layer removed
+! roughness length for scalar quantites z0h added
+! damping zone for the potential temperatur in case of non-cyclic lateral
+! boundaries added
+! initialization of ptdf_x, ptdf_y
+! initialization of c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l
 !
 ! Former revisions:
@@ -164,4 +169,5 @@
     USE control_parameters
     USE cpulog
+    USE grid_variables
     USE indices
     USE interfaces
@@ -211,5 +217,5 @@
               sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                &
               ts_value(dots_max,0:statistic_regions) )
-    ALLOCATE( km_damp_x(nxlg:nxrg), km_damp_y(nysg:nyng) )
+    ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
 
     ALLOCATE( rif_1(nysg:nyng,nxlg:nxrg), shf_1(nysg:nyng,nxlg:nxrg), &
@@ -218,5 +224,6 @@
               uswst_1(nysg:nyng,nxlg:nxrg),                               &
               vsws_1(nysg:nyng,nxlg:nxrg),                                &
-              vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg) )
+              vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg),      &
+              z0h(nysg:nyng,nxlg:nxrg) )
 
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
@@ -413,4 +420,9 @@
                  c_w(nzb:nzt+1,nxlg:nxrg) )
     ENDIF
+    IF ( outflow_l  .OR.  outflow_r .OR. outflow_s .OR. outflow_n )  THEN
+       ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) )                   
+       ALLOCATE( c_u_m(nzb:nzt+1), c_v_m(nzb:nzt+1), c_w_m(nzb:nzt+1) )
+    ENDIF
+
 
 !
@@ -847,4 +859,5 @@
 
           z0 = roughness_length
+          z0h = z0h_factor * z0
 
           IF ( .NOT. constant_heatflux )  THEN 
@@ -1533,61 +1546,41 @@
 
 !
-!-- Initialize diffusivities used within the outflow damping layer in case of
-!-- non-cyclic lateral boundaries. A linear increase is assumed over the first
-!-- half of the width of the damping layer
-    IF ( bc_lr_dirrad )  THEN
-
-       DO  i = nxlg, nxrg
-          IF ( i >= nx - outflow_damping_width )  THEN
-             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
-                            ( i - ( nx - outflow_damping_width ) ) /   &
-                            REAL( outflow_damping_width/2 )            &
-                                             )
-          ELSE
-             km_damp_x(i) = 0.0
+!-- Initialize damping zone for the potential temperature in case of
+!-- non-cyclic lateral boundaries. The damping zone has the maximum value 
+!-- at the inflow boundary and decreases to zero at pt_damping_width.
+    ptdf_x = 0.0
+    ptdf_y = 0.0
+    IF ( bc_lr_dirrad .OR. bc_lr_dirneu )  THEN
+       DO i = nxl, nxr
+          IF ( ( i * dx ) < pt_damping_width )  THEN
+             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5 *        &
+                         REAL( pt_damping_width - i * dx ) / (        &
+                         REAL( pt_damping_width )            ) ) )**2 
           ENDIF
        ENDDO
-
-    ELSEIF ( bc_lr_raddir )  THEN
-
-       DO  i = nxlg, nxrg
-          IF ( i <= outflow_damping_width )  THEN
-             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
-                            ( outflow_damping_width - i ) /            &
-                            REAL( outflow_damping_width/2 )            &
-                                             )
-          ELSE
-             km_damp_x(i) = 0.0
+    ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir )  THEN
+       DO i = nxl, nxr
+          IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) )  THEN
+             ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5 *                 &
+                         REAL( ( i - nx ) * dx + pt_damping_width ) / (        &
+                         REAL( pt_damping_width )                     ) ) )**2       
+          ENDIF
+       ENDDO 
+    ELSEIF ( bc_ns_dirrad .OR. bc_ns_dirneu )  THEN
+       DO j = nys, nyn
+          IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) )  THEN
+             ptdf_y(j) = pt_damping_factor * ( SIN( pi * 0.5 *                 &
+                         REAL( ( j - ny ) * dy + pt_damping_width ) / (        &
+                         REAL( pt_damping_width )                     ) ) )**2       
+          ENDIF
+       ENDDO 
+    ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir )  THEN
+       DO j = nys, nyn
+          IF ( ( j * dy ) < pt_damping_width )  THEN
+             ptdf_y(j) = pt_damping_factor * ( SIN( pi * 0.5 *        &
+                         REAL( pt_damping_width - j * dy ) / (        &
+                         REAL( pt_damping_width )            ) ) )**2       
           ENDIF
        ENDDO
-
-    ENDIF
-
-    IF ( bc_ns_dirrad )  THEN
-
-       DO  j = nysg, nyng
-          IF ( j >= ny - outflow_damping_width )  THEN
-             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
-                            ( j - ( ny - outflow_damping_width ) ) /   &
-                            REAL( outflow_damping_width/2 )            &
-                                             )
-          ELSE
-             km_damp_y(j) = 0.0
-          ENDIF
-       ENDDO
-
-    ELSEIF ( bc_ns_raddir )  THEN
-
-       DO  j = nysg, nyng
-          IF ( j <= outflow_damping_width )  THEN
-             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
-                            ( outflow_damping_width - j ) /            &
-                            REAL( outflow_damping_width/2 )            &
-                                             )
-          ELSE
-             km_damp_y(j) = 0.0
-          ENDIF
-       ENDDO
-
     ENDIF
 
Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 977)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 978)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: nzb_max is set to nzt at non-cyclic lateral boundaries
+! Bugfix: Set wall_flags_0 for inflow boundary
 !
 ! Former revisions:
@@ -608,5 +609,12 @@
 !-- Determine the maximum level of topography. Furthermore it is used for
 !-- steering the degradation of order of the applied advection scheme.
+!-- In case of non-cyclic lateral boundaries, the order of the advection
+!-- scheme is reduced at the lateral boundaries up to nzt.
     nzb_max = MAXVAL( nzb_local )
+    IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR.    &
+         inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s )  THEN
+         nzb_max = nzt
+    ENDIF
+
 !
 !-- Consistency checks and index array initialization are only required for
@@ -1079,10 +1087,12 @@
 !--             scalar - x-direction
 !--             WS1 (0), WS3 (1), WS5 (2)
-                IF ( k <= nzb_s_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) &
-                     .OR. ( outflow_r .AND. i == nxr ) )  THEN
+                IF ( k <= nzb_s_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )&
+                     .AND. i == nxl ) .OR. ( ( inflow_r .OR. outflow_r )       &
+                     .AND. i == nxr ) )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 )
                 ELSEIF ( k <= nzb_s_inner(j,i+2) .OR. k <= nzb_s_inner(j,i-1)  &
-                         .OR. ( outflow_r .AND. i == nxr-1 )                   &
-                         .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
+                         .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) &
+                         .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu  ) &
+                       )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 )
                 ELSE
@@ -1092,11 +1102,13 @@
 !--             scalar - y-direction
 !--             WS1 (3), WS3 (4), WS5 (5)
-                IF ( k <= nzb_s_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) &
-                     .OR. ( outflow_n .AND. j == nyn ) )  THEN
+                IF ( k <= nzb_s_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )&
+                     .AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n )       &
+                     .AND. j == nyn ) )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 )
 !--             WS3
                 ELSEIF ( k <= nzb_s_inner(j+2,i) .OR. k <= nzb_s_inner(j-1,i)  &
-                        .OR. ( outflow_s .AND. j == nysv )                     &
-                        .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
+                         .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv  ) &
+                         .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) &
+                       )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 )
 !--             WS5
@@ -1135,22 +1147,28 @@
 !--             WS1 (9), WS3 (10), WS5 (11)
                 IF ( k <= nzb_u_inner(j,i+1)                                  &
-                     .OR. ( outflow_r .AND. i == nxr ) )  THEN
+                     .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu )     &
+                     .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr  )     &
+                   )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 )
                 ELSEIF ( k <= nzb_u_inner(j,i+2) .OR. k <= nzb_u_inner(j,i-1) &
-                         .OR. ( outflow_r .AND. i == nxr-1 )                  &
-                         .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
+                         .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 )&
+                         .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu+1)&
+                       )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 )
                 ELSE
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 )
                 ENDIF
+
 !
 !--             u component - y-direction
 !--             WS1 (12), WS3 (13), WS5 (14)
-                IF ( k <= nzb_u_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) &
-                    .OR. ( outflow_n .AND. j == nyn ) )  THEN
+                IF ( k <= nzb_u_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )&
+                     .AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n )       &
+                     .AND. j == nyn ) )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
                 ELSEIF ( k <= nzb_u_inner(j+2,i) .OR. k <= nzb_u_inner(j-1,i)  &
-                    .OR. ( outflow_s .AND. j == nysv )                         &
-                    .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
+                         .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv  ) &
+                         .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) &
+                       )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 )
                 ELSE
@@ -1181,11 +1199,13 @@
 !--             v component - x-direction
 !--             WS1 (18), WS3 (19), WS5 (20)
-                IF ( k <= nzb_v_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) &
-                    .OR. ( outflow_r .AND. i == nxr ) )  THEN
+                IF ( k <= nzb_v_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )&
+                     .AND. i == nxl ) .OR. (( inflow_r .OR. outflow_r )        &
+                     .AND. i == nxr ) )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
 !--             WS3
                 ELSEIF ( k <= nzb_v_inner(j,i+2) .OR. k <= nzb_v_inner(j,i-1)  &
-                    .OR. ( outflow_r .AND. i == nxr-1 )                        &
-                    .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
+                         .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) &
+                         .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu  ) &
+                       )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 )
                 ELSE
@@ -1196,9 +1216,12 @@
 !--             WS1 (21), WS3 (22), WS5 (23)
                 IF ( k <= nzb_v_inner(j+1,i)                                   &
-                     .OR. ( outflow_n .AND. j == nyn ) )  THEN
+                     .OR. ( ( inflow_s .OR. outflow_s ) .AND. i == nysv )      &
+                     .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn  )      &
+                   )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
                 ELSEIF ( k <= nzb_v_inner(j+2,i) .OR. k <= nzb_v_inner(j-1,i)  &
-                         .OR. ( outflow_s .AND. j == nysv )                    &
-                         .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
+                         .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv+1 )&
+                         .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1  )&
+                       )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 )
                 ELSE
@@ -1228,10 +1251,12 @@
 !--             w component - x-direction
 !--             WS1 (27), WS3 (28), WS5 (29)
-                IF ( k <= nzb_w_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) &
-                    .OR. ( outflow_r .AND. i == nxr ) )  THEN
+                IF ( k <= nzb_w_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )&
+                     .AND. i == nxl ) .OR. ( ( inflow_r .OR. outflow_r )       &
+                     .AND. i == nxr ) )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
                 ELSEIF ( k <= nzb_w_inner(j,i+2) .OR. k <= nzb_w_inner(j,i-1)  &
-                         .OR. ( outflow_r .AND. i == nxr-1 )                   &
-                         .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
+                         .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) &
+                         .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu  ) &
+                       )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 )
                 ELSE
@@ -1241,10 +1266,12 @@
 !--             w component - y-direction
 !--             WS1 (30), WS3 (31), WS5 (32)
-                IF ( k <= nzb_w_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) &
-                     .OR. ( outflow_n .AND. j == nyn ) )  THEN
+                IF ( k <= nzb_w_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )&
+                     .AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n )       &
+                     .AND. j == nyn ) )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
                 ELSEIF ( k <= nzb_w_inner(j+2,i) .OR. k <= nzb_w_inner(j-1,i)  &
-                         .OR. ( outflow_s .AND. j == nysv )                    &
-                         .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
+                         .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv  ) &
+                         .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) &
+                       )  THEN
                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 )
                 ELSE
Index: /palm/trunk/SOURCE/init_masks.f90
===================================================================
--- /palm/trunk/SOURCE/init_masks.f90	(revision 977)
+++ /palm/trunk/SOURCE/init_masks.f90	(revision 978)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! +z0h*
 !
 ! Former revisions:
@@ -257,5 +257,5 @@
                 unit = 'psu'
 
-             CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*' )
+             CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*', 'z0h*' )
                 WRITE ( message_string, * ) 'illegal value for data_',&
                      'output: "', TRIM( var ), '" is only allowed', &
Index: /palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- /palm/trunk/SOURCE/init_pegrid.f90	(revision 977)
+++ /palm/trunk/SOURCE/init_pegrid.f90	(revision 978)
@@ -4,4 +4,6 @@
 ! Current revisions:
 ! -----------------
+! dirichlet/neumann and neumann/dirichlet added
+! nxlu and nysv are also calculated for inflow boundary
 !
 ! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
@@ -1108,7 +1110,7 @@
 !-- horizontal boundary conditions.
     IF ( pleft == MPI_PROC_NULL )  THEN
-       IF ( bc_lr == 'dirichlet/radiation' )  THEN
+       IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' )  THEN
           inflow_l  = .TRUE.
-       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
+       ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' )  THEN
           outflow_l = .TRUE.
        ENDIF
@@ -1116,7 +1118,7 @@
 
     IF ( pright == MPI_PROC_NULL )  THEN
-       IF ( bc_lr == 'dirichlet/radiation' )  THEN
+       IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' )  THEN
           outflow_r = .TRUE.
-       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
+       ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' )  THEN
           inflow_r  = .TRUE.
        ENDIF
@@ -1124,7 +1126,7 @@
 
     IF ( psouth == MPI_PROC_NULL )  THEN
-       IF ( bc_ns == 'dirichlet/radiation' )  THEN
+       IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' )  THEN
           outflow_s = .TRUE.
-       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
+       ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' )  THEN
           inflow_s  = .TRUE.
        ENDIF
@@ -1132,7 +1134,7 @@
 
     IF ( pnorth == MPI_PROC_NULL )  THEN
-       IF ( bc_ns == 'dirichlet/radiation' )  THEN
+       IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' )  THEN
           inflow_n  = .TRUE.
-       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
+       ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' )  THEN
           outflow_n = .TRUE.
        ENDIF
@@ -1164,16 +1166,16 @@
 
 #elif ! defined ( __parallel )
-    IF ( bc_lr == 'dirichlet/radiation' )  THEN
+    IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' )  THEN
        inflow_l  = .TRUE.
        outflow_r = .TRUE.
-    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
+    ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' )  THEN
        outflow_l = .TRUE.
        inflow_r  = .TRUE.
     ENDIF
 
-    IF ( bc_ns == 'dirichlet/radiation' )  THEN
+    IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' )  THEN
        inflow_n  = .TRUE.
        outflow_s = .TRUE.
-    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
+    ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' )  THEN
        outflow_n = .TRUE.
        inflow_s  = .TRUE.
@@ -1182,12 +1184,12 @@
 
 !
-!-- At the outflow, u or v, respectively, have to be calculated for one more
-!-- grid point.
-    IF ( outflow_l )  THEN
+!-- At the inflow or outflow, u or v, respectively, have to be calculated for
+!-- one more grid point.
+    IF ( inflow_l .OR. outflow_l )  THEN
        nxlu = nxl + 1
     ELSE
        nxlu = nxl
     ENDIF
-    IF ( outflow_s )  THEN
+    IF ( inflow_s .OR. outflow_s )  THEN
        nysv = nys + 1
     ELSE
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 977)
+++ /palm/trunk/SOURCE/modules.f90	(revision 978)
@@ -4,5 +4,9 @@
 ! Current revisions:
 ! -----------------
-!
+! +c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l,
+! +bc_lr_dirneu, bc_lr_neudir, bc_ns_dirneu, bc_ns_neudir
+! -km_damp_x, km_damp_y, km_damp_max, outflow_damping_width
+! +z0h, z0h_av, z0h_factor, z0h1d
+! +ptdf_x, ptdf_y, pt_damping_width, pt_damping_factor 
 !
 ! Former revisions:
@@ -324,15 +328,15 @@
 
     REAL, DIMENSION(:), ALLOCATABLE ::                                         &
-          ddzu, ddzu_pres, dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor,  &
-          km_damp_x, km_damp_y, lad, l_grid, pt_init, q_init, rdf, rdf_sc,     &
-          sa_init, ug, u_init, u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, &
-          w_subs, zu, zw
+          c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l, ddzu, ddzu_pres,     &
+          dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, lad, l_grid,      &
+          ptdf_x, ptdf_y, pt_init, q_init, rdf, rdf_sc, sa_init, ug, u_init,   &
+          u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw 
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
-          c_u, c_v, c_w, dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg,                  &
-          mean_inflow_profiles, pt_slope_ref, qs, qswst_remote, ts, us, z0,    &
-          flux_s_pt, diss_s_pt, flux_s_e, diss_s_e, flux_s_q, diss_s_q,        &
-          flux_s_sa, diss_s_sa, flux_s_u, flux_s_v, flux_s_w, diss_s_u,        &
-          diss_s_v, diss_s_w, total_2d_o, total_2d_a
+          c_u, c_v, c_w, diss_s_e, diss_s_pt, diss_s_q, diss_s_sa, diss_s_u,   &
+          diss_s_v, diss_s_w, dzu_mg, dzw_mg, flux_s_e, flux_s_pt, flux_s_q,   &
+          flux_s_sa, flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,        &
+          mean_inflow_profiles, pt_slope_ref, qs, qswst_remote, total_2d_a,    &
+          total_2d_o, ts, us, z0, z0h
 
     REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
@@ -347,10 +351,10 @@
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
-          canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, lad_s, lad_u,   &
-          lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, u_m_l, u_m_n,      &
-          u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r,       &
-          w_m_s, flux_l_pt, diss_l_pt, flux_l_e, diss_l_e, flux_l_q, diss_l_q, &
-          flux_l_sa, diss_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_l_u,        &
-          diss_l_v, diss_l_w
+          canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, diss_l_e,       &
+          diss_l_pt, diss_l_q, diss_l_sa, diss_l_u, diss_l_v, diss_l_w,        &
+          flux_l_e, flux_l_pt, flux_l_q, flux_l_sa, flux_l_u, flux_l_v,        &
+          flux_l_w, lad_s, lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls,  &
+          tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l, &
+          w_m_n, w_m_r, w_m_s
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
@@ -388,5 +392,6 @@
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::  lwp_av, precipitation_rate_av,       &
-                                          qsws_av, shf_av,ts_av, us_av, z0_av
+                                          qsws_av, shf_av,ts_av, us_av, z0_av, &
+                                          z0h_av
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
@@ -549,6 +554,6 @@
                 netcdf_data_format = 2, ngsrb = 2, nsor = 20, &
                 nsor_ini = 100, n_sor, normalizing_region = 0, &
-                nz_do3d = -9999, outflow_damping_width = -1, &
-                pch_index = 0, prt_time_count = 0, recycling_plane, runnr = 0, &
+                nz_do3d = -9999, pch_index = 0, prt_time_count = 0, &
+                recycling_plane, runnr = 0, &
                 skip_do_avs = 0, subdomain_size, terminate_coupled = 0, &
                 terminate_coupled_remote = 0, timestep_count = 0
@@ -582,7 +587,9 @@
 
     LOGICAL ::  adjust_mixing_length = .FALSE., avs_output = .FALSE., &
-                bc_lr_cyc =.TRUE., bc_lr_dirrad = .FALSE., &
+                bc_lr_cyc =.TRUE., bc_lr_dirneu = .FALSE., &
+                bc_lr_dirrad = .FALSE., bc_lr_neudir = .FALSE., &
                 bc_lr_raddir = .FALSE., bc_ns_cyc = .TRUE., &
-                bc_ns_dirrad = .FALSE., bc_ns_raddir = .FALSE., &
+                bc_ns_dirneu = .FALSE., bc_ns_dirrad = .FALSE., &
+                bc_ns_neudir = .FALSE., bc_ns_raddir = .FALSE., &
                 call_psolver_at_all_substeps = .TRUE., &
                 cloud_droplets = .FALSE., cloud_physics = .FALSE., &
@@ -659,5 +666,5 @@
              f = 0.0, fs = 0.0, g = 9.81, inflow_damping_height = 9999999.9, &
              inflow_damping_width = 9999999.9, kappa = 0.4, km_constant = -1.0,&
-             km_damp_max = -1.0, lad_surface = 0.0,  &
+             lad_surface = 0.0,  &
              leaf_surface_concentration = 0.0, long_filter_factor = 0.0, &
              mask_scale_x = 1.0, mask_scale_y = 1.0, mask_scale_z = 1.0, &
@@ -670,4 +677,5 @@
              phi = 55.0, prandtl_number = 1.0, &
              precipitation_amount_interval = 9999999.9, prho_reference, &
+             pt_damping_factor = 0.0, pt_damping_width = 0.0, &
              pt_reference = 9999999.9, pt_slope_offset = 0.0, &
              pt_surface = 300.0, pt_surface_initial_change = 0.0, &
@@ -700,5 +708,5 @@
              ups_limit_v = 0.0, ups_limit_w = 0.0, vg_surface = 0.0, &
              v_bulk = 0.0, v_gtrans = 0.0, wall_adjustment_factor = 1.8, &
-             z_max_do2d = -1.0
+             z_max_do2d = -1.0, z0h_factor = 1.0
 
     REAL ::  do2d_xy_last_time(0:1) = -1.0, do2d_xz_last_time(0:1) = -1.0, &
@@ -1057,5 +1065,5 @@
                 qs1d, simulated_time_1d = 0.0, time_pr_1d = 0.0, &
                 time_run_control_1d = 0.0, ts1d, us1d, usws1d, usws1d_m, &
-                vsws1d, vsws1d_m, z01d
+                vsws1d, vsws1d_m, z01d, z0h1d
 
 
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 977)
+++ /palm/trunk/SOURCE/parin.f90	(revision 978)
@@ -4,5 +4,7 @@
 ! Current revisions:
 ! -----------------
-!
+! -km_damp_max, outflow_damping_width
+! +pt_damping_factor, pt_damping_width
+! +z0h_factor
 !
 ! Former revisions:
@@ -168,5 +170,5 @@
              inflow_damping_height, inflow_damping_width, &
              inflow_disturbance_begin, inflow_disturbance_end, &
-             initializing_actions, km_constant, km_damp_max, lad_surface, &
+             initializing_actions, km_constant, lad_surface, &
              lad_vertical_gradient, lad_vertical_gradient_level, &
              leaf_surface_concentration, long_filter_factor, &
@@ -175,8 +177,9 @@
              netcdf_precision, neutral, ngsrb, nsor, &
              nsor_ini, nx, ny, nz, ocean, omega, omega_sor, &
-             outflow_damping_width, overshoot_limit_e, overshoot_limit_pt, &
+             overshoot_limit_e, overshoot_limit_pt, &
              overshoot_limit_u, overshoot_limit_v, overshoot_limit_w, &
              passive_scalar, pch_index, phi, plant_canopy, prandtl_layer, &
-             prandtl_number, precipitation, psolver, pt_reference, pt_surface, &
+             prandtl_number, precipitation, psolver, pt_damping_factor, &
+             pt_damping_width, pt_reference, pt_surface, &
              pt_surface_initial_change, pt_vertical_gradient, &
              pt_vertical_gradient_level, q_surface, q_surface_initial_change, &
@@ -200,5 +203,5 @@
              uv_heights, u_bulk, u_profile, vg_surface, vg_vertical_gradient, &
              vg_vertical_gradient_level, v_bulk, v_profile, wall_adjustment, &
-             wall_heatflux, wall_humidityflux, wall_scalarflux
+             wall_heatflux, wall_humidityflux, wall_scalarflux, z0h_factor
 
 
Index: /palm/trunk/SOURCE/poismg.f90
===================================================================
--- /palm/trunk/SOURCE/poismg.f90	(revision 977)
+++ /palm/trunk/SOURCE/poismg.f90	(revision 978)
@@ -7,5 +7,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! bc_lr/ns_dirneu/neudir added
+!
 !
 ! Former revisions:
@@ -1229,10 +1230,10 @@
 !--       outflow conditions have to be used on all PEs after the switch,
 !--       because then they have the total domain.
-          IF ( bc_lr_dirrad )  THEN
+          IF ( bc_lr_dirrad .OR. bc_lr_dirneu )  THEN
              inflow_l  = .TRUE.
              inflow_r  = .FALSE.
              outflow_l = .FALSE.
              outflow_r = .TRUE.
-          ELSEIF ( bc_lr_raddir )  THEN
+          ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir )  THEN
              inflow_l  = .FALSE.
              inflow_r  = .TRUE.
@@ -1241,10 +1242,10 @@
           ENDIF
 
-          IF ( bc_ns_dirrad )  THEN
+          IF ( bc_ns_dirrad .OR. bc_ns_dirneu )  THEN
              inflow_n  = .TRUE.
              inflow_s  = .FALSE.
              outflow_n = .FALSE.
              outflow_s = .TRUE.
-          ELSEIF ( bc_ns_raddir )  THEN
+          ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir )  THEN
              inflow_n  = .FALSE.
              inflow_s  = .TRUE.
@@ -1316,7 +1317,7 @@
 
           IF ( pleft == MPI_PROC_NULL )  THEN
-             IF ( bc_lr_dirrad )  THEN
+             IF ( bc_lr_dirrad .OR. bc_lr_dirneu )  THEN
                 inflow_l  = .TRUE.
-             ELSEIF ( bc_lr_raddir )  THEN
+             ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir )  THEN
                 outflow_l = .TRUE.
              ENDIF
@@ -1324,7 +1325,7 @@
 
           IF ( pright == MPI_PROC_NULL )  THEN
-             IF ( bc_lr_dirrad )  THEN
+             IF ( bc_lr_dirrad .OR. bc_lr_dirneu )  THEN
                 outflow_r = .TRUE.
-             ELSEIF ( bc_lr_raddir )  THEN
+             ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir )  THEN
                 inflow_r  = .TRUE.
              ENDIF
@@ -1332,7 +1333,7 @@
 
           IF ( psouth == MPI_PROC_NULL )  THEN
-             IF ( bc_ns_dirrad )  THEN
+             IF ( bc_ns_dirrad .OR. bc_ns_dirneu )  THEN
                 outflow_s = .TRUE.
-             ELSEIF ( bc_ns_raddir )  THEN
+             ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir )  THEN
                 inflow_s  = .TRUE.
              ENDIF
@@ -1340,7 +1341,7 @@
 
           IF ( pnorth == MPI_PROC_NULL )  THEN
-             IF ( bc_ns_dirrad )  THEN
+             IF ( bc_ns_dirrad .OR. bc_ns_dirneu )  THEN
                 inflow_n  = .TRUE.
-             ELSEIF ( bc_ns_raddir )  THEN
+             ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir )  THEN
                 outflow_n = .TRUE.
              ENDIF
Index: /palm/trunk/SOURCE/prandtl_fluxes.f90
===================================================================
--- /palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 977)
+++ /palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 978)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! roughness length for scalar quantities z0h added
 !
 ! Former revisions:
@@ -96,7 +96,7 @@
 !
 !--             Stable stratification
-                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
-                                  LOG( z_p / z0(j,i) ) +                   &
-                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
+                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (           &
+                                  LOG( z_p / z0h(j,i) ) +                   &
+                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
                                                                 )
              ELSE
@@ -104,8 +104,8 @@
 !--             Unstable stratification
                 a = SQRT( 1.0 - 16.0 * rif(j,i) )
-                b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p )
-
-                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (         &
-                          LOG( z_p / z0(j,i) ) -                           &
+                b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p )
+
+                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (          &
+                          LOG( z_p / z0h(j,i) ) -                           &
                           2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) )
              ENDIF
@@ -308,7 +308,7 @@
 !
 !--                Stable stratification
-                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
-                                  LOG( z_p / z0(j,i) ) +                   &
-                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
+                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
+                                  LOG( z_p / z0h(j,i) ) +                   &
+                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
                                                                  )
                 ELSE
@@ -316,8 +316,8 @@
 !--                Unstable stratification
                    a = SQRT( 1.0 - 16.0 * rif(j,i) ) 
-                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p ) 
+                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p ) 
  
-                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (   &
-                             LOG( z_p / z0(j,i) ) -                    &
+                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (        &
+                             LOG( z_p / z0h(j,i) ) -                        &
                               2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
                 ENDIF
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 977)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 978)
@@ -4,5 +4,8 @@
 ! Current revisions:
 ! -----------------
-!
+! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
+! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
+! boundary in case of non-cyclic lateral boundaries
+! Bugfix: first thread index changes for WS-scheme at the inflow
 !
 ! Former revisions:
@@ -201,9 +204,9 @@
           ENDIF
           IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
-             CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m,   &
-                               usws_m, uswst_m, v_m, w_m )
+             CALL diffusion_u( i, j, ddzu, ddzw, km_m, tend, u_m, usws_m, &
+                               uswst_m, v_m, w_m )
           ELSE
-             CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
-                               uswst, v, w )
+             CALL diffusion_u( i, j, ddzu, ddzw, km, tend, u, usws, uswst, &
+                               v, w )
           ENDIF
           CALL coriolis( i, j, 1 )
@@ -289,9 +292,9 @@
           ENDIF
           IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
-             CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, &
-                               v_m, vsws_m, vswst_m, w_m )
+             CALL diffusion_v( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, &
+                               vsws_m, vswst_m, w_m )
           ELSE
-             CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v,  &
-                               vsws, vswst, w )
+             CALL diffusion_v( i, j, ddzu, ddzw, km, tend, u, v, vsws, &
+                               vswst, w )
           ENDIF
           CALL coriolis( i, j, 2 )
@@ -373,9 +376,7 @@
           ENDIF
           IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
-             CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, km_damp_y,  &
-                               tend, u_m, v_m, w_m )
+             CALL diffusion_w( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, w_m )
           ELSE
-             CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y,    &
-                               tend, u, v, w )
+             CALL diffusion_w( i, j, ddzu, ddzw, km, tend, u, v, w )
           ENDIF
           CALL coriolis( i, j, 3 )
@@ -524,8 +525,9 @@
              DO  k = nzb_s_inner(j,i)+1, nzt
                 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
-                              dt_3d * (                                           &
+                              dt_3d * (                                        &
                                        sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
-                                      ) -                                         &
-                              tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
+                                      ) - &
+                              tsc(5) * ( pt(k,j,i) - pt_init(k) ) * &
+                              ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )  
              ENDDO
 
@@ -983,5 +985,5 @@
              IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
                 IF ( ws_scheme_mom )  THEN
-                   IF ( outflow_l .AND. i_omp_start == nxl )  THEN 
+                   IF ( ( inflow_l .OR. outflow_l ) .AND. i_omp_start == nxl )  THEN 
 !                     CALL local_diss( i, j, u)    ! dissipation control
                       CALL advec_u_ws( i, j, i_omp_start + 1, tn )
@@ -997,9 +999,9 @@
              IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
              THEN
-                CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend,  &
-                                  u_m, usws_m, uswst_m, v_m, w_m )
-             ELSE
-                CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
-                                  usws, uswst, v, w )
+                CALL diffusion_u( i, j, ddzu, ddzw, km_m, tend,  u_m, &
+                                  usws_m, uswst_m, v_m, w_m )
+             ELSE
+                CALL diffusion_u( i, j, ddzu, ddzw, km, tend, u, usws, &
+                                  uswst, v, w )
              ENDIF
              CALL coriolis( i, j, 1 )
@@ -1066,9 +1068,9 @@
              IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
              THEN
-                CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend,     &
-                                  u_m, v_m, vsws_m, vswst_m, w_m )
-             ELSE
-                CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
-                                  vsws, vswst, w )
+                CALL diffusion_v( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, &
+                                  vsws_m, vswst_m, w_m )
+             ELSE
+                CALL diffusion_v( i, j, ddzu, ddzw, km, tend, u, v, vsws, &
+                                  vswst, w )
              ENDIF
              CALL coriolis( i, j, 2 )
@@ -1130,9 +1132,8 @@
           IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
           THEN
-             CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x,          &
-                               km_damp_y, tend, u_m, v_m, w_m )
+             CALL diffusion_w( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, &
+                               w_m )
           ELSE
-             CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
-                               tend, u, v, w )
+             CALL diffusion_w( i, j, ddzu, ddzw, km, tend, u, v, w )
           ENDIF
           CALL coriolis( i, j, 3 )
@@ -1240,6 +1241,7 @@
                               dt_3d * (                                        &
                                   tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
-                                      ) -                                      &
-                              tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
+                                      ) - &
+                              tsc(5) * ( pt(k,j,i) - pt_init(k) ) * &
+                              ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )  
              ENDDO
 
@@ -1540,8 +1542,8 @@
     ENDIF
     IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
-       CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, &
-                         uswst_m, v_m, w_m )
+       CALL diffusion_u( ddzu, ddzw, km_m, tend, u_m, usws_m, uswst_m, &
+                         v_m, w_m )
     ELSE
-       CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, v, w )
+       CALL diffusion_u( ddzu, ddzw, km, tend, u, usws, uswst, v, w )
     ENDIF
     CALL coriolis( 1 )
@@ -1634,8 +1636,8 @@
     ENDIF
     IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
-       CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
-                         vswst_m, w_m )
+       CALL diffusion_v( ddzu, ddzw, km_m, tend, u_m, v_m, vsws_m, vswst_m, &
+                         w_m )
     ELSE
-       CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, vswst, w )
+       CALL diffusion_v( ddzu, ddzw, km, tend, u, v, vsws, vswst, w )
     ENDIF
     CALL coriolis( 2 )
@@ -1725,8 +1727,7 @@
     ENDIF
     IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
-       CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
-                         v_m, w_m )
+       CALL diffusion_w( ddzu, ddzw, km_m, tend, u_m, v_m, w_m )
     ELSE
-       CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
+       CALL diffusion_w( ddzu, ddzw, km, tend, u, v, w )
     ENDIF
     CALL coriolis( 3 )
@@ -1877,9 +1878,10 @@
           DO  j = nys, nyn
              DO  k = nzb_s_inner(j,i)+1, nzt
-                pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) +       &
-                              dt_3d * (                                           &
-                                       sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
-                                      ) -                                         &
-                              tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
+                pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) +    &
+                              dt_3d * (                                        &
+                                       sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i)&
+                                      ) - &
+                              tsc(5) * ( pt(k,j,i) - pt_init(k) ) * &
+                              ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )  
              ENDDO
           ENDDO
Index: /palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/read_3d_binary.f90	(revision 977)
+++ /palm/trunk/SOURCE/read_3d_binary.f90	(revision 978)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! +z0h, z0h_av
 !
 ! Former revisions:
@@ -963,4 +964,17 @@
                        tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+                CASE ( 'z0h' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   z0h(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                     tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'z0h_av' )
+                   IF ( .NOT. ALLOCATED( z0h_av ) )  THEN
+                      ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   z0h_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                       tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
                 CASE DEFAULT
                    WRITE( message_string, * ) 'unknown field named "', &
Index: /palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/read_var_list.f90	(revision 977)
+++ /palm/trunk/SOURCE/read_var_list.f90	(revision 978)
@@ -4,5 +4,7 @@
 ! Current revisions:
 ! ------------------
-!
+! -km_damp_max, outflow_damping_width
+! +pt_damping_factor, pt_damping_width
+! +z0h_factor
 !
 ! Former revisions:
@@ -385,6 +387,4 @@
           CASE ( 'km_constant' )
              READ ( 13 )  km_constant
-          CASE ( 'km_damp_max' )
-             READ ( 13 )  km_damp_max
           CASE ( 'lad' )
              READ ( 13 )  lad
@@ -446,6 +446,4 @@
           CASE ( 'omega_sor' )
              READ ( 13 )  omega_sor
-          CASE ( 'outflow_damping_width' )
-             READ ( 13 )  outflow_damping_width
           CASE ( 'output_for_t0' )
              READ (13)    output_for_t0
@@ -476,4 +474,8 @@
           CASE ( 'psolver' )
              READ ( 13 )  psolver
+          CASE ( 'pt_damping_factor' )
+             READ ( 13 )  pt_damping_factor
+          CASE ( 'pt_damping_width' )
+             READ ( 13 )  pt_damping_width
           CASE ( 'pt_init' )
              READ ( 13 )  pt_init
@@ -684,4 +686,6 @@
           CASE ( 'w_max_ijk' )
              READ ( 13 )  w_max_ijk
+          CASE ( 'z0h_factor' )
+             READ ( 13 )  z0h_factor
 
           CASE DEFAULT
Index: /palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 977)
+++ /palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 978)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! +z0h*
 !
 ! Former revisions:
@@ -227,4 +227,10 @@
                 ENDIF
                 z0_av = 0.0
+
+             CASE ( 'z0h*' )
+                IF ( .NOT. ALLOCATED( z0h_av ) )  THEN
+                   ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                z0h_av = 0.0
 
              CASE DEFAULT
@@ -494,4 +500,11 @@
              ENDDO
 
+          CASE ( 'z0h*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   z0h_av(j,i) = z0h_av(j,i) + z0h(j,i)
+                ENDDO
+             ENDDO
+
           CASE DEFAULT
 !
Index: /palm/trunk/SOURCE/timestep.f90
===================================================================
--- /palm/trunk/SOURCE/timestep.f90	(revision 977)
+++ /palm/trunk/SOURCE/timestep.f90	(revision 978)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! restriction of the outflow damping layer in the diffusion criterion removed
 !
 ! Former revisions:
@@ -185,14 +185,4 @@
        dt_diff = dt_diff_l
 #endif
-
-!
-!--    In case of non-cyclic lateral boundaries, the diffusion time step
-!--    may be further restricted by the lateral damping layer (damping only
-!--    along x and y)
-       IF ( .NOT. bc_lr_cyc )  THEN
-          dt_diff = MIN( dt_diff, 0.125 * dx2 / ( km_damp_max + 1E-20 ) )
-       ELSEIF ( .NOT. bc_ns_cyc )  THEN
-          dt_diff = MIN( dt_diff, 0.125 * dy2 / ( km_damp_max + 1E-20 ) )
-       ENDIF
 
 !
Index: /palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/write_3d_binary.f90	(revision 977)
+++ /palm/trunk/SOURCE/write_3d_binary.f90	(revision 978)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! +z0h, z0h_av
 !
 ! Former revisions:
@@ -298,4 +298,8 @@
        WRITE ( 14 )  'z0_av               ';  WRITE ( 14 )  z0_av
     ENDIF
+    WRITE ( 14 )  'z0h                 ';  WRITE ( 14 )  z0h
+    IF ( ALLOCATED( z0h_av ) )  THEN
+       WRITE ( 14 )  'z0h_av              ';  WRITE ( 14 )  z0h_av
+    ENDIF
 
 !
Index: /palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/write_var_list.f90	(revision 977)
+++ /palm/trunk/SOURCE/write_var_list.f90	(revision 978)
@@ -4,5 +4,7 @@
 ! Current revisions:
 ! -----------------
-!
+! -km_damp_max, outflow_damping_width
+! +pt_damping_factor, pt_damping_width
+! +z0h_factor
 !
 ! Former revisions:
@@ -300,6 +302,4 @@
     WRITE ( 14 )  'km_constant                   '
     WRITE ( 14 )  km_constant
-    WRITE ( 14 )  'km_damp_max                   '
-    WRITE ( 14 )  km_damp_max
     WRITE ( 14 )  'lad                           '
     WRITE ( 14 )  lad
@@ -358,6 +358,4 @@
     WRITE ( 14 )  'omega_sor                     '
     WRITE ( 14 )  omega_sor
-    WRITE ( 14 )  'outflow_damping_width         '
-    WRITE ( 14 )  outflow_damping_width
     WRITE ( 14 )  'output_for_t0                 '
     WRITE ( 14 )  output_for_t0
@@ -388,4 +386,8 @@
     WRITE ( 14 )  'psolver                       '
     WRITE ( 14 )  psolver
+    WRITE ( 14 )  'pt_damping_factor             '
+    WRITE ( 14 )  pt_damping_factor
+    WRITE ( 14 )  'pt_damping_width              '
+    WRITE ( 14 )  pt_damping_width
     WRITE ( 14 )  'pt_init                       '
     WRITE ( 14 )  pt_init
@@ -596,4 +598,6 @@
     WRITE ( 14 )  'w_max_ijk                     '
     WRITE ( 14 )  w_max_ijk
+    WRITE ( 14 )  'z0h_factor                    '
+    WRITE ( 14 )  z0h_factor
 
 !
