Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 73)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 75)
@@ -4,6 +4,7 @@
 ! 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
+! The "main" part sets conditions for time level t+dt instead of level t,
+! outflow boundary conditions changed from Neumann to radiation condition,
+! uxrp, vynp eliminated, moisture renamed humidity
 !
 ! Former revisions:
@@ -126,7 +127,7 @@
 !--    Boundary conditions for total water content or scalar,
 !--    bottom and surface boundary (see also temperature)
-       IF ( moisture  .OR.  passive_scalar )  THEN
-!
-!--       Surface conditions for constant_moisture_flux
+       IF ( humidity  .OR.  passive_scalar )  THEN
+!
+!--       Surface conditions for constant_humidity_flux
           IF ( ibc_q_b == 0 ) THEN
              DO  i = nxl-1, nxr+1
@@ -154,9 +155,9 @@
           v_p(:,nys,:) = v_p(:,nys-1,:)
        ELSEIF ( inflow_n ) THEN
-          v_p(:,nyn+vynp,:) = v_p(:,nyn+vynp+1,:)
+          v_p(:,nyn,:) = v_p(:,nyn+1,:)
        ELSEIF ( inflow_l ) THEN
           u_p(:,:,nxl) = u_p(:,:,nxl-1)
        ELSEIF ( inflow_r ) THEN
-          u_p(:,:,nxr+uxrp) = u_p(:,:,nxr+uxrp+1)
+          u_p(:,:,nxr) = u_p(:,:,nxr+1)
        ENDIF
 
@@ -166,254 +167,395 @@
           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,:)
+          IF ( humidity .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
        ELSEIF ( outflow_n )  THEN
           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,:)
+          IF ( humidity .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
        ELSEIF ( outflow_l )  THEN
           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)      
+          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)
        ELSEIF ( outflow_r )  THEN
           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)
+          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
        ENDIF
 
     ENDIF
 
-    IF ( range == 'outflow_uvw' ) THEN
-!
-!--    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  .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
+!
+!-- Radiation boundary condition for the velocities at the respective outflow
+    IF ( outflow_s  .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_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
-
-                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
+             ELSE
+                c_u = c_max
+             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
-
-                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
+             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
-
-!
-!--             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
-!
-!--       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  .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
+             ELSE
+                c_w = c_max
+             ENDIF
+
+!
+!--          Calculate the new velocities
+             u_p(k,-1,i) = u(k,-1,i) + dt_3d * c_u * &
+                                       ( 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 * &
+                                       ( 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
+       ENDDO
+
+!
+!--    Bottom boundary at the outflow
+       IF ( ibc_uv_b == 0 )  THEN
+          u_p(nzb,-1,:) = -u_p(nzb+1,-1,:) 
+          v_p(nzb,-1,:) = -v_p(nzb+1,-1,:)  
+       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
+
+!
+!--    Top boundary at the outflow
+       IF ( ibc_uv_t == 0 )  THEN
+          u_p(nzt+1,-1,:) = ug(nzt+1)
+          v_p(nzt+1,-1,:) = vg(nzt+1)
+       ELSE
+          u_p(nzt+1,-1,:) = u(nzt,-1,:)
+          v_p(nzt+1,-1,:) = v(nzt,-1,:)
+       ENDIF
+       w_p(nzt:nzt+1,-1,:) = 0.0
+
+    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) - 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
-
-                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
+             ELSE
+                c_u = c_max
+             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
-
-                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
+             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
-
-!
-!--             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
+             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(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(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(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
+       ENDDO
+
+!
+!--    Bottom boundary at the outflow
+       IF ( ibc_uv_b == 0 )  THEN
+          u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:) 
+          v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:)  
+       ELSE                    
+          u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
+          v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
+       ENDIF
+       w_p(nzb,ny+1,:) = 0.0
+
+!
+!--    Top boundary at the outflow
+       IF ( ibc_uv_t == 0 )  THEN
+          u_p(nzt+1,ny+1,:) = ug(nzt+1)
+          v_p(nzt+1,ny+1,:) = vg(nzt+1)
+       ELSE
+          u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
+          v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
+       ENDIF
+       w_p(nzt:nzt+1,ny+1,:) = 0.0
+
+    ENDIF
+
+    IF ( outflow_l  .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_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
+             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 * &
+                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
+
+             w_p(k,j,-1) = w(k,j,-1) + dt_3d * c_w * &
+                                       ( 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
+       ENDDO
+
+!
+!--    Bottom boundary at the outflow
+       IF ( ibc_uv_b == 0 )  THEN
+          u_p(nzb,:,-1) = -u_p(nzb+1,:,-1) 
+          v_p(nzb,:,-1) = -v_p(nzb+1,:,-1)  
+       ELSE                    
+          u_p(nzb,:,-1) =  u_p(nzb+1,:,-1)
+          v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
+       ENDIF
+       w_p(nzb,:,-1) = 0.0
+
+!
+!--    Top boundary at the outflow
+       IF ( ibc_uv_t == 0 )  THEN
+          u_p(nzt+1,:,-1) = ug(nzt+1)
+          v_p(nzt+1,:,-1) = vg(nzt+1)
+       ELSE
+          u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
+          v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
+       ENDIF
+       w_p(nzt:nzt+1,:,-1) = 0.0
+
+    ENDIF
+
+    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) - 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
+             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(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(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(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
+       ENDDO
+
+!
+!--    Bottom boundary at the outflow
+       IF ( ibc_uv_b == 0 )  THEN
+          u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1) 
+          v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1)  
+       ELSE                    
+          u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
+          v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
+       ENDIF
+       w_p(nzb,:,nx+1) = 0.0
+
+!
+!--    Top boundary at the outflow
+       IF ( ibc_uv_t == 0 )  THEN
+          u_p(nzt+1,:,nx+1) = ug(nzt+1)
+          v_p(nzt+1,:,nx+1) = vg(nzt+1)
+       ELSE
+          u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
+          v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
+       ENDIF
+       w(nzt:nzt+1,:,nx+1) = 0.0
 
     ENDIF
