Index: /palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 72)
+++ /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 73)
@@ -33,4 +33,6 @@
 Changed:
 -------
+General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d
+
 +age_m in particle_type
 
@@ -44,9 +46,11 @@
 Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.
 
+q is not allowed to become negative (prognostic_equations).
+
 __vtk directives removed from main program.
 
 The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.
 
-check_parameters, data_output_dvrp, data_output_ptseries, data_output_ts, flow_statistics, header, init_particles, init_3d_model, modules, palm, package_parin, parin, time_integration
+check_parameters, data_output_dvrp, data_output_ptseries, data_output_ts, exchange_horiz_2d, flow_statistics, header, init_particles, init_pegrid, init_3d_model, modules, palm, package_parin, parin, prognostic_equations, read_3d_binary, time_integration, write_3d_binary
 
 
@@ -58,4 +62,6 @@
 Bugfix in sample for reading user defined data from restart file (user_init)
 
+Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)
+
 Check for possible negative humidities in the initial humidity profile.
 
@@ -64,3 +70,3 @@
 
 Makefile
-check_parameters, init_1d_model, user_interface
+check_parameters, init_1d_model, init_3d_model, user_interface
Index: /palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/boundary_conds.f90	(revision 72)
+++ /palm/trunk/SOURCE/boundary_conds.f90	(revision 73)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! The "main" part sets conditions for time level t+dt insteat of level t,
+! outflow boundary conditions changed from Neumann to radiation condition
 !
 ! Former revisions:
@@ -51,4 +52,6 @@
     INTEGER ::  i, j, k
 
+    REAL    ::  c_max, c_u, c_v, c_w, denom
+
 
     IF ( range == 'main')  THEN
@@ -56,13 +59,16 @@
 !--    Bottom boundary
        IF ( ibc_uv_b == 0 )  THEN
-          u(nzb,:,:) = -u(nzb+1,:,:)  ! satisfying the Dirichlet condition with
-          v(nzb,:,:) = -v(nzb+1,:,:)  ! an extra layer below the surface where
-       ELSE                           ! the u and v component change their sign
-          u(nzb,:,:) = u(nzb+1,:,:)
-          v(nzb,:,:) = v(nzb+1,:,:)
+!
+!--       Satisfying the Dirichlet condition with an extra layer below the
+!--       surface where the u and v component change their sign
+          u_p(nzb,:,:) = -u_p(nzb+1,:,:)
+          v_p(nzb,:,:) = -v_p(nzb+1,:,:)
+       ELSE
+          u_p(nzb,:,:) = u_p(nzb+1,:,:)
+          v_p(nzb,:,:) = v_p(nzb+1,:,:)
        ENDIF
        DO  i = nxl-1, nxr+1
           DO  j = nys-1, nyn+1
-             w(nzb_w_inner(j,i),j,i) = 0.0
+             w_p(nzb_w_inner(j,i),j,i) = 0.0
           ENDDO
        ENDDO
@@ -71,33 +77,24 @@
 !--    Top boundary
        IF ( ibc_uv_t == 0 )  THEN
-          u(nzt+1,:,:) = ug(nzt+1)
-          v(nzt+1,:,:) = vg(nzt+1)
+          u_p(nzt+1,:,:) = ug(nzt+1)
+          v_p(nzt+1,:,:) = vg(nzt+1)
        ELSE
-          u(nzt+1,:,:) = u(nzt,:,:)
-          v(nzt+1,:,:) = v(nzt,:,:)
-       ENDIF
-       w(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
+          u_p(nzt+1,:,:) = u_p(nzt,:,:)
+          v_p(nzt+1,:,:) = v_p(nzt,:,:)
+       ENDIF
+       w_p(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
 
 !
 !--    Temperature at bottom boundary
        IF ( ibc_pt_b == 0 )  THEN
-          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
-                   pt(nzb_s_inner(j,i),j,i) = pt_m(nzb_s_inner(j,i),j,i)
-                ENDDO
-             ENDDO
-          ELSE
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
-                   pt(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i),j,i)  
-                                            ! pt_m is not used for Runge-Kutta
-                ENDDO
-             ENDDO
-          ENDIF
+          DO  i = nxl-1, nxr+1
+             DO  j = nys-1, nyn+1
+                pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
+             ENDDO
+          ENDDO
        ELSE
           DO  i = nxl-1, nxr+1
              DO  j = nys-1, nyn+1
-                pt(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i)+1,j,i)
+                pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
              ENDDO
           ENDDO
@@ -107,13 +104,9 @@
 !--    Temperature at top boundary
        IF ( ibc_pt_t == 0 )  THEN
-          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-             pt(nzt+1,:,:) = pt_m(nzt+1,:,:)
-          ELSE
-             pt(nzt+1,:,:) = pt_p(nzt+1,:,:)  ! pt_m not used for Runge-Kutta
-          ENDIF
+          pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
        ELSEIF ( ibc_pt_t == 1 )  THEN
-          pt(nzt+1,:,:) = pt(nzt,:,:)
+          pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
        ELSEIF ( ibc_pt_t == 2 )  THEN
-          pt(nzt+1,:,:) = pt(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
+          pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
        ENDIF
 
@@ -124,8 +117,8 @@
           DO  i = nxl-1, nxr+1
              DO  j = nys-1, nyn+1
-                e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i)
-             ENDDO
-          ENDDO
-          e(nzt+1,:,:) = e(nzt,:,:)
+                e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
+             ENDDO
+          ENDDO
+          e_p(nzt+1,:,:) = e_p(nzt,:,:)
        ENDIF
 
