Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 105)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 106)
@@ -13,5 +13,9 @@
 inipar parameter top_momentumflux_u|v.
 
-boundary_conds, check_open, check_parameters, diffusion_u, diffusion_v, flow_statistics, header, init_pegrid, init_3d_model, modules, palm, parin, prognostic_equations, read_var_list, read_3d_binary, swap_timelevel, time_integration, write_var_list, write_3d_binary
+Non-cyclic boundary conditions can be used along all horizontal directions.
+
+Quantities w*p* and w"e can be output as vertical profiles.
+
+boundary_conds, check_open, check_parameters, diffusion_u, diffusion_v, flow_statistics, header, init_pegrid, init_3d_model, modules, palm, parin, pres, prognostic_equations, read_var_list, read_3d_binary, swap_timelevel, time_integration, write_var_list, write_3d_binary
 
 New:
@@ -21,11 +25,31 @@
 Changed:
 -------
+Remaining variables iran changed to iran_part (advec_particles, init_particles).
+
+In case that the presure solver is not called for every Runge-Kutta substep
+(call_psolver_at_all_substeps = .F.), it is called after the first substep
+instead of the last. In that case, random perturbations are also added to the
+velocity field after the first substep.
+
+advec_particles, init_particles, time_integration
 
 
 Errors:
 ------
-+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface.f90)
+Bugs from code parts for non-cyclic boundary conditions are removed: loops for
+u and v are starting from index nxlu, nysv, respectively. The radiation boundary
+condition is used for every Runge-Kutta substep. Velocity phase speeds for
+the radiation boundary conditions are calculated for the first Runge-Kutta
+substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
+are defined for this purpose. Several index errors are removed from the
+radiation boundary condition code parts. Upper bounds for calculating
+u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
+values are not available in case of non-cyclic boundary conditions.
+
++dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)
+
+Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)
 
 Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).
 
-check_parameters, user_interface
+advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, boundary_conds, buoyancy, check_parameters, coriolis, diffusion_u, diffusion_v, init_pegrid, init_3d_model, modules, production_e, prognostic_equations, user_interface
Index: palm/trunk/SOURCE/advec_particles.f90
===================================================================
--- palm/trunk/SOURCE/advec_particles.f90	(revision 105)
+++ palm/trunk/SOURCE/advec_particles.f90	(revision 106)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! remaining variables iran changed to iran_part
 ! TEST: PRINT statements on unit 9 (commented out)
 !
@@ -1955,5 +1955,5 @@
                    THEN
                       particles(n)%x = particles(n)%x + &
-                                       ( random_function( iran ) - 0.5 ) * &
+                                       ( random_function( iran_part ) - 0.5 ) *&
                                        pdx(particles(n)%group)
                       IF ( particles(n)%x  <=  ( nxl - 0.5 ) * dx )  THEN
@@ -1966,5 +1966,5 @@
                    THEN
                       particles(n)%y = particles(n)%y + &
-                                       ( random_function( iran ) - 0.5 ) * &
+                                       ( random_function( iran_part ) - 0.5 ) *&
                                        pdy(particles(n)%group)
                       IF ( particles(n)%y  <=  ( nys - 0.5 ) * dy )  THEN
@@ -1977,5 +1977,5 @@
                    THEN
                       particles(n)%z = particles(n)%z + &
-                                       ( random_function( iran ) - 0.5 ) * &
+                                       ( random_function( iran_part ) - 0.5 ) *&
                                        pdz(particles(n)%group)
                    ENDIF
Index: palm/trunk/SOURCE/advec_u_pw.f90
===================================================================
--- palm/trunk/SOURCE/advec_u_pw.f90	(revision 105)
+++ palm/trunk/SOURCE/advec_u_pw.f90	(revision 106)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
 !
 ! Former revisions:
@@ -58,5 +58,5 @@
        gu = 2.0 * u_gtrans
        gv = 2.0 * v_gtrans
-       DO  i = nxl, nxr
+       DO  i = nxlu, nxr
           DO  j = nys, nyn
              DO  k = nzb_u_inner(j,i)+1, nzt
Index: palm/trunk/SOURCE/advec_u_up.f90
===================================================================
--- palm/trunk/SOURCE/advec_u_up.f90	(revision 105)
+++ palm/trunk/SOURCE/advec_u_up.f90	(revision 106)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
 !
 ! Former revisions:
@@ -57,5 +57,5 @@
 
 
-       DO  i = nxl, nxr
+       DO  i = nxlu, nxr
           DO  j = nys, nyn
              DO  k = nzb_u_inner(j,i)+1, nzt
Index: palm/trunk/SOURCE/advec_v_pw.f90
===================================================================
--- palm/trunk/SOURCE/advec_v_pw.f90	(revision 105)
+++ palm/trunk/SOURCE/advec_v_pw.f90	(revision 106)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! j loop is starting from nysv (needed for non-cyclic boundary conditions)
 !
 ! Former revisions:
@@ -60,5 +60,5 @@
        gv = 2.0 * v_gtrans
        DO  i = nxl, nxr
-          DO  j = nys, nyn
+          DO  j = nysv, nyn
              DO  k = nzb_v_inner(j,i)+1, nzt
                 tend(k,j,i) = tend(k,j,i) - 0.25 * (                           &
Index: palm/trunk/SOURCE/advec_v_up.f90
===================================================================
--- palm/trunk/SOURCE/advec_v_up.f90	(revision 105)
+++ palm/trunk/SOURCE/advec_v_up.f90	(revision 106)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! j loop is starting from nysv (needed for non-cyclic boundary conditions)
 !
 ! Former revisions:
@@ -57,5 +57,5 @@
 
        DO  i = nxl, nxr
-          DO  j = nys, nyn
+          DO  j = nysv, nyn
              DO  k = nzb_v_inner(j,i)+1, nzt
 !
Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 105)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 106)
@@ -4,5 +4,9 @@
 ! Actual revisions:
 ! -----------------
-! Boundary conditions for temperature adjusted for coupled runs
+! Boundary conditions for temperature adjusted for coupled runs,
+! bugfixes for the radiation boundary conditions at the outflow: radiation
+! conditions are used for every substep, phase speeds are calculated for the
+! first Runge-Kutta substep only and then reused, several index values changed
+! 
 !
 ! Former revisions:
@@ -59,5 +63,5 @@
     INTEGER ::  i, j, k
 
-    REAL    ::  c_max, c_u, c_v, c_w, denom
+    REAL    ::  c_max, denom
 
 
@@ -216,7 +220,5 @@
 !
 !-- Radiation boundary condition for the velocities at the respective outflow
-    IF ( outflow_s  .AND.                                                 &
-         intermediate_timestep_count == intermediate_timestep_count_max ) &
-    THEN
+    IF ( outflow_s )  THEN
 
        c_max = dy / dt_3d
@@ -226,59 +228,66 @@
 
 !
