Index: /palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/boundary_conds.f90	(revision 4267)
+++ /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)
Index: /palm/trunk/SOURCE/bulk_cloud_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/bulk_cloud_model_mod.f90	(revision 4267)
+++ /palm/trunk/SOURCE/bulk_cloud_model_mod.f90	(revision 4268)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Introducing bcm_boundary_conditions
+! 
+! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
 ! 
@@ -114,5 +117,5 @@
                bc_radiation_n,                                                 &
                bc_radiation_r,                                                 &
-               bc_radiation_s,                                                 &                                                                
+               bc_radiation_s,                                                 &
                debug_output,                                                   &
                dt_3d, dt_do2d_xy, intermediate_timestep_count,                 &
@@ -247,4 +250,5 @@
            bcm_exchange_horiz, &
            bcm_prognostic_equations, &
+           bcm_boundary_conditions, &
            bcm_3d_data_averaging, &
            bcm_data_output_2d, &
@@ -313,5 +317,5 @@
     INTERFACE bcm_exchange_horiz
        MODULE PROCEDURE bcm_exchange_horiz
-    END INTERFACE bcm_exchange_horiz    
+    END INTERFACE bcm_exchange_horiz
 
     INTERFACE bcm_prognostic_equations
@@ -319,4 +323,8 @@
        MODULE PROCEDURE bcm_prognostic_equations_ij
     END INTERFACE bcm_prognostic_equations
+
+    INTERFACE bcm_boundary_conditions
+       MODULE PROCEDURE bcm_boundary_conditions
+    END INTERFACE bcm_boundary_conditions
 
     INTERFACE bcm_swap_timelevel
@@ -2014,4 +2022,107 @@
 
 !------------------------------------------------------------------------------!
+! Description: Boundary conditions of the bulk cloud module variables
+!------------------------------------------------------------------------------!
+    SUBROUTINE bcm_boundary_conditions
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i !<
+       INTEGER(iwp) ::  j !<
+       INTEGER(iwp) ::  k !<
+       INTEGER(iwp) ::  m !<
+       INTEGER(iwp) ::  l !<
+
+       IF ( 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 ( 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
+
+!
+!--    Lateral boundary conditions for scalar quantities at the outflow.
+!--    Lateral oundary conditions for TKE and dissipation are set
+!--    in tcm_boundary_conds.
+       IF ( bc_radiation_s )  THEN
+          IF ( microphysics_morrison )  THEN
+             qc_p(:,nys-1,:) = qc_p(:,nys,:)
+             nc_p(:,nys-1,:) = nc_p(:,nys,:)
+          ENDIF
+          IF ( microphysics_seifert )  THEN
+             qr_p(:,nys-1,:) = qr_p(:,nys,:)
+             nr_p(:,nys-1,:) = nr_p(:,nys,:)
+          ENDIF
+       ELSEIF ( bc_radiation_n )  THEN
+          IF ( microphysics_morrison )  THEN
+             qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
+             nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
+          ENDIF
+          IF ( microphysics_seifert )  THEN
+             qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
+             nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
+          ENDIF
+       ELSEIF ( bc_radiation_l )  THEN
+          IF ( microphysics_morrison )  THEN
+             qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
+             nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
+          ENDIF
+          IF ( microphysics_seifert )  THEN
+             qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
+             nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
+          ENDIF
+       ELSEIF ( bc_radiation_r )  THEN
+          IF ( microphysics_morrison )  THEN
+             qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
+             nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
+          ENDIF
+          IF ( microphysics_seifert )  THEN
+             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
+             nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
+          ENDIF
+       ENDIF
+
+    END SUBROUTINE bcm_boundary_conditions
+
+!------------------------------------------------------------------------------!
 !
 ! Description:
Index: /palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4267)
+++ /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4268)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Moving module specific boundary conditions from time_integration to module
+! 
+! 4230 2019-09-11 13:58:14Z suehring
 ! Bugfix, initialize mean profiles also in restart runs. Also initialize 
 ! array used for Runge-Kutta tendecies in restart runs.  
@@ -380,4 +383,8 @@
     END INTERFACE chem_boundary_conds
 
+    INTERFACE chem_boundary_conditions
+       MODULE PROCEDURE chem_boundary_conditions
+    END INTERFACE chem_boundary_conditions
+
     INTERFACE chem_check_data_output
        MODULE PROCEDURE chem_check_data_output
@@ -537,5 +544,5 @@
 
 
-    PUBLIC chem_3d_data_averaging, chem_boundary_conds,                       &
+    PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_boundary_conditions, &
             chem_boundary_conds_decycle, chem_check_data_output,              &
          chem_check_data_output_pr, chem_check_parameters,                    &
