Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 4182)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 4268)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Removing bulk cloud variables to respective module
+! 
+! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
 ! 
@@ -64,11 +67,8 @@
     USE arrays_3d,                                                             &
         ONLY:  c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l,  &
-               dzu, nc_p, nr_p, pt, pt_init, pt_p, q,                          &
-               q_p, qc_p, qr_p, s, s_p, sa, sa_p, u, u_init, u_m_l, u_m_n,     &
+               dzu, pt, pt_init, pt_p, q,                                      &
+               q_p, s, s_p, sa, sa_p, u, u_init, u_m_l, u_m_n,                 &
                u_m_r, u_m_s, u_p, v, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,  &
                w, w_p, w_m_l, w_m_n, w_m_r, w_m_s
-
-    USE bulk_cloud_model_mod,                                                  &
-        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
 
     USE chemistry_model_mod,                                                   &
@@ -102,5 +102,5 @@
  
     USE salsa_mod,                                                             &
-        ONLY:  salsa_boundary_conds        
+        ONLY:  salsa_boundary_conds
 
     USE surface_mod,                                                           &
@@ -136,5 +136,5 @@
        !$ACC PRESENT(bc_h, w_p)
        DO  m = 1, bc_h(l)%ns
-          i = bc_h(l)%i(m)            
+          i = bc_h(l)%i(m)
           j = bc_h(l)%j(m)
           k = bc_h(l)%k(m)
@@ -174,5 +174,5 @@
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, bc_h(l)%ns
-                i = bc_h(l)%i(m)            
+                i = bc_h(l)%i(m)
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
@@ -223,5 +223,5 @@
           !$OMP PARALLEL DO PRIVATE( i, j, k )
           DO  m = 1, bc_h(l)%ns
-             i = bc_h(l)%i(m)            
+             i = bc_h(l)%i(m)
              j = bc_h(l)%j(m)
              k = bc_h(l)%k(m)
@@ -253,5 +253,5 @@
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, bc_h(l)%ns
-                i = bc_h(l)%i(m)            
+                i = bc_h(l)%i(m)
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
@@ -259,11 +259,11 @@
              ENDDO
           ENDDO
-          
+
        ELSE
-          
+
           DO  l = 0, 1
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, bc_h(l)%ns
-                i = bc_h(l)%i(m)            
+                i = bc_h(l)%i(m)
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
@@ -279,50 +279,4 @@
           q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
        ENDIF
-
-       IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-!             
-!--       Surface conditions cloud water (Dirichlet)
-!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
-!--       the k coordinate belongs to the atmospheric grid point, therefore, set
-!--       qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
-          DO  l = 0, 1
-          !$OMP PARALLEL DO PRIVATE( i, j, k )
-             DO  m = 1, bc_h(l)%ns
-                i = bc_h(l)%i(m)            
-                j = bc_h(l)%j(m)
-                k = bc_h(l)%k(m)
-                qc_p(k+bc_h(l)%koff,j,i) = 0.0_wp
-                nc_p(k+bc_h(l)%koff,j,i) = 0.0_wp
-             ENDDO
-          ENDDO
-!
-!--       Top boundary condition for cloud water (Dirichlet)
-          qc_p(nzt+1,:,:) = 0.0_wp
-          nc_p(nzt+1,:,:) = 0.0_wp
-           
-       ENDIF
-
-       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-!             
-!--       Surface conditions rain water (Dirichlet)
-!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
-!--       the k coordinate belongs to the atmospheric grid point, therefore, set
-!--       qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
-          DO  l = 0, 1
-          !$OMP PARALLEL DO PRIVATE( i, j, k )
-             DO  m = 1, bc_h(l)%ns
-                i = bc_h(l)%i(m)            
-                j = bc_h(l)%j(m)
-                k = bc_h(l)%k(m)
-                qr_p(k+bc_h(l)%koff,j,i) = 0.0_wp
-                nr_p(k+bc_h(l)%koff,j,i) = 0.0_wp
-             ENDDO
-          ENDDO
-!
-!--       Top boundary condition for rain water (Dirichlet)
-          qr_p(nzt+1,:,:) = 0.0_wp
-          nr_p(nzt+1,:,:) = 0.0_wp
-           
-       ENDIF
     ENDIF
 !