-!--          First calculate the phase speeds for u,v, and w
-             denom = u_m_s(k,0,i) - u_m_s(k,1,i)
-
-             IF ( denom /= 0.0 )  THEN
-                c_u = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
-                IF ( c_u < 0.0 )  THEN
-                   c_u = 0.0
-                ELSEIF ( c_u > c_max )  THEN
-                   c_u = c_max
-                ENDIF
-             ELSE
-                c_u = c_max
+!--          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
+
+                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
+                   IF ( c_u(k,i) < 0.0 )  THEN
+                      c_u(k,i) = 0.0
+                   ELSEIF ( c_u(k,i) > c_max )  THEN
+                      c_u(k,i) = c_max
+                   ENDIF
+                ELSE
+                   c_u(k,i) = c_max
+                ENDIF
+
+                denom = v_m_s(k,1,i) - v_m_s(k,2,i)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / denom
+                   IF ( c_v(k,i) < 0.0 )  THEN
+                      c_v(k,i) = 0.0
+                   ELSEIF ( c_v(k,i) > c_max )  THEN
+                      c_v(k,i) = c_max
+                   ENDIF
+                ELSE
+                   c_v(k,i) = c_max
+                ENDIF
+
+                denom = w_m_s(k,0,i) - w_m_s(k,1,i)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
+                   IF ( c_w(k,i) < 0.0 )  THEN
+                      c_w(k,i) = 0.0
+                   ELSEIF ( c_w(k,i) > c_max )  THEN
+                      c_w(k,i) = c_max
+                   ENDIF
+                ELSE
+                   c_w(k,i) = c_max
+                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
-             denom = v_m_s(k,0,i) - v_m_s(k,1,i)
-
-             IF ( denom /= 0.0 )  THEN
-                c_v = -c_max * ( v(k,0,i) - v_m_s(k,0,i) ) / denom
-                IF ( c_v < 0.0 )  THEN
-                   c_v = 0.0
-                ELSEIF ( c_v > c_max )  THEN
-                   c_v = c_max
-                ENDIF
-             ELSE
-                c_v = c_max
-             ENDIF
-
-             denom = w_m_s(k,0,i) - w_m_s(k,1,i)
-
-             IF ( denom /= 0.0 )  THEN
-                c_w = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
-                IF ( c_w < 0.0 )  THEN
-                   c_w = 0.0
-                ELSEIF ( c_w > c_max )  THEN
-                   c_w = c_max
-                ENDIF
-             ELSE
-                c_w = c_max
-             ENDIF
 
 !
 !--          Calculate the new velocities
-             u_p(k,-1,i) = u(k,-1,i) + dt_3d * c_u * &
+             u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u(k,i) * &
                                        ( u(k,-1,i) - u(k,0,i) ) * ddy
 
-             v_p(k,-1,i) = v(k,-1,i) + dt_3d * c_v * &
-                                       ( v(k,-1,i) - v_m_s(k,0,i) ) * ddy
-
-             w_p(k,-1,i) = w(k,-1,i) + dt_3d * c_w * &
+             v_p(k,-1,i) = v(k,-1,i) - dt_3d * tsc(2) * c_v(k,i) * &
+                                       ( 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(k,-1,i) - w(k,0,i) ) * ddy
-
-!
-!--          Save old timelevels for the next timestep
-             u_m_s(k,:,i) = u(k,-1:1,i)
-             v_m_s(k,:,i) = v(k,-1:1,i)
-             w_m_s(k,:,i) = w(k,-1:1,i)
 
           ENDDO
@@ -289,10 +298,10 @@
        IF ( ibc_uv_b == 0 )  THEN
           u_p(nzb,-1,:) = -u_p(nzb+1,-1,:) 
-          v_p(nzb,-1,:) = -v_p(nzb+1,-1,:)  
+          v_p(nzb,0,:)  = -v_p(nzb+1,0,:)  
        ELSE                    
           u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
-          v_p(nzb,-1,:) =  v_p(nzb+1,-1,:)
-       ENDIF
-       w_p(nzb,ny+1,:) = 0.0
+          v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
+       ENDIF
+       w_p(nzb,-1,:) = 0.0
 
 !
@@ -300,8 +309,8 @@
        IF ( ibc_uv_t == 0 )  THEN
           u_p(nzt+1,-1,:) = ug(nzt+1)
-          v_p(nzt+1,-1,:) = vg(nzt+1)
+          v_p(nzt+1,0,:)  = vg(nzt+1)
        ELSE
           u_p(nzt+1,-1,:) = u(nzt,-1,:)
-          v_p(nzt+1,-1,:) = v(nzt,-1,:)
+          v_p(nzt+1,0,:)  = v(nzt,0,:)
        ENDIF
        w_p(nzt:nzt+1,-1,:) = 0.0
@@ -309,7 +318,5 @@
     ENDIF
 
-    IF ( outflow_n  .AND.                                                 &
-         intermediate_timestep_count == intermediate_timestep_count_max ) &
-    THEN
+    IF ( outflow_n )  THEN
 
        c_max = dy / dt_3d
@@ -319,60 +326,66 @@
 
 !
-!--          First calculate the phase speeds for u,v, and w
-             denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
-
-             IF ( denom /= 0.0 )  THEN
-                c_u = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
-                IF ( c_u < 0.0 )  THEN
-                   c_u = 0.0
-                ELSEIF ( c_u > c_max )  THEN
-                   c_u = c_max
-                ENDIF
-             ELSE
-                c_u = c_max
+!--          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
+
+                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
+                   IF ( c_u(k,i) < 0.0 )  THEN
+                      c_u(k,i) = 0.0
+                   ELSEIF ( c_u(k,i) > c_max )  THEN
+                      c_u(k,i) = c_max
+                   ENDIF
+                ELSE
+                   c_u(k,i) = c_max
+                ENDIF
+
+                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
+                   IF ( c_v(k,i) < 0.0 )  THEN
+                      c_v(k,i) = 0.0
+                   ELSEIF ( c_v(k,i) > c_max )  THEN
+                      c_v(k,i) = c_max
+                   ENDIF
+                ELSE
+                   c_v(k,i) = c_max
+                ENDIF
+
+                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
+                   IF ( c_w(k,i) < 0.0 )  THEN
+                      c_w(k,i) = 0.0
+                   ELSEIF ( c_w(k,i) > c_max )  THEN
+                      c_w(k,i) = c_max
+                   ENDIF
+                ELSE
+                   c_w(k,i) = c_max
+                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
 
-             denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
-
-             IF ( denom /= 0.0 )  THEN
-                c_v = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
-                IF ( c_v < 0.0 )  THEN
-                   c_v = 0.0
-                ELSEIF ( c_v > c_max )  THEN
-                   c_v = c_max
-                ENDIF
-             ELSE
-                c_v = c_max
-             ENDIF
-
-             denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
-
-             IF ( denom /= 0.0 )  THEN
-                c_w = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
-                IF ( c_w < 0.0 )  THEN
-                   c_w = 0.0
-                ELSEIF ( c_w > c_max )  THEN
-                   c_w = c_max
-                ENDIF
-             ELSE
-                c_w = c_max
-             ENDIF
-
 !
 !--          Calculate the new velocities
-             u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * c_u * &
+             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 * c_v * &
+             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 * c_w * &
+             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
-
-!
-!--          Swap timelevels for the next timestep
-             u_m_n(k,:,i) = u(k,ny-1:ny+1,i)
-             v_m_n(k,:,i) = v(k,ny-1:ny+1,i)
-             w_m_n(k,:,i) = w(k,ny-1:ny+1,i)
 
           ENDDO
@@ -403,7 +416,5 @@
     ENDIF
 
-    IF ( outflow_l  .AND.                                                 &
-         intermediate_timestep_count == intermediate_timestep_count_max ) &
-    THEN
+    IF ( outflow_l )  THEN
 
        c_max = dx / dt_3d
@@ -413,60 +424,66 @@
 
 !