@@ -137,21 +130,13 @@
 !--       Surface conditions for constant_moisture_flux
           IF ( ibc_q_b == 0 ) THEN
-             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-                DO  i = nxl-1, nxr+1
-                   DO  j = nys-1, nyn+1
-                      q(nzb_s_inner(j,i),j,i) = q_m(nzb_s_inner(j,i),j,i)
-                   ENDDO
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
                 ENDDO
-             ELSE
-                DO  i = nxl-1, nxr+1
-                   DO  j = nys-1, nyn+1
-                      q(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i),j,i)
-                   ENDDO                      ! q_m is not used for Runge-Kutta
-                ENDDO
-             ENDIF
+             ENDDO
           ELSE
              DO  i = nxl-1, nxr+1
                 DO  j = nys-1, nyn+1
-                   q(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i)+1,j,i)
+                   q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
                 ENDDO
              ENDDO
@@ -159,5 +144,5 @@
 !
 !--       Top boundary
-          q(nzt+1,:,:) = q(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
+          q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
        ENDIF
 
@@ -167,11 +152,11 @@
 !--    divergence. Dirichlet conditions are used for all other quantities.
        IF ( inflow_s )  THEN
-          v(:,nys,:) = v(:,nys-1,:)
+          v_p(:,nys,:) = v_p(:,nys-1,:)
        ELSEIF ( inflow_n ) THEN
-          v(:,nyn+vynp,:) = v(:,nyn+vynp+1,:)
+          v_p(:,nyn+vynp,:) = v_p(:,nyn+vynp+1,:)
        ELSEIF ( inflow_l ) THEN
-          u(:,:,nxl) = u(:,:,nxl-1)
+          u_p(:,:,nxl) = u_p(:,:,nxl-1)
        ELSEIF ( inflow_r ) THEN
-          u(:,:,nxr+uxrp) = u(:,:,nxr+uxrp+1)
+          u_p(:,:,nxr+uxrp) = u_p(:,:,nxr+uxrp+1)
        ENDIF
 
@@ -179,19 +164,19 @@
 !--    Lateral boundary conditions for scalar quantities at the outflow
        IF ( outflow_s )  THEN
-          pt(:,nys-1,:)     = pt(:,nys,:)
-          IF ( .NOT. constant_diffusion     )  e(:,nys-1,:) = e(:,nys,:)
-          IF ( moisture .OR. passive_scalar )  q(:,nys-1,:) = q(:,nys,:)
+          pt_p(:,nys-1,:)     = pt_p(:,nys,:)
+          IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
+          IF ( moisture .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
        ELSEIF ( outflow_n )  THEN
-          pt(:,nyn+1,:)     = pt(:,nyn,:)
-          IF ( .NOT. constant_diffusion     )  e(:,nyn+1,:) = e(:,nyn,:)
-          IF ( moisture .OR. passive_scalar )  q(:,nyn+1,:) = q(:,nyn,:)
+          pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
+          IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
+          IF ( moisture .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
        ELSEIF ( outflow_l )  THEN
-          pt(:,:,nxl-1)     = pt(:,:,nxl)
-          IF ( .NOT. constant_diffusion     )  e(:,:,nxl-1) = e(:,:,nxl)
-          IF ( moisture .OR. passive_scalar )  q(:,:,nxl-1) = q(:,:,nxl)      
+          pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
+          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
+          IF ( moisture .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)      
        ELSEIF ( outflow_r )  THEN
-          pt(:,:,nxr+1)     = pt(:,:,nxr)
-          IF ( .NOT. constant_diffusion     )  e(:,:,nxr+1) = e(:,:,nxr)
-          IF ( moisture .OR. passive_scalar )  q(:,:,nxr+1) = q(:,:,nxr)
+          pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
+          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
+          IF ( moisture .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
        ENDIF
 
@@ -200,81 +185,234 @@
     IF ( range == 'outflow_uvw' ) THEN
 !
-!--    Horizontal boundary conditions for the velocities at the outflow.
-!--    A Neumann condition is used for the wall normal velocity. The vertical
-!--    velocity is assumed as zero and a horizontal average along the wall is
-!--    used for the wall parallel horizontal velocity. The combination of all
-!--    three conditions ensures that the velocity field is free of divergence
-!--    at the outflow (uvmean_outflow_l is calculated in pres).
+!--    Radiation boundary condition for the velocities at the respective outflow
        IF ( outflow_s ) THEN
-          v(:,nys-1,:) = v(:,nys,:)
-          w(:,nys-1,:) = 0.0
-!
-!--       Compute the mean horizontal wind parallel to and within the outflow
-!--       wall and use this as boundary condition for u
-#if defined( __parallel )
-          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
-                              MPI_REAL, MPI_SUM, comm1dx, ierr )
-          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
-#else
-          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
-#endif
-          DO  k = nzb, nzt+1
-             u(k,nys-1,:) = uvmean_outflow(k)
-          ENDDO
-       ENDIF
-
-       IF ( outflow_n ) THEN
-          v(:,nyn+vynp+1,:) = v(:,nyn+vynp,:)
-          w(:,nyn+1,:)      = 0.0
-!
-!--       Compute the mean horizontal wind parallel to and within the outflow
-!--       wall and use this as boundary condition for u
-#if defined( __parallel )
-          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
-                              MPI_REAL, MPI_SUM, comm1dx, ierr )
-          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
-#else
-          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
-#endif
-          DO  k = nzb, nzt+1
-             u(k,nyn+1,:) = uvmean_outflow(k)
-          ENDDO
+!          v(:,nys-1,:) = v(:,nys,:)
+!          w(:,nys-1,:) = 0.0
+!!
+!!--       Compute the mean horizontal wind parallel to and within the outflow
+!!--       wall and use this as boundary condition for u
+!#if defined( __parallel )
+!          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
+!                              MPI_REAL, MPI_SUM, comm1dx, ierr )
+!          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
+!#else
+!          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
+!#endif
+!          DO  k = nzb, nzt+1
+!             u(k,nys-1,:) = uvmean_outflow(k)
+!          ENDDO
+       ENDIF
+
+       IF ( outflow_n  .AND.                                                 &
+            intermediate_timestep_count == intermediate_timestep_count_max ) &
+       THEN
+
+          c_max = dy / dt_3d
+
+          DO i = nxl-1, nxr+1
+             DO k = nzb+1, nzt+1
+
+!
+!--             First calculate the phase speeds for u,v, and w
+                denom = u_m_n(k,ny,i,-2) - u_m_n(k,ny-1,i,-2)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_u = -c_max * ( u_m_n(k,ny,i,-1)-u_m_n(k,ny,i,-2) ) / 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
+                ENDIF
+
+                denom = v_m_n(k,ny,i,-2) - v_m_n(k,ny-1,i,-2)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_v = -c_max * ( v_m_n(k,ny,i,-1)-v_m_n(k,ny,i,-2) ) / 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,-2) - w_m_n(k,ny-1,i,-2)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_w = -c_max * ( w_m_n(k,ny,i,-1)-w_m_n(k,ny,i,-2) ) / 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(k,ny+1,i) = u_m_n(k,ny+1,i,-1) - dt_3d * c_u * &
+                              ( u_m_n(k,ny+1,i,-1) - u_m_n(k,ny,i,-1) ) * ddy
+
+                v(k,ny+1,i) = v_m_n(k,ny+1,i,-1) - dt_3d * c_v * &
+                              ( v_m_n(k,ny+1,i,-1) - v_m_n(k,ny,i,-1) ) * ddy
+
+                w(k,ny+1,i) = w_m_n(k,ny+1,i,-1) - dt_3d * c_w * &
+                              ( w_m_n(k,ny+1,i,-1) - w_m_n(k,ny,i,-1) ) * ddy
+
+!
+!--             Swap timelevels for the next timestep
+                u_m_n(k,:,i,-2) = u_m_n(k,:,i,-1)
+                u_m_n(k,:,i,-1) = u(k,ny-1:ny+1,i)
+                v_m_n(k,:,i,-2) = v_m_n(k,:,i,-1)
+                v_m_n(k,:,i,-1) = v(k,ny-1:ny+1,i)
+                w_m_n(k,:,i,-2) = w_m_n(k,:,i,-1)
+                w_m_n(k,:,i,-1) = w(k,ny-1:ny+1,i)
+
+             ENDDO
+          ENDDO
+
+!
+!--       Bottom boundary at the outflow
+          IF ( ibc_uv_b == 0 )  THEN
+             u(nzb,ny+1,:) = -u(nzb+1,ny+1,:) 
+             v(nzb,ny+1,:) = -v(nzb+1,ny+1,:)  
+          ELSE                    
+             u(nzb,ny+1,:) =  u(nzb+1,ny+1,:)
+             v(nzb,ny+1,:) =  v(nzb+1,ny+1,:)
+          ENDIF
+          w(nzb,ny+1,:) = 0.0
+
+!
+!--       Top boundary at the outflow
+          IF ( ibc_uv_t == 0 )  THEN
+             u(nzt+1,ny+1,:) = ug(nzt+1)
+             v(nzt+1,ny+1,:) = vg(nzt+1)
+          ELSE
+             u(nzt+1,ny+1,:) = u(nzt,nyn+1,:)
+             v(nzt+1,ny+1,:) = v(nzt,nyn+1,:)
+          ENDIF
+          w(nzt:nzt+1,ny+1,:) = 0.0
+
        ENDIF
 
        IF ( outflow_l ) THEN