@@ -336,9 +290,9 @@
 !--    s_p at k-1
        IF ( ibc_s_b == 0 ) THEN
-          
+
           DO  l = 0, 1
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, bc_h(l)%ns
-                i = bc_h(l)%i(m)            
+                i = bc_h(l)%i(m)
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
@@ -346,11 +300,11 @@
              ENDDO
           ENDDO
-          
+
        ELSE
-          
+
           DO  l = 0, 1
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, bc_h(l)%ns
-                i = bc_h(l)%i(m)            
+                i = bc_h(l)%i(m)
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
@@ -411,12 +365,4 @@
        IF ( humidity )  THEN
           q_p(:,nys-1,:) = q_p(:,nys,:)
-          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-             qc_p(:,nys-1,:) = qc_p(:,nys,:)
-             nc_p(:,nys-1,:) = nc_p(:,nys,:)
-          ENDIF
-          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-             qr_p(:,nys-1,:) = qr_p(:,nys,:)
-             nr_p(:,nys-1,:) = nr_p(:,nys,:)
-          ENDIF
        ENDIF
        IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
@@ -425,12 +371,4 @@
        IF ( humidity )  THEN
           q_p(:,nyn+1,:) = q_p(:,nyn,:)
-          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-             qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
-             nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
-          ENDIF
-          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-             qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
-             nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
-          ENDIF
        ENDIF
        IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
@@ -439,12 +377,4 @@
        IF ( humidity )  THEN
           q_p(:,:,nxl-1) = q_p(:,:,nxl)
-          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-             qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
-             nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
-          ENDIF
-          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-             qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
-             nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
-          ENDIF
        ENDIF
        IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
@@ -453,12 +383,4 @@
        IF ( humidity )  THEN
           q_p(:,:,nxr+1) = q_p(:,:,nxr)
-          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-             qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
-             nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
-          ENDIF
-          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
-             nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
-          ENDIF
        ENDIF
        IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
@@ -479,5 +401,5 @@
           u_p(:,-1,:) = u(:,0,:)
           v_p(:,0,:)  = v(:,1,:)
-          w_p(:,-1,:) = w(:,0,:)          
+          w_p(:,-1,:) = w(:,0,:)
        ELSEIF ( .NOT. use_cmax )  THEN
 
@@ -528,4 +450,143 @@
                 IF ( denom /= 0.0_wp )  THEN
                    c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