-!--          First calculate the phase speeds for u,v, and w
-             denom = u_m_l(k,j,0) - u_m_l(k,j,1)
-
-             IF ( denom /= 0.0 )  THEN
-                c_u = -c_max * ( u(k,j,0) - u_m_r(k,j,0) ) / denom
-                IF ( c_u > 0.0 )  THEN
-                   c_u = 0.0
-                ELSEIF ( c_u < -c_max )  THEN
-                   c_u = -c_max
-                ENDIF
-             ELSE
-                c_u = -c_max
+!--          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
+
+                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
+                   IF ( c_u(k,j) > 0.0 )  THEN
+                      c_u(k,j) = 0.0
+                   ELSEIF ( c_u(k,j) < -c_max )  THEN
+                      c_u(k,j) = -c_max
+                   ENDIF
+                ELSE
+                   c_u(k,j) = -c_max
+                ENDIF
+
+                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
+                   IF ( c_v(k,j) < 0.0 )  THEN
+                      c_v(k,j) = 0.0
+                   ELSEIF ( c_v(k,j) > c_max )  THEN
+                      c_v(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_v(k,j) = c_max
+                ENDIF
+
+                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
+                   IF ( c_w(k,j) < 0.0 )  THEN
+                      c_w(k,j) = 0.0
+                   ELSEIF ( c_w(k,j) > c_max )  THEN
+                      c_w(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_w(k,j) = c_max
+                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
 
-             denom = v_m_l(k,j,0) - v_m_l(k,j,1)
-
-             IF ( denom /= 0.0 )  THEN
-                c_v = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
-                IF ( c_v < 0.0 )  THEN
-                   c_v = 0.0
-                ELSEIF ( c_v > c_max )  THEN
-                   c_v = c_max
-                ENDIF
-             ELSE
-                c_v = c_max
-             ENDIF
-
-             denom = w_m_l(k,j,0) - w_m_l(k,j,1)
-
-             IF ( denom /= 0.0 )  THEN
-                c_w = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
-                IF ( c_w < 0.0 )  THEN
-                   c_w = 0.0
-                ELSEIF ( c_w > c_max )  THEN
-                   c_w = c_max
-                ENDIF
-             ELSE
-                c_w = c_max
-             ENDIF
-
 !
 !--          Calculate the new velocities
-             u_p(k,j,-1) = u(k,j,-1) + dt_3d * c_u * &
-                                       ( u(k,j,-1) - u(k,j,0) ) * ddx
-
-             v_p(k,j,-1) = v(k,j,-1) + dt_3d * c_v * &
+             u_p(k,j,0)  = u(k,j,0)  - dt_3d * tsc(2) * c_u(k,j) * &
+                                       ( 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(k,j,-1) - v(k,j,0) ) * ddx
 
-             w_p(k,j,-1) = w(k,j,-1) + dt_3d * c_w * &
+             w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w(k,j) * &
                                        ( w(k,j,-1) - w(k,j,0) ) * ddx
-
-!
-!--          Swap timelevels for the next timestep
-             u_m_l(k,j,:) = u(k,j,-1:1)
-             v_m_l(k,j,:) = v(k,j,-1:1)
-             w_m_l(k,j,:) = w(k,j,-1:1)
 
           ENDDO
@@ -497,7 +514,5 @@
     ENDIF
 
-    IF ( outflow_r  .AND.                                                 &
-         intermediate_timestep_count == intermediate_timestep_count_max ) &
-    THEN
+    IF ( outflow_r )  THEN
 
        c_max = dx / dt_3d
@@ -507,60 +522,66 @@
 
 !
-!--          First calculate the phase speeds for u,v, and w
-             denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
-
-             IF ( denom /= 0.0 )  THEN
-                c_u = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
-                IF ( c_u < 0.0 )  THEN
-                   c_u = 0.0
-                ELSEIF ( c_u > c_max )  THEN
-                   c_u = c_max
-                ENDIF
-             ELSE
-                c_u = c_max
+!--          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
+
+                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
+                   IF ( c_u(k,j) < 0.0 )  THEN
+                      c_u(k,j) = 0.0
+                   ELSEIF ( c_u(k,j) > c_max )  THEN
+                      c_u(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_u(k,j) = c_max
+                ENDIF
+
+                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
+                   IF ( c_v(k,j) < 0.0 )  THEN
+                      c_v(k,j) = 0.0
+                   ELSEIF ( c_v(k,j) > c_max )  THEN
+                      c_v(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_v(k,j) = c_max
+                ENDIF
+
+                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
+                   IF ( c_w(k,j) < 0.0 )  THEN
+                      c_w(k,j) = 0.0
+                   ELSEIF ( c_w(k,j) > c_max )  THEN
+                      c_w(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_w(k,j) = c_max
+                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
 
-             denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
-
-             IF ( denom /= 0.0 )  THEN
-                c_v = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
-                IF ( c_v < 0.0 )  THEN
-                   c_v = 0.0
-                ELSEIF ( c_v > c_max )  THEN
-                   c_v = c_max
-                ENDIF
-             ELSE
-                c_v = c_max
-             ENDIF
-
-             denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
-
-             IF ( denom /= 0.0 )  THEN
-                c_w = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
-                IF ( c_w < 0.0 )  THEN
-                   c_w = 0.0
-                ELSEIF ( c_w > c_max )  THEN
-                   c_w = c_max
-                ENDIF
-             ELSE
-                c_w = c_max
-             ENDIF
-
 !
 !--          Calculate the new velocities
-             u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * c_u * &
+             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 * c_v * &
+             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 * c_w * &
+             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
-
-!
-!--          Swap timelevels for the next timestep
-             u_m_r(k,j,:) = u(k,j,nx-1:nx+1)
-             v_m_r(k,j,:) = v(k,j,nx-1:nx+1)
-             w_m_r(k,j,:) = w(k,j,nx-1:nx+1)
 
           ENDDO
Index: palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- palm/trunk/SOURCE/buoyancy.f90	(revision 105)
+++ palm/trunk/SOURCE/buoyancy.f90	(revision 106)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! i loop for u-component (sloping surface) is starting from nxlu (needed for
+! non-cyclic boundary conditions)
 !
 ! Former revisions:
@@ -105,5 +106,5 @@
           IF ( wind_component == 1 )  THEN
 
-             DO  i = nxl, nxr
+             DO  i = nxlu, nxr
                 DO  j = nys, nyn
                    DO  k = nzb_s_inner(j,i)+1, nzt-1
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 105)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 106)
@@ -5,5 +5,6 @@
 ! -----------------
 ! Check coupling_mode and set default (obligatory) values (like boundary
-! conditions for temperature and fluxes) in case of coupled runs
+! conditions for temperature and fluxes) in case of coupled runs.
+! +profiles for w*p* and w"e
 ! Bugfix: Error message concerning output of particle concentration (pc)
 ! modified
@@ -1968,5 +1969,5 @@
              dopr_index(i) = 56
              dopr_unit(i)  = 'm2/s3'
-             hom(:,2,56,:) = SPREAD( zu, 2, statistic_regions+1 )
+             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'w"e/dz' )
@@ -2051,4 +2052,14 @@
                 hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
              ENDIF
+
+          CASE ( 'w*p*' )
+             dopr_index(i) = 68
+             dopr_unit(i)  = 'm3/s3'
+             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'w"e' )
+             dopr_index(i) = 69
+             dopr_unit(i)  = 'm3/s3'
+             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
 
 
Index: palm/trunk/SOURCE/coriolis.f90
===================================================================
--- palm/trunk/SOURCE/coriolis.f90	(revision 105)
+++ palm/trunk/SOURCE/coriolis.f90	(revision 106)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! loops for u and v are starting from index nxlu, nysv, respectively (needed
+! for non-cyclic boundary conditions)
 !
 ! Former revisions:
@@ -60,5 +61,5 @@
 !--       u-component
           CASE ( 1 )
-             DO  i = nxl, nxr
+             DO  i = nxlu, nxr
                 DO  j = nys, nyn
                    DO  k = nzb_u_inner(j,i)+1, nzt
@@ -78,5 +79,5 @@
           CASE ( 2 )
              DO  i = nxl, nxr