-          u(:,:,nxl-1) = u(:,:,nxl)
-          w(:,:,nxl-1) = 0.0
+!          u(:,:,nxl-1) = u(:,:,nxl)
+!          w(:,:,nxl-1) = 0.0
 !
 !--       Compute the mean horizontal wind parallel to and within the outflow
 !--       wall and use this as boundary condition for v
-#if defined( __parallel )
-          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
-                              MPI_REAL, MPI_SUM, comm1dy, ierr )
-          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
-#else
-          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
-#endif
-          DO  k = nzb, nzt+1
-             v(k,:,nxl-1) = uvmean_outflow(k)
-          ENDDO    
-
-       ENDIF
-
-       IF ( outflow_r ) THEN
-          u(:,:,nxr+uxrp+1) = u(:,:,nxr+uxrp)
-          w(:,:,nxr+1) = 0.0
-!
-!--       Compute the mean horizontal wind parallel to and within the outflow
-!--       wall and use this as boundary condition for v
-#if defined( __parallel )
-          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
-                              MPI_REAL, MPI_SUM, comm1dy, ierr )
-          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
-#else
-          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
-#endif
-          DO  k = nzb, nzt+1
-             v(k,:,nxr+1) = uvmean_outflow(k)
-          ENDDO    
+!#if defined( __parallel )
+!          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
+!                              MPI_REAL, MPI_SUM, comm1dy, ierr )
+!          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
+!#else
+!          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
+!#endif
+!          DO  k = nzb, nzt+1
+!             v(k,:,nxl-1) = uvmean_outflow(k)
+!          ENDDO    
+!
+       ENDIF
+
+       IF ( outflow_r  .AND.                                                 &
+            intermediate_timestep_count == intermediate_timestep_count_max ) &
+       THEN
+
+          c_max = dx / dt_3d
+
+          DO j = nys-1, nyn+1
+             DO k = nzb+1, nzt+1
+
+!
+!--             First calculate the phase speeds for u,v, and w
+                denom = u_m_r(k,j,nx,-2) - u_m_r(k,j,nx-1,-2)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_u = -c_max * ( u_m_r(k,j,nx,-1)-u_m_r(k,j,nx,-2) ) / 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
+                ENDIF
+
+                denom = v_m_r(k,j,nx,-2) - v_m_r(k,j,nx-1,-2)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_v = -c_max * ( v_m_r(k,j,nx,-1)-v_m_r(k,j,nx,-2) ) / 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,-2) - w_m_r(k,j,nx-1,-2)
+
+                IF ( denom /= 0.0 )  THEN
+                   c_w = -c_max * ( w_m_r(k,j,nx,-1)-w_m_n(k,j,nx,-2) ) / 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(k,j,nx+1) = u_m_r(k,j,nx+1,-1) - dt_3d * c_u * &
+                              ( u_m_r(k,j,nx+1,-1) - u_m_r(k,j,nx,-1) ) * ddx
+
+                v(k,j,nx+1) = v_m_r(k,j,nx+1,-1) - dt_3d * c_v * &
+                              ( v_m_r(k,j,nx+1,-1) - v_m_r(k,j,nx,-1) ) * ddx
+
+                w(k,j,nx+1) = w_m_r(k,j,nx+1,-1) - dt_3d * c_w * &
+                              ( w_m_r(k,j,nx+1,-1) - w_m_r(k,j,nx,-1) ) * ddx
+
+!
+!--             Swap timelevels for the next timestep
+                u_m_r(k,j,:,-2) = u_m_r(k,j,:,-1)
+                u_m_r(k,j,:,-1) = u(k,j,nx-1:nx+1)
+                v_m_r(k,j,:,-2) = v_m_r(k,j,:,-1)
+                v_m_r(k,j,:,-1) = v(k,j,nx-1:nx+1)
+                w_m_r(k,j,:,-2) = w_m_r(k,j,:,-1)
+                w_m_r(k,j,:,-1) = w(k,j,nx-1:nx+1)
+
+             ENDDO
+          ENDDO
+
+!
+!--       Bottom boundary at the outflow
+          IF ( ibc_uv_b == 0 )  THEN
+             u(nzb,ny+1,:) = -u(nzb+1,ny+1,:) 
+             v(nzb,ny+1,:) = -v(nzb+1,ny+1,:)  
+          ELSE                    
+             u(nzb,ny+1,:) =  u(nzb+1,ny+1,:)
+             v(nzb,ny+1,:) =  v(nzb+1,ny+1,:)
+          ENDIF
+          w(nzb,ny+1,:) = 0.0
+
+!
+!--       Top boundary at the outflow
+          IF ( ibc_uv_t == 0 )  THEN
+             u(nzt+1,ny+1,:) = ug(nzt+1)
+             v(nzt+1,ny+1,:) = vg(nzt+1)
+          ELSE
+             u(nzt+1,ny+1,:) = u(nzt,nyn+1,:)
+             v(nzt+1,ny+1,:) = v(nzt,nyn+1,:)
+          ENDIF
+          w(nzt:nzt+1,ny+1,:) = 0.0
+
        ENDIF
 