+                   IF ( c_w(k,i) < 0.0_wp )  THEN
+                      c_w(k,i) = 0.0_wp
+                   ELSEIF ( c_w(k,i) > c_max )  THEN
+                      c_w(k,i) = c_max
+                   ENDIF
+                ELSE
+                   c_w(k,i) = c_max
+                ENDIF
+
+                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
+                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
+                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
+
+             ENDDO
+          ENDDO
+
+#if defined( __parallel )
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
+          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dx, ierr )   
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
+          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dx, ierr )  
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
+          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dx, ierr )  
+#else
+          c_u_m = c_u_m_l
+          c_v_m = c_v_m_l
+          c_w_m = c_w_m_l
+#endif
+
+          c_u_m = c_u_m / (nx+1)
+          c_v_m = c_v_m / (nx+1)
+          c_w_m = c_w_m / (nx+1)
+
+!
+!--       Save old timelevels for the next timestep
+          IF ( intermediate_timestep_count == 1 )  THEN
+             u_m_s(:,:,:) = u(:,0:1,:)
+             v_m_s(:,:,:) = v(:,1:2,:)
+             w_m_s(:,:,:) = w(:,0:1,:)
+          ENDIF
+
+!
+!--       Calculate the new velocities
+          DO  k = nzb+1, nzt+1
+             DO  i = nxlg, nxrg
+                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *          &
+                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
+
+                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *          &
+                                       ( v(k,0,i) - v(k,1,i) ) * ddy
+
+                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *          &
+                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
+             ENDDO
+          ENDDO
+
+!
+!--       Bottom boundary at the outflow
+          IF ( ibc_uv_b == 0 )  THEN
+             u_p(nzb,-1,:) = 0.0_wp 
+             v_p(nzb,0,:)  = 0.0_wp  
+          ELSE
+             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
+             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
+          ENDIF
+          w_p(nzb,-1,:) = 0.0_wp
+
+!
+!--       Top boundary at the outflow
+          IF ( ibc_uv_t == 0 )  THEN
+             u_p(nzt+1,-1,:) = u_init(nzt+1)
+             v_p(nzt+1,0,:)  = v_init(nzt+1)
+          ELSE
+             u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
+             v_p(nzt+1,0,:)  = v_p(nzt,0,:)
+          ENDIF
+          w_p(nzt:nzt+1,-1,:) = 0.0_wp
+
+       ENDIF
+
+    ENDIF
+
+    IF ( bc_radiation_n )  THEN
+
+       IF ( use_cmax )  THEN
+          u_p(:,ny+1,:) = u(:,ny,:)
+          v_p(:,ny+1,:) = v(:,ny,:)
+          w_p(:,ny+1,:) = w(:,ny,:)
+       ELSEIF ( .NOT. use_cmax )  THEN
+
+          c_max = dy / dt_3d
+
+          c_u_m_l = 0.0_wp 
+          c_v_m_l = 0.0_wp
+          c_w_m_l = 0.0_wp
+
+          c_u_m = 0.0_wp 
+          c_v_m = 0.0_wp
+          c_w_m = 0.0_wp
+
+!
+!--       Calculate the phase speeds for u, v, and w, first local and then
+!--       average along the outflow boundary.
+          DO  k = nzb+1, nzt+1
+             DO  i = nxl, nxr
+
+                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
+
+                IF ( denom /= 0.0_wp )  THEN
+                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
+                   IF ( c_u(k,i) < 0.0_wp )  THEN
+                      c_u(k,i) = 0.0_wp
+                   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_wp )  THEN
+                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
+                   IF ( c_v(k,i) < 0.0_wp )  THEN
+                      c_v(k,i) = 0.0_wp
+                   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_wp )  THEN
+                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
                    IF ( c_w(k,i) < 0.0_wp )  THEN
                       c_w(k,i) = 0.0_wp
@@ -567,143 +628,4 @@
 !--       Save old timelevels for the next timestep
           IF ( intermediate_timestep_count == 1 )  THEN