-                DO  j = nys, nyn
+                DO  j = nysv, nyn
                    DO  k = nzb_v_inner(j,i)+1, nzt
                       tend(k,j,i) = tend(k,j,i) - f *     ( 0.25 *          &
Index: palm/trunk/SOURCE/data_output_dvrp.f90
===================================================================
--- palm/trunk/SOURCE/data_output_dvrp.f90	(revision 105)
+++ palm/trunk/SOURCE/data_output_dvrp.f90	(revision 106)
@@ -151,5 +151,6 @@
 !
 !--       Definition of characteristics of particle material
-          CALL DVRP_MATERIAL_RGB( m-1, 1, 0.1, 0.7, 0.1, 0.0 )
+!          CALL DVRP_MATERIAL_RGB( m-1, 1, 0.1, 0.7, 0.1, 0.0 )
+          CALL DVRP_MATERIAL_RGB( m-1, 1, 0.0, 0.0, 0.0, 0.0 )
 
 !
@@ -318,13 +319,14 @@
                    ENDDO
                 ENDDO
-                DO  k = nzb, nz_do3d
-                   DO  j = nys+1, nyn
-                      DO  i = nxl, nxr+1
-                         local_pf(i,j,k) = 0.25 * local_pf(i,j-1,k) + &
-                                           0.50 * local_pf(i,j,k)   + &
-                                           0.25 * local_pf(i,j+1,k)
-                      ENDDO
-                   ENDDO
-                ENDDO
+! Averaging for Langmuir circulation
+!                DO  k = nzb, nz_do3d
+!                   DO  j = nys+1, nyn
+!                      DO  i = nxl, nxr+1
+!                         local_pf(i,j,k) = 0.25 * local_pf(i,j-1,k) + &
+!                                           0.50 * local_pf(i,j,k)   + &
+!                                           0.25 * local_pf(i,j+1,k)
+!                      ENDDO
+!                   ENDDO
+!                ENDDO
 
              CASE ( 'p', 'p_xy', 'p_xz', 'p_yz' )
@@ -475,5 +477,13 @@
              CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
                              cyclic_dvrp, cyclic_dvrp, cyclic_dvrp )
-             CALL DVRP_SLICER( m-1, section_mode, slicer_position )
+!             CALL DVRP_SLICER( m-1, section_mode, slicer_position )
+             CALL DVRP_SLICER( m-1, 2, 1.0 )
+             WRITE (9,*) 'nx_dvrp=', nx_dvrp
+             WRITE (9,*) 'ny_dvrp=', ny_dvrp
+             WRITE (9,*) 'nz_dvrp=', nz_dvrp
+             WRITE (9,*) 'section_mode=', section_mode
+             WRITE (9,*) 'slicer_position=', slicer_position
+             CALL local_flush( 9 )
+
              CALL DVRP_VISUALIZE( m-1, 2, dvrp_filecount )
 
Index: palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_u.f90	(revision 105)
+++ palm/trunk/SOURCE/diffusion_u.f90	(revision 106)
@@ -4,6 +4,7 @@
 ! Actual revisions:
 ! -----------------
-! Momentumflux at top (uswst) included as boundary condition
-! 
+! Momentumflux at top (uswst) included as boundary condition,
+! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
+!
 ! Former revisions:
 ! -----------------
@@ -79,5 +80,5 @@
        ENDIF
 
-       DO  i = nxl, nxr
+       DO  i = nxlu, nxr
           DO  j = nys,nyn
 !
Index: palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_v.f90	(revision 105)
+++ palm/trunk/SOURCE/diffusion_v.f90	(revision 106)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! Momentumflux at top (vswst) included as boundary condition
+! Momentumflux at top (vswst) included as boundary condition,
+! j loop is starting from nysv (needed for non-cyclic boundary conditions)
 !
 ! Former revisions:
@@ -78,5 +79,5 @@
 
        DO  i = nxl, nxr
-          DO  j = nys, nyn
+          DO  j = nysv, nyn
 !
 !--          Compute horizontal diffusion
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 105)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 106)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! Prescribed momentum fluxes at the top surface are used
+! Prescribed momentum fluxes at the top surface are used,
+! profiles for w*p* and w"e are calculated
 !
 ! Former revisions:
@@ -620,7 +621,8 @@
 !
 !--    Divergence of vertical flux of resolved scale energy and pressure
-!--    fluctuations. First calculate the products, then the divergence.
+!--    fluctuations as well as flux of pressure fluctuation itself (68).
+!--    First calculate the products, then the divergence.
 !--    Calculation is time consuming. Do it only, if profiles shall be plotted.
-       IF ( hom(nzb+1,2,55,0) /= 0.0 )  THEN
+       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
 
           sums_ll = 0.0  ! local array
@@ -654,13 +656,15 @@
              sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
              sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
+             sums_l(k,68,tn) = sums_ll(k,2)
           ENDDO
           sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
           sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
-
-       ENDIF
-
-!
-!--    Divergence of vertical flux of SGS TKE
-       IF ( hom(nzb+1,2,57,0) /= 0.0 )  THEN
+          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
+
+       ENDIF
+
+!
+!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
+       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
 
           !$OMP DO
@@ -669,8 +673,12 @@
                 DO  k = nzb_s_outer(j,i)+1, nzt
 
-                   sums_l(k,57,tn) = sums_l(k,57,tn) + (                       &
+                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
                  - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
-                                                 ) * ddzw(k)
+                                                             ) * ddzw(k)
+
+                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
+                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
+                                                              )
 
                 ENDDO
@@ -678,4 +686,5 @@
           ENDDO
           sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
+          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
 
        ENDIF
@@ -832,5 +841,5 @@
        hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
        hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
-       hom(:,1,57,sr) = sums(:,57)     ! w"e/dz
+       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
        hom(:,1,58,sr) = sums(:,58)     ! u"pt"
        hom(:,1,59,sr) = sums(:,59)     ! u*pt*
@@ -843,4 +852,6 @@
        hom(:,1,66,sr) = sums(:,66)     ! w*sa*
        hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
+       hom(:,1,68,sr) = sums(:,68)     ! w*p*
+       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
 
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 105)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 106)
@@ -7,5 +7,6 @@
 ! Actual revisions:
 ! -----------------
-! Flux initialization in case of coupled runs, +momentum fluxes at top boundary
+! Flux initialization in case of coupled runs, +momentum fluxes at top boundary,
+! +arrays for phase speed c_u, c_v, c_w, indices for u|v|w_m_l|r changed
 !
 ! Former revisions:
@@ -233,25 +234,33 @@
 
 !
-!-- Arrays to store velocity data from t-dt needed for radiation boundary
-!-- conditions
+!-- Arrays to store velocity data from t-dt and the phase speeds which
+!-- are needed for radiation boundary conditions
     IF ( outflow_l )  THEN
-       ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1), &
-                 v_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1), &
-                 w_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1) )
+       ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,1:2), &
+                 v_m_l(nzb:nzt+1,nys-1:nyn+1,0:1), &
+                 w_m_l(nzb:nzt+1,nys-1:nyn+1,0:1) )
     ENDIF
     IF ( outflow_r )  THEN
-       ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1), &
-                 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1), &
-                 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1) )
+       ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
+                 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
+                 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx) )
+    ENDIF
+    IF ( outflow_l  .OR.  outflow_r )  THEN
+       ALLOCATE( c_u(nzb:nzt+1,nys-1:nyn+1), c_v(nzb:nzt+1,nys-1:nyn+1), &
+                 c_w(nzb:nzt+1,nys-1:nyn+1) )
     ENDIF
     IF ( outflow_s )  THEN
-       ALLOCATE( u_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1), &
-                 v_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1), &
-                 w_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1) )
+       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1), &
+                 v_m_s(nzb:nzt+1,1:2,nxl-1:nxr+1), &
+                 w_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1) )
     ENDIF
     IF ( outflow_n )  THEN
-       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1), &
-                 v_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1), &
-                 w_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1) )
+       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
+                 v_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
+                 w_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1) )
+    ENDIF
+    IF ( outflow_s  .OR.  outflow_n )  THEN
+       ALLOCATE( c_u(nzb:nzt+1,nxl-1:nxr+1), c_v(nzb:nzt+1,nxl-1:nxr+1), &
+                 c_w(nzb:nzt+1,nxl-1:nxr+1) )
     ENDIF
 
@@ -808,22 +817,22 @@
 !--    Initialize old timelevels needed for radiation boundary conditions
        IF ( outflow_l )  THEN