Index: /palm/trunk/SOURCE/calc_precipitation.f90
===================================================================
--- /palm/trunk/SOURCE/calc_precipitation.f90	(revision 72)
+++ /palm/trunk/SOURCE/calc_precipitation.f90	(revision 73)
@@ -71,12 +71,11 @@
                 tend(k,j,i) = tend(k,j,i) - dqdt_precip
 !
-!--             Precipitation rate in (kg * 0.001) / m**2 / s (because 1kg
-!--             gives 1 mm)
+!--             Precipitation rate in kg / m**2 / s (= mm/s)
                 precipitation_rate(j,i) = precipitation_rate(j,i) + &
-                                          dqdt_precip * dzw(k) * 0.001
+                                          dqdt_precip * dzw(k)
 
              ENDDO
 !
-!--          Sum up the precipitation amount (unit kg * 0.001 / m**2)
+!--          Sum up the precipitation amount, unit kg / m**2 (= mm)
              IF ( intermediate_timestep_count ==         &
                   intermediate_timestep_count_max  .AND. &
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 72)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 73)
@@ -8,5 +8,6 @@
 ! and pt_reference are checked,
 ! output of precipitation amount/rate and roughnes length + check
-! possible negative humidities are avoided in initial profile
+! possible negative humidities are avoided in initial profile,
+! dirichlet/neumann changed to dirichlet/radiation, etc.
 !
 ! Former revisions:
@@ -705,6 +706,6 @@
 !-- Check boundary conditions and set internal variables:
 !-- Lateral boundary conditions
-    IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/neumann'  .AND. &
-         bc_lr /= 'neumann/dirichlet' )  THEN
+    IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/radiation'  .AND. &
+         bc_lr /= 'radiation/dirichlet' )  THEN
        IF ( myid == 0 )  THEN
           PRINT*, '+++ check_parameters:'
@@ -713,6 +714,6 @@
        CALL local_stop
     ENDIF
-    IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/neumann'  .AND. &
-         bc_ns /= 'neumann/dirichlet' )  THEN
+    IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/radiation'  .AND. &
+         bc_ns /= 'radiation/dirichlet' )  THEN
        IF ( myid == 0 )  THEN
           PRINT*, '+++ check_parameters:'
@@ -2284,14 +2285,15 @@
     ENDIF
 
-    IF ( outflow_l )  THEN
+    IF ( bc_lr == 'radiation/dirichlet' )  THEN
        dist_nxr    = nx - inflow_disturbance_begin
        dist_nxl(1) = nx - inflow_disturbance_end
-    ELSEIF ( outflow_r )  THEN
+    ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
        dist_nxl    = inflow_disturbance_begin
        dist_nxr(1) = inflow_disturbance_end
-    ELSEIF ( outflow_s )  THEN
+    ENDIF
+    IF ( bc_ns == 'dirichlet/radiation' )  THEN
        dist_nyn    = ny - inflow_disturbance_begin
        dist_nys(1) = ny - inflow_disturbance_end
-    ELSEIF ( outflow_n )  THEN
+    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
        dist_nys    = inflow_disturbance_begin
        dist_nyn(1) = inflow_disturbance_end
Index: /palm/trunk/SOURCE/exchange_horiz_2d.f90
===================================================================
--- /palm/trunk/SOURCE/exchange_horiz_2d.f90	(revision 72)
+++ /palm/trunk/SOURCE/exchange_horiz_2d.f90	(revision 73)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary
+! conditions
 !
 ! Former revisions:
@@ -98,4 +99,12 @@
 #endif
 
+!
+!-- Neumann-conditions at inflow/outflow in case of non-cyclic boundary
+!-- conditions
+    IF ( inflow_l .OR. outflow_l )  ar(:,nxl-1) = ar(:,nxl)
+    IF ( inflow_r .OR. outflow_r )  ar(:,nxr+1) = ar(:,nxr)
+    IF ( inflow_s .OR. outflow_s )  ar(nys-1,:) = ar(nys,:)
+    IF ( inflow_n .OR. outflow_n )  ar(nyn+1,:) = ar(nyn,:)
+
     CALL cpu_log( log_point_s(13), 'exchange_horiz_2d', 'stop' )
 
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 72)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 73)
@@ -7,4 +7,7 @@
 ! Actual revisions:
 ! -----------------
+! Arrays for radiation boundary conditions are allocated (u_m_l, u_m_r, etc.),
+! bugfix for cases with the outflow damping layer extending over more than one
+! subdomain,
 ! New initializing action "by_user" calls user_init_3d_model,
 ! precipitation_amount/rate, ts_value are allocated, +module netcdf_control,
@@ -206,4 +209,28 @@
 
 !
+!-- Arrays to store velocity data from t-dt needed for radiation boundary
+!-- conditions
+    IF ( outflow_l )  THEN
+       ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1,-2:-1), &
+                 v_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1,-2:-1), &
+                 w_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1,-2:-1) )
+    ENDIF
+    IF ( outflow_r )  THEN
+       ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1,-2:-1), &
+                 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1,-2:-1), &
+                 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1,-2:-1) )
+    ENDIF
+    IF ( outflow_s )  THEN
+       ALLOCATE( u_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1,-2:-1), &
+                 v_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1,-2:-1), &
+                 w_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1,-2:-1) )
+    ENDIF
+    IF ( outflow_n )  THEN
+       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1,-2:-1), &
+                 v_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1,-2:-1), &
+                 w_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1,-2:-1) )
+    ENDIF
+
+!
 !-- Initial assignment of the pointers
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
@@ -726,4 +753,39 @@
        ENDIF
 