@@ -799,4 +806,33 @@
  END SUBROUTINE chem_boundary_conds
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine for boundary conditions
+!------------------------------------------------------------------------------!
+ SUBROUTINE chem_boundary_conditions
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  lsp             !<
+    INTEGER(iwp) ::  lsp_usr         !<
+
+!
+!--       Boundary conditions for prognostic quantitites of other modules:
+!--       Here, only decycling is carried out
+
+          DO  lsp = 1, nvar
+             lsp_usr = 1
+             DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
+                IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
+                   CALL chem_boundary_conds( chem_species(lsp)%conc_p,                          &
+                                             chem_species(lsp)%conc_pr_init )
+                ENDIF
+                lsp_usr = lsp_usr + 1
+             ENDDO
+          ENDDO
+
+
+ END SUBROUTINE chem_boundary_conditions
 
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/module_interface.f90
===================================================================
--- /palm/trunk/SOURCE/module_interface.f90	(revision 4267)
+++ /palm/trunk/SOURCE/module_interface.f90	(revision 4268)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Introduction of module_interface_boundary_conditions
+! 
+! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
 ! 
@@ -245,4 +248,5 @@
                bcm_exchange_horiz,                                             &
                bcm_prognostic_equations,                                       &
+               bcm_boundary_conditions,                                        &
                bcm_swap_timelevel,                                             &
                bcm_3d_data_averaging,                                          &
@@ -266,4 +270,5 @@
               chem_non_advective_processes,                                    &
               chem_prognostic_equations,                                       &
+              chem_boundary_conditions,                                        &
               chem_swap_timelevel,                                             &
               chem_3d_data_averaging,                                          &
@@ -418,4 +423,5 @@
                salsa_exchange_horiz_bounds,                                    &
                salsa_prognostic_equations,                                     &
+               salsa_boundary_conditions,                                      &
                salsa_swap_timelevel,                                           &
                salsa_3d_data_averaging,                                        &
@@ -520,4 +526,5 @@
        module_interface_exchange_horiz,                                        &
        module_interface_prognostic_equations,                                  &
+       module_interface_boundary_conditions,                                   &
        module_interface_swap_timelevel,                                        &
        module_interface_3d_data_averaging,                                     &
@@ -598,4 +605,8 @@
        MODULE PROCEDURE module_interface_swap_timelevel
     END INTERFACE module_interface_swap_timelevel
+
+    INTERFACE module_interface_boundary_conditions
+       MODULE PROCEDURE module_interface_boundary_conditions
+    END INTERFACE module_interface_boundary_conditions
 
     INTERFACE module_interface_3d_data_averaging
@@ -1273,4 +1284,22 @@
  END SUBROUTINE module_interface_prognostic_equations_ij
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute module-specific boundary conditions
+!------------------------------------------------------------------------------!
+ SUBROUTINE module_interface_boundary_conditions
+
+
+    IF ( debug_output_timestep )  CALL debug_message( 'module-specific boundary_conditions', 'start' )
+
+    IF ( bulk_cloud_model    )  CALL bcm_boundary_conditions
+    IF ( air_chemistry       )  CALL chem_boundary_conditions
+    IF ( salsa               )  CALL salsa_boundary_conditions
+
+    IF ( debug_output_timestep )  CALL debug_message( 'module-specific boundary_conditions', 'end' )
+
+
+ END SUBROUTINE module_interface_boundary_conditions
 
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- /palm/trunk/SOURCE/salsa_mod.f90	(revision 4267)
+++ /palm/trunk/SOURCE/salsa_mod.f90	(revision 4268)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! Moving module specific boundary conditions from time_integration to module
+! 
+! 4256 2019-10-07 10:08:52Z monakurppa
 ! Document previous changes: use global variables nx, ny and nz in salsa_header
 ! 
@@ -773,4 +776,8 @@
     END INTERFACE salsa_boundary_conds
 
+    INTERFACE salsa_boundary_conditions
+       MODULE PROCEDURE salsa_boundary_conditions
+    END INTERFACE salsa_boundary_conditions
+
     INTERFACE salsa_check_data_output
        MODULE PROCEDURE salsa_check_data_output
@@ -866,5 +873,5 @@
 !
 !-- Public functions:
-    PUBLIC salsa_boundary_conds, salsa_check_data_output, salsa_check_parameters,                  &
+    PUBLIC salsa_boundary_conds, salsa_check_data_output, salsa_check_parameters, salsa_boundary_conditions, &
            salsa_3d_data_averaging, salsa_data_output_2d, salsa_data_output_3d,                    &
            salsa_data_output_mask, salsa_define_netcdf_grid, salsa_diagnostics, salsa_driver,      &
@@ -7954,4 +7961,41 @@
  END SUBROUTINE salsa_tendency
 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Boundary conditions for prognostic variables in SALSA from module interface