-          u_m_l(:,:,:) = u(:,:,-1:1)
-          v_m_l(:,:,:) = v(:,:,-1:1)
-          w_m_l(:,:,:) = w(:,:,-1:1)
+          u_m_l(:,:,:) = u(:,:,1:2)
+          v_m_l(:,:,:) = v(:,:,0:1)
+          w_m_l(:,:,:) = w(:,:,0:1)
        ENDIF
        IF ( outflow_r )  THEN
-          u_m_r(:,:,:) = u(:,:,nx-1:nx+1)
-          v_m_r(:,:,:) = v(:,:,nx-1:nx+1)
-          w_m_r(:,:,:) = w(:,:,nx-1:nx+1)
+          u_m_r(:,:,:) = u(:,:,nx-1:nx)
+          v_m_r(:,:,:) = v(:,:,nx-1:nx)
+          w_m_r(:,:,:) = w(:,:,nx-1:nx)
        ENDIF
        IF ( outflow_s )  THEN
-          u_m_s(:,:,:) = u(:,-1:1,:)
-          v_m_s(:,:,:) = v(:,-1:1,:)
-          w_m_s(:,:,:) = w(:,-1:1,:)
+          u_m_s(:,:,:) = u(:,0:1,:)
+          v_m_s(:,:,:) = v(:,1:2,:)
+          w_m_s(:,:,:) = w(:,0:1,:)
        ENDIF
        IF ( outflow_n )  THEN
-          u_m_n(:,:,:) = u(:,ny-1:ny+1,:)
-          v_m_n(:,:,:) = v(:,ny-1:ny+1,:)
-          w_m_n(:,:,:) = w(:,ny-1:ny+1,:)
+          u_m_n(:,:,:) = u(:,ny-1:ny,:)
+          v_m_n(:,:,:) = v(:,ny-1:ny,:)
+          w_m_n(:,:,:) = w(:,ny-1:ny,:)
        ENDIF
 
Index: palm/trunk/SOURCE/init_particles.f90
===================================================================
--- palm/trunk/SOURCE/init_particles.f90	(revision 105)
+++ palm/trunk/SOURCE/init_particles.f90	(revision 106)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! variable iran replaced by iran_part
 !
 ! Former revisions:
@@ -354,11 +354,8 @@
        IF ( random_start_position )  THEN
 
-          iran = iran + myid    ! Random positions should be different on
-                                ! different PEs
-
           DO  n = 1, number_of_initial_particles
              IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
                 particles(n)%x = particles(n)%x + &
-                                 ( random_function( iran ) - 0.5 ) * &
+                                 ( random_function( iran_part ) - 0.5 ) * &
                                  pdx(particles(n)%group)
                 IF ( particles(n)%x  <=  ( nxl - 0.5 ) * dx )  THEN
@@ -370,5 +367,5 @@
              IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
                 particles(n)%y = particles(n)%y + &
-                                 ( random_function( iran ) - 0.5 ) * &
+                                 ( random_function( iran_part ) - 0.5 ) * &
                                  pdy(particles(n)%group)
                 IF ( particles(n)%y  <=  ( nys - 0.5 ) * dy )  THEN
@@ -380,5 +377,5 @@
              IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
                 particles(n)%z = particles(n)%z + &
-                                 ( random_function( iran ) - 0.5 ) * &
+                                 ( random_function( iran_part ) - 0.5 ) * &
                                  pdz(particles(n)%group)
              ENDIF
Index: palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- palm/trunk/SOURCE/init_pegrid.f90	(revision 105)
+++ palm/trunk/SOURCE/init_pegrid.f90	(revision 106)
@@ -5,5 +5,7 @@
 ! -----------------
 ! Intercommunicator (comm_inter) and derived data type (type_xy) for
-! coupled model runs created
+! coupled model runs created,
+! indices nxlu and nysv are calculated (needed for non-cyclic boundary
+! conditions)
 !
 ! Former revisions:
@@ -818,6 +820,5 @@
 !
 !-- Setting of flags for inflow/outflow conditions in case of non-cyclic
-!-- horizontal boundary conditions. Set variables for extending array u (v)
-!-- by one gridpoint on the left/rightmost (northest/southest) processor
+!-- horizontal boundary conditions.
     IF ( pleft == MPI_PROC_NULL )  THEN
        IF ( bc_lr == 'dirichlet/radiation' )  THEN
@@ -869,4 +870,17 @@
     ENDIF
 #endif
+!
+!-- At the outflow, u or v, respectively, have to be calculated for one grid
+!-- point less.
+    IF ( outflow_l )  THEN
+       nxlu = nxl + 1
+    ELSE
+       nxlu = nxl
+    ENDIF
+    IF ( outflow_s )  THEN
+       nysv = nys + 1
+    ELSE
+       nysv = nys
+    ENDIF
 
     IF ( psolver == 'poisfft_hybrid' )  THEN
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 105)
+++ palm/trunk/SOURCE/modules.f90	(revision 106)
@@ -6,6 +6,6 @@
 ! -----------------
 ! +comm_inter, constant_top_momentumflux, coupling_char, coupling_mode,
-! dt_coupling, ngp_xy, port_name, time_coupling, top_momentumflux_u|v,
-! type_xy, uswst*, vswst*
+! c_u, c_v, c_w, dt_coupling, ngp_xy, nxlu, nysv, port_name, time_coupling,
+! top_momentumflux_u|v, type_xy, uswst*, vswst*
 !
 ! Former revisions:
@@ -105,5 +105,6 @@
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
-          dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg, pt_slope_ref, qs, ts, us, z0
+          c_u, c_v, c_w, dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg, pt_slope_ref,    &
+          qs, ts, us, z0
 
     REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
@@ -565,7 +566,7 @@
 !------------------------------------------------------------------------------!
 
-    INTEGER ::  ngp_sums, nnx, nx = 0, nxa, nxl, nxr, nxra, nny, ny = 0, nya,  &
-                nyn, nyna, nys, nnz, nz = 0, nza, nzb, nzb_diff, nzt, nzta,    &
-                nzt_diff
+    INTEGER ::  ngp_sums, nnx, nx = 0, nxa, nxl, nxlu, nxr, nxra, nny, ny = 0, &
+                nya, nyn, nyna, nys, nysv, nnz, nz = 0, nza, nzb, nzb_diff,    &
+                nzt, nzta, nzt_diff
 
     INTEGER, DIMENSION(:), ALLOCATABLE ::                                      &
Index: palm/trunk/SOURCE/pres.f90
===================================================================
--- palm/trunk/SOURCE/pres.f90	(revision 105)
+++ palm/trunk/SOURCE/pres.f90	(revision 106)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Volume flow conservation added for the remaining three outflow boundaries
 !
 ! Former revisions:
@@ -74,19 +74,68 @@
 !-- Conserve the volume flow at the outflow in case of non-cyclic lateral
 !-- boundary conditions
-    IF ( conserve_volume_flow  .AND.  bc_ns == 'radiation/dirichlet')  THEN
+!-- WARNING: so far, this conservation does not work at the left/south
+!--          boundary if the topography at the inflow differs from that at the
+!--          outflow! For this case, volume_flow_area needs adjustment!
+!
+!-- Left/right
+    IF ( conserve_volume_flow  .AND.  ( outflow_l  .OR. outflow_r ) )  THEN
+
+       volume_flow(1)   = 0.0
+       volume_flow_l(1) = 0.0
+
+       IF ( outflow_l )  THEN
+          i = 0
+       ELSEIF ( outflow_r )  THEN
+          i = nx+1
+       ENDIF
+
+       DO  j = nys, nyn
+!
+!--       Sum up the volume flow through the south/north boundary
+          DO  k = nzb_2d(j,i) + 1, nzt
+             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzu(k)
+          ENDDO
+       ENDDO
+
+#if defined( __parallel )   
+       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
+                           MPI_SUM, comm1dy, ierr )   
+#else
+       volume_flow = volume_flow_l  
+#endif
+       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) )    &
+                               / volume_flow_area(1)
+
+       DO  j = nys, nyn
+          DO  k = nzb_v_inner(j,i) + 1, nzt
+             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
+          ENDDO
+       ENDDO
+
+       CALL exchange_horiz( u )
+
+    ENDIF
+
+!
+!-- South/north
+    IF ( conserve_volume_flow  .AND.  ( outflow_n  .OR. outflow_s ) )  THEN
 
        volume_flow(2)   = 0.0
        volume_flow_l(2) = 0.0
 