+!
+!--    Initialize old timelevels needed for radiation boundary conditions
+       IF ( outflow_l )  THEN
+          u_m_l(:,:,:,-2) = u(:,:,-1:1)
+          v_m_l(:,:,:,-2) = v(:,:,-1:1)
+          w_m_l(:,:,:,-2) = w(:,:,-1:1)
+          u_m_l(:,:,:,-1) = u(:,:,-1:1)
+          v_m_l(:,:,:,-1) = v(:,:,-1:1)
+          w_m_l(:,:,:,-1) = w(:,:,-1:1)
+       ENDIF
+       IF ( outflow_r )  THEN
+          u_m_r(:,:,:,-2) = u(:,:,nx-1:nx+1)
+          v_m_r(:,:,:,-2) = v(:,:,nx-1:nx+1)
+          w_m_r(:,:,:,-2) = w(:,:,nx-1:nx+1)
+          u_m_r(:,:,:,-1) = u(:,:,nx-1:nx+1)
+          v_m_r(:,:,:,-1) = v(:,:,nx-1:nx+1)
+          w_m_r(:,:,:,-1) = w(:,:,nx-1:nx+1)
+       ENDIF
+       IF ( outflow_s )  THEN
+          u_m_s(:,:,:,-2) = u(:,-1:1,:)
+          v_m_s(:,:,:,-2) = v(:,-1:1,:)
+          w_m_s(:,:,:,-2) = w(:,-1:1,:)
+          u_m_s(:,:,:,-1) = u(:,-1:1,:)
+          v_m_s(:,:,:,-1) = v(:,-1:1,:)
+          w_m_s(:,:,:,-1) = w(:,-1:1,:)
+       ENDIF
+       IF ( outflow_n )  THEN
+          u_m_n(:,:,:,-2) = u(:,ny-1:ny+1,:)
+          v_m_n(:,:,:,-2) = v(:,ny-1:ny+1,:)
+          w_m_n(:,:,:,-2) = w(:,ny-1:ny+1,:)
+          u_m_n(:,:,:,-1) = u(:,ny-1:ny+1,:)
+          v_m_n(:,:,:,-1) = v(:,ny-1:ny+1,:)
+          w_m_n(:,:,:,-1) = w(:,ny-1:ny+1,:)
+       ENDIF
+
     ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) &
     THEN
@@ -788,64 +850,58 @@
 !-- non-cyclic lateral boundaries. A linear increase is assumed over the first
 !-- half of the width of the damping layer
-    IF ( bc_lr /= 'cyclic' )  THEN
+    IF ( bc_lr == 'dirichlet/radiation' )  THEN
 
        DO  i = nxl-1, nxr+1
-
-          IF ( outflow_r )  THEN
-
-             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
-             ENDIF
-
-          ELSEIF ( outflow_l )  THEN
-
-             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
-             ENDIF
-
-          ENDIF
-
+          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
+          ENDIF
        ENDDO
-    ENDIF
-
-    IF ( bc_ns /= 'cyclic' )  THEN
+
+    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
+
+       DO  i = nxl-1, nxr+1
+          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
+          ENDIF
+       ENDDO
+
+    ENDIF
+
+    IF ( bc_ns == 'dirichlet/radiation' )  THEN
 
        DO  j = nys-1, nyn+1
-
-          IF ( outflow_n )  THEN
-
-             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
-
-          ELSEIF ( outflow_s )  THEN
-
-             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
-
-          ENDIF
-
+          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 == 'radiation/dirichlet' )  THEN
+
+       DO  j = nys-1, nyn+1
+          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_pegrid.f90
===================================================================
--- /palm/trunk/SOURCE/init_pegrid.f90	(revision 72)
+++ /palm/trunk/SOURCE/init_pegrid.f90	(revision 73)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! uxrp, vynp eliminated,
+! dirichlet/neumann changed to dirichlet/radiation, etc.
 !
 ! Former revisions:
@@ -495,6 +496,6 @@
 !
 !-- For non-cyclic boundaries extend array u (v) by one gridpoint
-    IF ( bc_lr /= 'cyclic' )  uxrp = 1
-    IF ( bc_ns /= 'cyclic' )  vynp = 1
+!    IF ( bc_lr /= 'cyclic' )  uxrp = 1
+!    IF ( bc_ns /= 'cyclic' )  vynp = 1
 
 !
@@ -743,7 +744,7 @@
 !-- by one gridpoint on the left/rightmost (northest/southest) processor
     IF ( pleft == MPI_PROC_NULL )  THEN
-       IF ( bc_lr == 'dirichlet/neumann' )  THEN
+       IF ( bc_lr == 'dirichlet/radiation' )  THEN
           inflow_l  = .TRUE.
-       ELSEIF ( bc_lr == 'neumann/dirichlet' )  THEN
+       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
           outflow_l = .TRUE.
        ENDIF
@@ -751,16 +752,16 @@
 
     IF ( pright == MPI_PROC_NULL )  THEN
-       IF ( bc_lr == 'dirichlet/neumann' )  THEN
+       IF ( bc_lr == 'dirichlet/radiation' )  THEN
           outflow_r = .TRUE.
-       ELSEIF ( bc_lr == 'neumann/dirichlet' )  THEN
+       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
           inflow_r  = .TRUE.
        ENDIF
-       uxrp      = 1
+!       uxrp      = 1
     ENDIF
 
     IF ( psouth == MPI_PROC_NULL )  THEN
-       IF ( bc_ns == 'dirichlet/neumann' )  THEN
+       IF ( bc_ns == 'dirichlet/radiation' )  THEN
           outflow_s = .TRUE.
-       ELSEIF ( bc_ns == 'neumann/dirichlet' )  THEN
+       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
           inflow_s  = .TRUE.
        ENDIF
@@ -768,10 +769,10 @@
 
     IF ( pnorth == MPI_PROC_NULL )  THEN
-       IF ( bc_ns == 'dirichlet/neumann' )  THEN
+       IF ( bc_ns == 'dirichlet/radiation' )  THEN
           inflow_n  = .TRUE.
-       ELSEIF ( bc_ns == 'neumann/dirichlet' )  THEN
+       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
           outflow_n = .TRUE.
        ENDIF
-       vynp      = 1
+!       vynp      = 1
     ENDIF
 
@@ -787,20 +788,20 @@
     ENDIF
 #else