-             u_m_s(:,:,:) = u(:,0:1,:)
-             v_m_s(:,:,:) = v(:,1:2,:)
-             w_m_s(:,:,:) = w(:,0:1,:)
-          ENDIF
-
-!
-!--       Calculate the new velocities
-          DO  k = nzb+1, nzt+1
-             DO  i = nxlg, nxrg
-                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *          &
-                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
-
-                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *          &
-                                       ( v(k,0,i) - v(k,1,i) ) * ddy
-
-                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *          &
-                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
-             ENDDO
-          ENDDO
-
-!
-!--       Bottom boundary at the outflow
-          IF ( ibc_uv_b == 0 )  THEN
-             u_p(nzb,-1,:) = 0.0_wp 
-             v_p(nzb,0,:)  = 0.0_wp  
-          ELSE                    
-             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
-             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
-          ENDIF
-          w_p(nzb,-1,:) = 0.0_wp
-
-!
-!--       Top boundary at the outflow
-          IF ( ibc_uv_t == 0 )  THEN
-             u_p(nzt+1,-1,:) = u_init(nzt+1)
-             v_p(nzt+1,0,:)  = v_init(nzt+1)
-          ELSE
-             u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
-             v_p(nzt+1,0,:)  = v_p(nzt,0,:)
-          ENDIF
-          w_p(nzt:nzt+1,-1,:) = 0.0_wp
-
-       ENDIF
-
-    ENDIF
-
-    IF ( bc_radiation_n )  THEN
-
-       IF ( use_cmax )  THEN
-          u_p(:,ny+1,:) = u(:,ny,:)
-          v_p(:,ny+1,:) = v(:,ny,:)
-          w_p(:,ny+1,:) = w(:,ny,:)          
-       ELSEIF ( .NOT. use_cmax )  THEN
-
-          c_max = dy / dt_3d
-
-          c_u_m_l = 0.0_wp 
-          c_v_m_l = 0.0_wp
-          c_w_m_l = 0.0_wp
-
-          c_u_m = 0.0_wp 
-          c_v_m = 0.0_wp
-          c_w_m = 0.0_wp
-
-!
-!--       Calculate the phase speeds for u, v, and w, first local and then
-!--       average along the outflow boundary.
-          DO  k = nzb+1, nzt+1
-             DO  i = nxl, nxr
-
-                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
-
-                IF ( denom /= 0.0_wp )  THEN
-                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
-                   IF ( c_u(k,i) < 0.0_wp )  THEN
-                      c_u(k,i) = 0.0_wp
-                   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_wp )  THEN
-                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
-                   IF ( c_v(k,i) < 0.0_wp )  THEN
-                      c_v(k,i) = 0.0_wp
-                   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_wp )  THEN
-                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
-                   IF ( c_w(k,i) < 0.0_wp )  THEN
-                      c_w(k,i) = 0.0_wp
-                   ELSEIF ( c_w(k,i) > c_max )  THEN
-                      c_w(k,i) = c_max
-                   ENDIF
-                ELSE
-                   c_w(k,i) = c_max
-                ENDIF
-
-                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
-                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
-                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
-
-             ENDDO
-          ENDDO
-
-#if defined( __parallel )   
-          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
-          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
-                              MPI_SUM, comm1dx, ierr )   
-          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
-          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
-                              MPI_SUM, comm1dx, ierr )  
-          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
-          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
-                              MPI_SUM, comm1dx, ierr )  
-#else
-          c_u_m = c_u_m_l
-          c_v_m = c_v_m_l
-          c_w_m = c_w_m_l
-#endif
-
-          c_u_m = c_u_m / (nx+1)
-          c_v_m = c_v_m / (nx+1)
-          c_w_m = c_w_m / (nx+1)
-
-!
-!--       Save old timelevels for the next timestep
-          IF ( intermediate_timestep_count == 1 )  THEN
                 u_m_n(:,:,:) = u(:,ny-1:ny,:)
                 v_m_n(:,:,:) = v(:,ny-1:ny,:)
@@ -731,5 +653,5 @@
              u_p(nzb,ny+1,:) = 0.0_wp
              v_p(nzb,ny+1,:) = 0.0_wp   
-          ELSE                    
+          ELSE
              u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
              v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
@@ -757,5 +679,5 @@
           u_p(:,:,0)  = u(:,:,1)
           v_p(:,:,-1) = v(:,:,0)
-          w_p(:,:,-1) = w(:,:,0)          
+          w_p(:,:,-1) = w(:,:,0)
        ELSEIF ( .NOT. use_cmax )  THEN
 
@@ -870,5 +792,5 @@
              u_p(nzb,:,0)  = 0.0_wp 
              v_p(nzb,:,-1) = 0.0_wp
-          ELSE                    
+          ELSE
              u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
              v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
@@ -896,5 +818,5 @@
           u_p(:,:,nx+1) = u(:,:,nx)
           v_p(:,:,nx+1) = v(:,:,nx)
-          w_p(:,:,nx+1) = w(:,:,nx)          
+          w_p(:,:,nx+1) = w(:,:,nx)
        ELSEIF ( .NOT. use_cmax )  THEN
 
@@ -1009,5 +931,5 @@
              u_p(nzb,:,nx+1) = 0.0_wp
              v_p(nzb,:,nx+1) = 0.0_wp 
-          ELSE                    
+          ELSE
              u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
              v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