-       IF ( nyn == ny )  THEN
+       IF ( outflow_s )  THEN
+          j = 0
+       ELSEIF ( outflow_n )  THEN
           j = ny+1
-          DO  i = nxl, nxr
-!
-!--          Sum up the volume flow through the north boundary
-             DO  k = nzb_2d(j,i) + 1, nzt
-                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k)
-             ENDDO
-          ENDDO
-       ENDIF
+       ENDIF
+
+       DO  i = nxl, nxr
+!
+!--       Sum up the volume flow through the south/north boundary
+          DO  k = nzb_2d(j,i) + 1, nzt
+             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k)
+          ENDDO
+       ENDDO
+
 #if defined( __parallel )   
        CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
@@ -96,14 +145,11 @@
 #endif
        volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) )    &
-                               / volume_flow_area(2)                          
-
-       IF ( outflow_n )  THEN
-          j = nyn+1
-          DO  i = nxl, nxr
-             DO  k = nzb_v_inner(j,i) + 1, nzt
-                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
-             ENDDO
-          ENDDO
-       ENDIF
+                               / volume_flow_area(2)
+
+       DO  i = nxl, nxr
+          DO  k = nzb_v_inner(j,i) + 1, nzt
+             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
+          ENDDO
+       ENDDO
 
        CALL exchange_horiz( v )
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 105)
+++ palm/trunk/SOURCE/production_e.f90	(revision 106)
@@ -4,5 +4,8 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Bugfix: wrong sign removed from the buoyancy production term in the case
+! use_reference = .T.,
+! u_0 and v_0 are calculated for nxr+1, nyn+1 also (otherwise these values are
+! not available in case of non-cyclic boundary conditions)
 !
 ! Former revisions:
@@ -386,5 +389,5 @@
                    DO  j = nys, nyn
                       DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                         tend(k,j,i) = tend(k,j,i) +                  &
+                         tend(k,j,i) = tend(k,j,i) -                  &
                                        kh(k,j,i) * g / pt_reference * &
                                        ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
@@ -976,9 +979,11 @@
 !--       (otherwise the timestep may be significantly reduced by large
 !--       surface winds).
+!--       Upper bounds are nxr+1 and nyn+1 because otherwise these values are
+!--       not available in case of non-cyclic boundary conditions.
 !--       WARNING: the exact analytical solution would require the determination
 !--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
           !$OMP PARALLEL DO PRIVATE( ku, kv )
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
+          DO  i = nxl, nxr+1
+             DO  j = nys, nyn+1
 
                 ku = nzb_u_inner(j,i)+1
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 105)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 106)
@@ -4,5 +4,7 @@
 ! Actual revisions:
 ! -----------------
-! +uswst, vswst as arguments in calls of diffusion_u|v
+! +uswst, vswst as arguments in calls of diffusion_u|v,
+! loops for u and v are starting from index nxlu, nysv, respectively (needed
+! for non-cyclic boundary conditions)
 !
 ! Former revisions:
@@ -133,5 +135,5 @@
 !
 !-- u-tendency terms with no communication
-    DO  i = nxl, nxr
+    DO  i = nxlu, nxr
        DO  j = nys, nyn
 !
@@ -202,5 +204,5 @@
 !-- v-tendency terms with no communication
     DO  i = nxl, nxr
-       DO  j = nys, nyn
+       DO  j = nysv, nyn
 !
 !--       Tendency terms
@@ -806,5 +808,5 @@
 !
 !--       Tendency terms for u-velocity component
-          IF ( j < nyn+1 )  THEN
+          IF ( .NOT. outflow_l  .OR.  i > nxl )  THEN
 
              tend(:,j,i) = 0.0
@@ -857,5 +859,5 @@
 !
 !--       Tendency terms for v-velocity component
-          IF ( i < nxr+1 )  THEN
+          IF ( .NOT. outflow_s  .OR.  j > nys )  THEN
 
              tend(:,j,i) = 0.0
@@ -906,96 +908,139 @@
 !
 !--       Tendency terms for w-velocity component
-          IF ( i < nxr+1  .AND.  j < nyn+1 )  THEN
-
+          tend(:,j,i) = 0.0
+          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
+             CALL advec_w_pw( i, j )
+          ELSE
+             CALL advec_w_up( i, j )
+          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 )
+          ELSE
+             CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
+                               tend, u, v, w )
+          ENDIF
+          CALL coriolis( i, j, 3 )
+          IF ( ocean )  THEN
+             CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
+          ELSE
+             IF ( .NOT. humidity )  THEN
+                CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
+             ELSE
+                CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
+             ENDIF
+          ENDIF
+          CALL user_actions( i, j, 'w-tendency' )
+
+!
+!--       Prognostic equation for w-velocity component
+          DO  k = nzb_w_inner(j,i)+1, nzt-1
+             w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
+                          dt_3d * (                                         &
+                                tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
+                           - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
+                                  ) -                                       &
+                          tsc(5) * rdf(k) * w(k,j,i)
+          ENDDO
+
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  k = nzb_w_inner(j,i)+1, nzt-1
+                   tw_m(k,j,i) = tend(k,j,i)
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  k = nzb_w_inner(j,i)+1, nzt-1
+                   tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
+                ENDDO
+             ENDIF
+          ENDIF
+
+!
+!--       Tendency terms for potential temperature
+          tend(:,j,i) = 0.0
+          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
+             CALL advec_s_pw( i, j, pt )
+          ELSE
+             CALL advec_s_up( i, j, pt )
+          ENDIF
+          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
+          THEN
+             CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
+                               tswst_m, tend )
+          ELSE
+             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
+          ENDIF
+
+!
+!--       If required compute heating/cooling due to long wave radiation
+!--       processes
+          IF ( radiation )  THEN
+             CALL calc_radiation( i, j )
+          ENDIF
+
+!
+!--       If required compute impact of latent heat due to precipitation
+          IF ( precipitation )  THEN
+             CALL impact_of_latent_heat( i, j )
+          ENDIF
+          CALL user_actions( i, j, 'pt-tendency' )
+
+!
+!--       Prognostic equation for potential temperature
+          DO  k = nzb_s_inner(j,i)+1, nzt
+             pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
+                           dt_3d * (                                        &
+                               tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
+                                   ) -                                      &
+                           tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
+          ENDDO
+
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   tpt_m(k,j,i) = tend(k,j,i)
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                   5.3125 * tpt_m(k,j,i)
+                ENDDO
+             ENDIF
+          ENDIF
+
+!
+!--       If required, compute prognostic equation for salinity
+          IF ( ocean )  THEN
+
+!
+!--          Tendency-terms for salinity
              tend(:,j,i) = 0.0
-             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
-                CALL advec_w_pw( i, j )
-             ELSE
-                CALL advec_w_up( i, j )
-             ENDIF
-             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
+             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
              THEN