-    IF ( bc_lr == 'dirichlet/neumann' )  THEN
+    IF ( bc_lr == 'dirichlet/radiation' )  THEN
        inflow_l  = .TRUE.
        outflow_r = .TRUE.
-       uxrp      = 1
-    ELSEIF ( bc_lr == 'neumann/dirichlet' )  THEN
+!       uxrp      = 1
+    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
        outflow_l = .TRUE.
        inflow_r  = .TRUE.
     ENDIF
 
-    IF ( bc_ns == 'dirichlet/neumann' )  THEN
+    IF ( bc_ns == 'dirichlet/radiation' )  THEN
        inflow_n  = .TRUE.
        outflow_s = .TRUE.
-    ELSEIF ( bc_ns == 'neumann/dirichlet' )  THEN
+    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
        outflow_n = .TRUE.
        inflow_s  = .TRUE.
-       vynp      = 1
+!       vynp      = 1
     ENDIF
 #endif
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 72)
+++ /palm/trunk/SOURCE/modules.f90	(revision 73)
@@ -6,5 +6,5 @@
 ! -----------------
 ! +arrays precipitation_amount, precipitation_rate, precipitation_rate_av,
-! rif_wall, z0_av
+! rif_wall, z0_av, +arrays u_m_l, u_m_r, etc. for radiation boundary conditions,
 ! +loop_optimization, netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*,
 ! id_var_zwwi_*, ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc, pt_reference,
@@ -111,5 +111,6 @@
           vpt, vpt_m, w, w_m, w_p
 
-    REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall
+    REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall, 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
                                                     
 
Index: /palm/trunk/SOURCE/pres.f90
===================================================================
--- /palm/trunk/SOURCE/pres.f90	(revision 72)
+++ /palm/trunk/SOURCE/pres.f90	(revision 73)
@@ -396,29 +396,29 @@
 !--       of non-cyclic boundary conditions). The respective mean velocity
 !--       is calculated from this in routine boundary_conds.
-          IF ( outflow_l  .AND.  i == nxl )  THEN
-             !$OMP CRITICAL
-             DO  k = nzb, nzt+1
-                uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxl)
-             ENDDO
-             !$OMP END CRITICAL
-          ELSEIF ( outflow_r  .AND.  i == nxr )  THEN
-             !$OMP CRITICAL
-             DO  k = nzb, nzt+1
-                uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxr)
-             ENDDO
-             !$OMP END CRITICAL
-          ELSEIF ( outflow_s  .AND.  j == nys )  THEN
-             !$OMP CRITICAL
-             DO  k = nzb, nzt+1
-                uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nys,i)
-             ENDDO
-             !$OMP END CRITICAL
-          ELSEIF ( outflow_n  .AND.  j == nyn )  THEN
-             !$OMP CRITICAL
-             DO  k = nzb, nzt+1
-                uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nyn,i)
-             ENDDO
-             !$OMP END CRITICAL
-          ENDIF
+!          IF ( outflow_l  .AND.  i == nxl )  THEN
+!             !$OMP CRITICAL
+!             DO  k = nzb, nzt+1
+!                uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxl)
+!             ENDDO
+!             !$OMP END CRITICAL
+!          ELSEIF ( outflow_r  .AND.  i == nxr )  THEN
+!             !$OMP CRITICAL
+!             DO  k = nzb, nzt+1
+!                uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxr)
+!             ENDDO
+!             !$OMP END CRITICAL
+!          ELSEIF ( outflow_s  .AND.  j == nys )  THEN
+!             !$OMP CRITICAL
+!             DO  k = nzb, nzt+1
+!                uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nys,i)
+!             ENDDO
+!             !$OMP END CRITICAL
+!          ELSEIF ( outflow_n  .AND.  j == nyn )  THEN
+!             !$OMP CRITICAL
+!             DO  k = nzb, nzt+1
+!                uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nyn,i)
+!             ENDDO
+!             !$OMP END CRITICAL
+!          ENDIF
 
 !
Index: /palm/trunk/SOURCE/production_e.f90
===================================================================
--- /palm/trunk/SOURCE/production_e.f90	(revision 72)
+++ /palm/trunk/SOURCE/production_e.f90	(revision 73)
@@ -910,6 +910,6 @@
 !--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
           !$OMP PARALLEL DO PRIVATE( ku, kv )
-          DO  i = nxl, nxr+1
-             DO  j = nys, nyn+1
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
 
                 ku = nzb_u_inner(j,i)+1
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 72)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 73)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! checking for negative q and limiting for positive values,
 ! z0 removed from arguments in calls of diffusion_u/v/w,
 ! subroutine names changed to .._noopt, .._cache, and .._vector
@@ -481,4 +482,5 @@
                                      ) -                                       &
                              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
 
@@ -915,4 +917,5 @@
                                         ) -                                    &
                                 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
 
@@ -1419,4 +1422,5 @@
                                      ) -                                       &
                              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
           ENDDO
Index: /palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/read_3d_binary.f90	(revision 72)
+++ /palm/trunk/SOURCE/read_3d_binary.f90	(revision 73)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! +precipitation_amount, precipitation_rate_av, rif_wall, z0_av
+! +precipitation_amount, precipitation_rate_av, rif_wall, u_m_l, u_m_r, etc.,
+! z0_av
 !
 ! Former revisions:
@@ -302,4 +303,12 @@
           CASE ( 'u_m' )
              READ ( 13 )  u_m