+!------------------------------------------------------------------------------!
+ SUBROUTINE salsa_boundary_conditions
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  ib              !< index for aerosol size bins
+    INTEGER(iwp) ::  ic              !< index for aerosol mass bins
+    INTEGER(iwp) ::  icc             !< additional index for aerosol mass bins
+    INTEGER(iwp) ::  ig              !< index for salsa gases
+
+!
+!-- Boundary conditions for prognostic quantitites of other modules:
+!-- Here, only decycling is carried out
+    IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
+
+       DO  ib = 1, nbins_aerosol
+          CALL salsa_boundary_conds( aerosol_number(ib)%conc_p, aerosol_number(ib)%init )
+          DO  ic = 1, ncomponents_mass
+             icc = ( ic - 1 ) * nbins_aerosol + ib
+             CALL salsa_boundary_conds( aerosol_mass(icc)%conc_p, aerosol_mass(icc)%init )
+          ENDDO
+       ENDDO
+       IF ( .NOT. salsa_gases_from_chem )  THEN
+          DO  ig = 1, ngases_salsa
+             CALL salsa_boundary_conds( salsa_gas(ig)%conc_p, salsa_gas(ig)%init )
+          ENDDO
+       ENDIF
+
+    ENDIF
+
+ END SUBROUTINE salsa_boundary_conditions
+
 !------------------------------------------------------------------------------!
 ! Description:
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 4267)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 4268)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Removing module specific boundary conditions an put them into their modules
+! 
+! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
 ! 
@@ -254,5 +257,6 @@
 
     USE module_interface,                                                                          &
-        ONLY:  module_interface_actions, module_interface_swap_timelevel
+        ONLY:  module_interface_actions, module_interface_swap_timelevel,                          &
+               module_interface_boundary_conditions
 
     USE multi_agent_system_mod,                                                                    &
@@ -583,5 +587,5 @@
 !--    Execute all other module actions routunes
        CALL module_interface_actions( 'before_timestep' )
-       
+
 !
 !--    Start of intermediate step loop
@@ -750,40 +754,7 @@
 !--       velocities at the outflow in case of a non-cyclic lateral wall)
           CALL boundary_conds
-
-!
-!--       Boundary conditions for prognostic quantitites of other modules:
-!--       Here, only decycling is carried out
-          IF ( air_chemistry )  THEN
-
-             DO  lsp = 1, nvar
-                lsp_usr = 1
-                DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
-                   IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
-                      CALL chem_boundary_conds( chem_species(lsp)%conc_p,                          &
-                                                chem_species(lsp)%conc_pr_init )
-                   ENDIF
-                   lsp_usr = lsp_usr + 1
-                ENDDO
-             ENDDO
-
-          ENDIF
-
-          IF ( salsa  .AND.  time_since_reference_point >= skip_time_do_salsa )  THEN
-
-             DO  ib = 1, nbins_aerosol
-                CALL salsa_boundary_conds( aerosol_number(ib)%conc_p, aerosol_number(ib)%init )
-                DO  ic = 1, ncomponents_mass
-                   icc = ( ic - 1 ) * nbins_aerosol + ib
-                   CALL salsa_boundary_conds( aerosol_mass(icc)%conc_p, aerosol_mass(icc)%init )
-                ENDDO
-             ENDDO
-             IF ( .NOT. salsa_gases_from_chem )  THEN
-                DO  ig = 1, ngases_salsa
-                   CALL salsa_boundary_conds( salsa_gas(ig)%conc_p, salsa_gas(ig)%init )
-                ENDDO
-             ENDIF
-
-          ENDIF
-
+!
+!--       Boundary conditions for module-specific variables
+          CALL module_interface_boundary_conditions
 !
 !--       Incrementing timestep counter
@@ -1009,5 +980,5 @@
 
              ELSE
-!                
+!
 !--             Mass (volume) flux correction to ensure global mass conservation for child domains.
                 IF ( child_domain )  THEN
@@ -1018,5 +989,5 @@
                    ENDIF
                 ENDIF
-                
+
                 CALL pres
 
@@ -1084,5 +1055,5 @@
                 CALL lsm_energy_balance( .FALSE., 3 )
                 CALL lsm_soil_model( .FALSE., 3, .TRUE. )
-                
+
 !
 !--             At the end, set boundary conditons for potential temperature
@@ -1091,5 +1062,5 @@
                 CALL lsm_boundary_condition
 
-                
+
                 CALL cpu_log( log_point(54), 'land_surface', 'stop' )
              ENDIF
@@ -1099,5 +1070,5 @@
              IF (urban_surface) THEN
                 CALL cpu_log( log_point(74), 'urban_surface', 'start' )
-                
+
                 CALL usm_surface_energy_balance( .FALSE. )
                 IF ( usm_material_model )  THEN