-                CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x,          &
-                                  km_damp_y, 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 )
-             ENDIF
-             CALL coriolis( i, j, 3 )
-             IF ( ocean )  THEN
-                CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
-             ELSE
-                IF ( .NOT. humidity )  THEN
-                   CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
-                ELSE
-                   CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
-                ENDIF
-             ENDIF
-             CALL user_actions( i, j, 'w-tendency' )
-
-!
-!--          Prognostic equation for w-velocity component
-             DO  k = nzb_w_inner(j,i)+1, nzt-1
-                w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
-                             dt_3d * (                                         &
-                                   tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
-                              - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
-                                     ) -                                       &
-                             tsc(5) * rdf(k) * w(k,j,i)
-             ENDDO
-
-!
-!--          Calculate tendencies for the next Runge-Kutta step
-             IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                IF ( intermediate_timestep_count == 1 )  THEN
-                   DO  k = nzb_w_inner(j,i)+1, nzt-1
-                      tw_m(k,j,i) = tend(k,j,i)
-                   ENDDO
-                ELSEIF ( intermediate_timestep_count < &
-                         intermediate_timestep_count_max )  THEN
-                   DO  k = nzb_w_inner(j,i)+1, nzt-1
-                      tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
-                   ENDDO
-                ENDIF
-             ENDIF
-
-!
-!--          Tendency terms for potential temperature
-             tend(:,j,i) = 0.0
-             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
-                CALL advec_s_pw( i, j, pt )
-             ELSE
-                CALL advec_s_up( i, j, pt )
-             ENDIF
-             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
-             THEN
-                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
-                                  tswst_m, tend )
-             ELSE
-                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
-             ENDIF
-
-!
-!--          If required compute heating/cooling due to long wave radiation
-!--          processes
-             IF ( radiation )  THEN
-                CALL calc_radiation( i, j )
-             ENDIF
-
-!
-!--          If required compute impact of latent heat due to precipitation
-             IF ( precipitation )  THEN
-                CALL impact_of_latent_heat( i, j )
-             ENDIF
-             CALL user_actions( i, j, 'pt-tendency' )
-
-!
-!--          Prognostic equation for potential temperature
+                CALL advec_s_pw( i, j, sa )
+             ELSE
+                CALL advec_s_up( i, j, sa )
+             ENDIF
+             CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
+                               tend )
+
+             CALL user_actions( i, j, 'sa-tendency' )
+
+!
+!--          Prognostic equation for salinity
              DO  k = nzb_s_inner(j,i)+1, nzt
-                pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
-                              dt_3d * (                                        &
-                                  tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
-                                      ) -                                      &
-                              tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
+                sa_p(k,j,i) = tsc(1) * sa(k,j,i) +                          &
+                              dt_3d * (                                     &
+                               tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
+                                      ) -                                   &
+                             tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
+                IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
              ENDDO
 
@@ -1005,203 +1050,156 @@
                 IF ( intermediate_timestep_count == 1 )  THEN
                    DO  k = nzb_s_inner(j,i)+1, nzt
-                      tpt_m(k,j,i) = tend(k,j,i)
+                      tsa_m(k,j,i) = tend(k,j,i)
                    ENDDO
                 ELSEIF ( intermediate_timestep_count < &
                          intermediate_timestep_count_max )  THEN
                    DO  k = nzb_s_inner(j,i)+1, nzt
-                      tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
-                                      5.3125 * tpt_m(k,j,i)
-                   ENDDO
-                ENDIF
-             ENDIF
-
-!
-!--          If required, compute prognostic equation for salinity
-             IF ( ocean )  THEN
-
-!
-!--             Tendency-terms for salinity
-                tend(:,j,i) = 0.0
-                IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
-                THEN
-                   CALL advec_s_pw( i, j, sa )
+                      tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                      5.3125 * tsa_m(k,j,i)
+                   ENDDO
+                ENDIF
+             ENDIF
+
+!
+!--          Calculate density by the equation of state for seawater
+             CALL eqn_state_seawater( i, j )
+
+          ENDIF
+
+!
+!--       If required, compute prognostic equation for total water content /
+!--       scalar
+          IF ( humidity  .OR.  passive_scalar )  THEN
+
+!
+!--          Tendency-terms for total water content / scalar
+             tend(:,j,i) = 0.0
+             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
+             THEN
+                CALL advec_s_pw( i, j, q )
+             ELSE
+                CALL advec_s_up( i, j, q )
+             ENDIF
+             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
+             THEN
+                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
+                                  qswst_m, tend )
+             ELSE
+                CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
+                                  tend )
+             ENDIF
+       
+!
+!--          If required compute decrease of total water content due to
+!--          precipitation
+             IF ( precipitation )  THEN
+                CALL calc_precipitation( i, j )
+             ENDIF
+             CALL user_actions( i, j, 'q-tendency' )
+
+!
+!--          Prognostic equation for total water content / scalar
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
+                             dt_3d * (                                      &
+                                tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
+                                     ) -                                    &
+                             tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
+                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
+             ENDDO
+
+!
+!--          Calculate tendencies for the next Runge-Kutta step
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( intermediate_timestep_count == 1 )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tq_m(k,j,i) = tend(k,j,i)
+                   ENDDO
+                ELSEIF ( intermediate_timestep_count < &
+                         intermediate_timestep_count_max )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                     5.3125 * tq_m(k,j,i)
+                   ENDDO
+                ENDIF
+             ENDIF
+
+          ENDIF
+
+!
+!--       If required, compute prognostic equation for turbulent kinetic 
+!--       energy (TKE)
+          IF ( .NOT. constant_diffusion )  THEN
+
+!
+!--          Tendency-terms for TKE
+             tend(:,j,i) = 0.0
+             IF ( ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  &
+                  .AND.  .NOT. use_upstream_for_tke )  THEN
+                CALL advec_s_pw( i, j, e )
+             ELSE
+                CALL advec_s_up( i, j, e )
+             ENDIF
+             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
+             THEN
+                IF ( .NOT. humidity )  THEN
+                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
+                                     km_m, l_grid, pt_m, pt_reference,   &
+                                     rif_m, tend, zu, zw )
                 ELSE
-                   CALL advec_s_up( i, j, sa )
-                ENDIF
-                CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
-                                  tend )
-       
-                CALL user_actions( i, j, 'sa-tendency' )
-
-!
-!--             Prognostic equation for salinity
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   sa_p(k,j,i) = tsc(1) * sa(k,j,i) +                          &
-                                 dt_3d * (                                     &
-                                  tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
-                                         ) -                                   &
-                                tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
-                   IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
-                ENDDO
-
-!
-!--             Calculate tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      DO  k = nzb_s_inner(j,i)+1, nzt
-                         tsa_m(k,j,i) = tend(k,j,i)
-                      ENDDO
-                   ELSEIF ( intermediate_timestep_count < &
-                            intermediate_timestep_count_max )  THEN
-                      DO  k = nzb_s_inner(j,i)+1, nzt
-                         tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
-                                         5.3125 * tsa_m(k,j,i)
-                      ENDDO
-                   ENDIF
-                ENDIF
-
-!
-!--             Calculate density by the equation of state for seawater
-                CALL eqn_state_seawater( i, j )
-
-             ENDIF
-
-!
-!--          If required, compute prognostic equation for total water content /
-!--          scalar
-             IF ( humidity  .OR.  passive_scalar )  THEN
-
-!
-!--             Tendency-terms for total water content / scalar
-                tend(:,j,i) = 0.0
-                IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
-                THEN
-                   CALL advec_s_pw( i, j, q )
-                ELSE
-                   CALL advec_s_up( i, j, q )
-                ENDIF
-                IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
-                THEN
-                   CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
-                                     qswst_m, tend )
-                ELSE
-                   CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
-                                     tend )
-                ENDIF
-       
-!
-!--             If required compute decrease of total water content due to
-!--             precipitation
-                IF ( precipitation )  THEN
-                   CALL calc_precipitation( i, j )
-                ENDIF
-                CALL user_actions( i, j, 'q-tendency' )
-
-!
-!--             Prognostic equation for total water content / scalar
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
-                                dt_3d * (                                      &
-                                   tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
-                                        ) -                                    &
-                                tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
-                   IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
-                ENDDO
-
-!
-!--             Calculate tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      DO  k = nzb_s_inner(j,i)+1, nzt
-                         tq_m(k,j,i) = tend(k,j,i)
-                      ENDDO
-                   ELSEIF ( intermediate_timestep_count < &
-                            intermediate_timestep_count_max )  THEN
-                      DO  k = nzb_s_inner(j,i)+1, nzt
-                         tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
-                                        5.3125 * tq_m(k,j,i)
-                      ENDDO
-                   ENDIF
-                ENDIF
-
-             ENDIF
-
-!
-!--          If required, compute prognostic equation for turbulent kinetic 
-!--          energy (TKE)
-             IF ( .NOT. constant_diffusion )  THEN
-
-!
-!--             Tendency-terms for TKE
-                tend(:,j,i) = 0.0
-                IF ( ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  &
-                     .AND.  .NOT. use_upstream_for_tke )  THEN
-                   CALL advec_s_pw( i, j, e )
-                ELSE
-                   CALL advec_s_up( i, j, e )
-                ENDIF
-                IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
-                THEN
-                   IF ( .NOT. humidity )  THEN
-                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
-                                        km_m, l_grid, pt_m, pt_reference,   &
-                                        rif_m, tend, zu, zw )
+                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
+                                     km_m, l_grid, vpt_m, pt_reference,  &
+                                     rif_m, tend, zu, zw )
+                ENDIF
+             ELSE
+                IF ( .NOT. humidity )  THEN
+                   IF ( ocean )  THEN
+                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
+                                        km, l_grid, rho, prho_reference,  &
+                                        rif, tend, zu, zw )
                    ELSE