+          CASE ( 'u_m_l' )
+             READ ( 13 )  u_m_l
+          CASE ( 'u_m_n' )
+             READ ( 13 )  u_m_n
+          CASE ( 'u_m_r' )
+             READ ( 13 )  u_m_r
+          CASE ( 'u_m_s' )
+             READ ( 13 )  u_m_s
           CASE ( 'us' )
              READ ( 13 )  us
@@ -322,4 +331,12 @@
           CASE ( 'v_m' )
              READ (13 )   v_m
+          CASE ( 'v_m_l' )
+             READ ( 13 )  v_m_l
+          CASE ( 'v_m_n' )
+             READ ( 13 )  v_m_n
+          CASE ( 'v_m_r' )
+             READ ( 13 )  v_m_r
+          CASE ( 'v_m_s' )
+             READ ( 13 )  v_m_s
           CASE ( 'vpt' )
              READ ( 13 )  vpt
@@ -340,4 +357,12 @@
           CASE ( 'w_m' )
              READ ( 13 )  w_m
+          CASE ( 'w_m_l' )
+             READ ( 13 )  w_m_l
+          CASE ( 'w_m_n' )
+             READ ( 13 )  w_m_n
+          CASE ( 'w_m_r' )
+             READ ( 13 )  w_m_r
+          CASE ( 'w_m_s' )
+             READ ( 13 )  w_m_s
           CASE ( 'z0' )
              READ ( 13 )  z0
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 72)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 73)
@@ -158,11 +158,11 @@
 
 !
-!--       Swap the time levels in preparation for the next time step.
-          CALL swap_timelevel
-
-!
 !--       Boundary conditions for the prognostic quantities (except of the
 !--       velocities at the outflow in case of a non-cyclic lateral wall)
           CALL boundary_conds( 'main' )
+
+!
+!--       Swap the time levels in preparation for the next time step.
+          CALL swap_timelevel
 
 !
Index: /palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/write_3d_binary.f90	(revision 72)
+++ /palm/trunk/SOURCE/write_3d_binary.f90	(revision 73)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! +precipitation_amount, precipitation_rate_av, rif_wall, z0_av
+! +precipitation_amount, precipitation_rate_av, rif_wall, u_m_l, u_m_r, etc.,
+! z0_av
 !
 ! Former revisions:
@@ -169,4 +170,16 @@
     ENDIF
     WRITE ( 14 )  'u_m                 ';  WRITE ( 14 )  u_m
+    IF ( ALLOCATED( u_m_l ) )  THEN
+       WRITE ( 14 )  'u_m_l               ';  WRITE ( 14 )  u_m_l
+    ENDIF
+    IF ( ALLOCATED( u_m_n ) )  THEN
+       WRITE ( 14 )  'u_m_n               ';  WRITE ( 14 )  u_m_n
+    ENDIF
+    IF ( ALLOCATED( u_m_r ) )  THEN
+       WRITE ( 14 )  'u_m_r               ';  WRITE ( 14 )  u_m_r
+    ENDIF
+    IF ( ALLOCATED( u_m_s ) )  THEN
+       WRITE ( 14 )  'u_m_s               ';  WRITE ( 14 )  u_m_s
+    ENDIF
     WRITE ( 14 )  'us                  ';  WRITE ( 14 )  us
     WRITE ( 14 )  'usws                ';  WRITE ( 14 )  usws
@@ -184,4 +197,16 @@
     ENDIF
     WRITE ( 14 )  'v_m                 ';  WRITE ( 14 )  v_m
+    IF ( ALLOCATED( v_m_l ) )  THEN
+       WRITE ( 14 )  'v_m_l               ';  WRITE ( 14 )  v_m_l
+    ENDIF
+    IF ( ALLOCATED( v_m_n ) )  THEN
+       WRITE ( 14 )  'v_m_n               ';  WRITE ( 14 )  v_m_n
+    ENDIF
+    IF ( ALLOCATED( v_m_r ) )  THEN
+       WRITE ( 14 )  'v_m_r               ';  WRITE ( 14 )  v_m_r
+    ENDIF
+    IF ( ALLOCATED( v_m_s ) )  THEN
+       WRITE ( 14 )  'v_m_s               ';  WRITE ( 14 )  v_m_s
+    ENDIF
     IF ( moisture )  THEN
        WRITE ( 14 )  'vpt                 ';  WRITE ( 14 )  vpt
@@ -202,4 +227,16 @@
     ENDIF
     WRITE ( 14 )  'w_m                 ';  WRITE ( 14 )  w_m
+    IF ( ALLOCATED( w_m_l ) )  THEN
+       WRITE ( 14 )  'w_m_l               ';  WRITE ( 14 )  w_m_l
+    ENDIF
+    IF ( ALLOCATED( w_m_n ) )  THEN
+       WRITE ( 14 )  'w_m_n               ';  WRITE ( 14 )  w_m_n
+    ENDIF
+    IF ( ALLOCATED( w_m_r ) )  THEN
+       WRITE ( 14 )  'w_m_r               ';  WRITE ( 14 )  w_m_r
+    ENDIF
+    IF ( ALLOCATED( w_m_s ) )  THEN
+       WRITE ( 14 )  'w_m_s               ';  WRITE ( 14 )  w_m_s
+    ENDIF
     WRITE ( 14 )  'z0                  ';  WRITE ( 14 )  z0
     IF ( ALLOCATED( z0_av ) )  THEN