-                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
-                                        km_m, l_grid, vpt_m, pt_reference,  &
-                                        rif_m, tend, zu, zw )
+                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
+                                        km, l_grid, pt, pt_reference, rif, &
+                                        tend, zu, zw )
                    ENDIF
                 ELSE
-                   IF ( .NOT. humidity )  THEN
-                      IF ( ocean )  THEN
-                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
-                                           km, l_grid, rho, prho_reference,  &
-                                           rif, tend, zu, zw )
-                      ELSE
-                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
-                                           km, l_grid, pt, pt_reference, rif, &
-                                           tend, zu, zw )
-                      ENDIF
-                   ELSE
-                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                        l_grid, vpt, pt_reference, rif, tend, &
-                                        zu, zw )
-                   ENDIF
-                ENDIF
-                CALL production_e( i, j )
-                CALL user_actions( i, j, 'e-tendency' )
-
-!
-!--             Prognostic equation for TKE.
-!--             Eliminate negative TKE values, which can occur due to numerical
-!--             reasons in the course of the integration. In such cases the old
-!--             TKE value is reduced by 90%.
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
-                                dt_3d * (                                      &
-                                   tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
-                                        )
-                   IF ( e_p(k,j,i) < 0.0 )  e_p(k,j,i) = 0.1 * e(k,j,i)
-                ENDDO
-
-!
-!--             Calculate tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      DO  k = nzb_s_inner(j,i)+1, nzt
-                         te_m(k,j,i) = tend(k,j,i)
-                      ENDDO
-                   ELSEIF ( intermediate_timestep_count < &
-                            intermediate_timestep_count_max )  THEN
-                      DO  k = nzb_s_inner(j,i)+1, nzt
-                         te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
-                                        5.3125 * te_m(k,j,i)
-                      ENDDO
-                   ENDIF
-                ENDIF
-
-             ENDIF   ! TKE equation
-
-          ENDIF   ! Gridpoints excluding the non-cyclic wall
+                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
+                                     l_grid, vpt, pt_reference, rif, tend, &
+                                     zu, zw )
+                ENDIF
+             ENDIF
+             CALL production_e( i, j )
+             CALL user_actions( i, j, 'e-tendency' )
+
+!
+!--          Prognostic equation for TKE.
+!--          Eliminate negative TKE values, which can occur due to numerical
+!--          reasons in the course of the integration. In such cases the old
+!--          TKE value is reduced by 90%.
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
+                             dt_3d * (                                      &
+                                tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
+                                     )
+                IF ( e_p(k,j,i) < 0.0 )  e_p(k,j,i) = 0.1 * e(k,j,i)
+             ENDDO
+
+!
+!--          Calculate tendencies for the next Runge-Kutta step
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( intermediate_timestep_count == 1 )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      te_m(k,j,i) = tend(k,j,i)
+                   ENDDO
+                ELSEIF ( intermediate_timestep_count < &
+                         intermediate_timestep_count_max )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                     5.3125 * te_m(k,j,i)
+                   ENDDO
+                ENDIF
+             ENDIF
+
+          ENDIF   ! TKE equation
 
        ENDDO
@@ -1268,5 +1266,5 @@
 !
 !-- Prognostic equation for u-velocity component
-    DO  i = nxl, nxr
+    DO  i = nxlu, nxr
        DO  j = nys, nyn
           DO  k = nzb_u_inner(j,i)+1, nzt
@@ -1285,5 +1283,5 @@
     IF ( timestep_scheme(1:5) == 'runge' )  THEN
        IF ( intermediate_timestep_count == 1 )  THEN
-          DO  i = nxl, nxr
+          DO  i = nxlu, nxr
              DO  j = nys, nyn
                 DO  k = nzb_u_inner(j,i)+1, nzt
@@ -1294,5 +1292,5 @@
        ELSEIF ( intermediate_timestep_count < &
                 intermediate_timestep_count_max )  THEN
-          DO  i = nxl, nxr
+          DO  i = nxlu, nxr
              DO  j = nys, nyn
                 DO  k = nzb_u_inner(j,i)+1, nzt
@@ -1340,5 +1338,5 @@
 !-- Prognostic equation for v-velocity component
     DO  i = nxl, nxr
-       DO  j = nys, nyn
+       DO  j = nysv, nyn
           DO  k = nzb_v_inner(j,i)+1, nzt
              v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) +    &
@@ -1357,5 +1355,5 @@
        IF ( intermediate_timestep_count == 1 )  THEN
           DO  i = nxl, nxr
-             DO  j = nys, nyn
+             DO  j = nysv, nyn
                 DO  k = nzb_v_inner(j,i)+1, nzt
                    tv_m(k,j,i) = tend(k,j,i)
@@ -1366,5 +1364,5 @@
                 intermediate_timestep_count_max )  THEN
           DO  i = nxl, nxr
-             DO  j = nys, nyn
+             DO  j = nysv, nyn
                 DO  k = nzb_v_inner(j,i)+1, nzt
                    tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 105)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 106)
@@ -4,5 +4,9 @@
 ! Actual revisions:
 ! -----------------
-! Call of new routine surface_coupler
+! Call of new routine surface_coupler,
+! presure solver is called after the first Runge-Kutta substep instead of the
+! last in case that call_psolver_at_all_substeps = .F.; for this case, the
+! random perturbation has to be added to the velocity fields also after the
+! first substep
 !
 ! Former revisions:
@@ -194,6 +198,9 @@
 !
 !--       Impose a random perturbation on the horizontal velocity field
-          IF ( create_disturbances  .AND.  &
+          IF ( create_disturbances  .AND.                                      &
+               ( call_psolver_at_all_substeps  .AND.                           &
                intermediate_timestep_count == intermediate_timestep_count_max )&
+          .OR. ( .NOT. call_psolver_at_all_substeps  .AND.                     &
+               intermediate_timestep_count == 1 ) )                            &
           THEN
              time_disturb = time_disturb + dt_3d
@@ -218,15 +225,8 @@
 !--       Reduce the velocity divergence via the equation for perturbation
 !--       pressure.
-          IF ( intermediate_timestep_count == intermediate_timestep_count_max &
-               .OR.  call_psolver_at_all_substeps )  THEN
+          IF ( intermediate_timestep_count == 1  .OR. &
+                call_psolver_at_all_substeps )  THEN
              CALL pres
           ENDIF
-
-!
-!--       In case of a non-cyclic lateral wall, set the boundary conditions for
-!--       the velocities at the outflow
-!          IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
-!             CALL boundary_conds( 'outflow_uvw' )
-!          ENDIF
 
 !
