Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 4280)
+++ palm/trunk/SOURCE/Makefile	(revision 4281)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# delete boundary_conds, added missing dependencies
+# 
+# 4270 2019-10-23 10:46:20Z monakurppa
 # - Implement offline nesting for salsa and add dependency of nesting_offl_mod
 #   on salsa_mod
@@ -127,5 +130,4 @@
 	basic_constants_and_equations_mod.f90 \
 	biometeorology_mod.f90 \
-	boundary_conds.f90 \
 	buoyancy.f90 \
 	calc_mean_profile.f90 \
@@ -345,12 +347,4 @@
 	palm_date_time_mod.o \
 	radiation_model_mod.o
-boundary_conds.o: \
-	bulk_cloud_model_mod.o \
-	chemistry_model_mod.o \
-	mod_kinds.o \
-	modules.o \
-	salsa_mod.o \
-	surface_mod.o \
-	turbulence_closure_mod.o
 bulk_cloud_model_mod.o: \
 	basic_constants_and_equations_mod.o \
@@ -553,4 +547,6 @@
 dynamics_mod.o: \
 	mod_kinds.o \
+        surface_mod.o \
+        pmc_interface_mod.o \
 	modules.o
 exchange_horiz.o: \
Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 4280)
+++ 	(revision )
@@ -1,917 +1,0 @@
-!> @file boundary_conds.f90
-!------------------------------------------------------------------------------!
-! This file is part of the PALM model system.
-!
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 1997-2019 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
-!
-! Current revisions:
-! -----------------
-! 
-! 
-! Former revisions:
-! -----------------
-! $Id$
-! Further modularization of boundary conditions: moved boundary conditions to
-! respective modules
-! 
-! 4268 2019-10-17 11:29:38Z schwenkel
-! Removing bulk cloud variables to respective module
-! 
-! 4182 2019-08-22 15:20:23Z scharf
-! Corrected "Former revisions" section
-! 
-! 4102 2019-07-17 16:00:03Z suehring
-! - Revise setting for boundary conditions at horizontal walls, use the offset
-!   index that belongs to the data structure instead of pre-calculating 
-!   the offset index for each facing.
-! - Set boundary conditions for bulk-cloud quantities also at downward facing
-!   surfaces
-! 
-! 4087 2019-07-11 11:35:23Z gronemeier
-! Add comment
-! 
-! 4086 2019-07-11 05:55:44Z gronemeier
-! Bugfix: use constant-flux layer condition for e in all rans modes
-! 
-! 3879 2019-04-08 20:25:23Z knoop
-! Bugfix, do not set boundary conditions for potential temperature in neutral
-! runs.
-! 
-! 3655 2019-01-07 16:51:22Z knoop
-! OpenACC port for SPEC
-!
-! Revision 1.1  1997/09/12 06:21:34  raasch
-! Initial revision
-!
-!
-! Description:
-! ------------
-!> Boundary conditions for the prognostic quantities.
-!> One additional bottom boundary condition is applied for the TKE (=(u*)**2)
-!> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
-!> handled in routine exchange_horiz. Pressure boundary conditions are
-!> explicitly set in routines pres, poisfft, poismg and sor.
-!------------------------------------------------------------------------------!
- SUBROUTINE boundary_conds
- 
-
-    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, pt, pt_init, pt_p, q,                                      &
-               q_p, s, s_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 control_parameters,                                                    &
-        ONLY:  bc_dirichlet_l,                                  &
-               bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
-               bc_radiation_s, bc_pt_t_val, bc_q_t_val, bc_s_t_val,            &
-               child_domain, coupling_mode, dt_3d,                             &
-               humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b,        &
-               ibc_s_t, ibc_uv_b, ibc_uv_t, intermediate_timestep_count,       &
-               nesting_offline, neutral, nudging, passive_scalar,  &
-               tsc, use_cmax
-
-    USE grid_variables,                                                        &
-        ONLY:  ddx, ddy, dx, dy
-
-    USE indices,                                                               &
-        ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
-
-    USE kinds
-
-    USE pegrid
-
-    USE pmc_interface,                                                         &
-        ONLY : nesting_mode
- 
-    USE surface_mod,                                                           &
-        ONLY :  bc_h
-
-    USE turbulence_closure_mod,                                                &
-        ONLY:  tcm_boundary_conds
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  i  !< grid index x direction
-    INTEGER(iwp) ::  j  !< grid index y direction
-    INTEGER(iwp) ::  k  !< grid index z direction
-    INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
-    INTEGER(iwp) ::  m  !< running index surface elements
-
-    REAL(wp)    ::  c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition
-    REAL(wp)    ::  denom !< horizontal gradient of velocity component normal to the outflow boundary 
-
-
-!
-!-- Bottom boundary 
-    IF ( ibc_uv_b == 1 )  THEN
-       u_p(nzb,:,:) = u_p(nzb+1,:,:)
-       v_p(nzb,:,:) = v_p(nzb+1,:,:)
-    ENDIF
-!
-!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
-!-- of downward-facing surfaces. 
-    DO  l = 0, 1
-       !$OMP PARALLEL DO PRIVATE( i, j, k )
-       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
-       !$ACC PRESENT(bc_h, w_p)
-       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)
-          w_p(k+bc_h(l)%koff,j,i) = 0.0_wp
-       ENDDO
-    ENDDO
-
-!
-!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
-    IF ( ibc_uv_t == 0 )  THEN
-        !$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init)
-        u_p(nzt+1,:,:) = u_init(nzt+1)
-        v_p(nzt+1,:,:) = v_init(nzt+1)
-        !$ACC END KERNELS
-    ELSEIF ( ibc_uv_t == 1 )  THEN
-        u_p(nzt+1,:,:) = u_p(nzt,:,:)
-        v_p(nzt+1,:,:) = v_p(nzt,:,:)
-    ENDIF
-
-!
-!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
-    IF (  .NOT.  child_domain  .AND.  .NOT.  nesting_offline  .AND.            &
-                 TRIM(coupling_mode) /= 'vnested_fine' )  THEN
-       !$ACC KERNELS PRESENT(w_p)
-       w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
-       !$ACC END KERNELS
-    ENDIF
-
-!
-!-- Temperature at bottom and top boundary.
-!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
-!-- the sea surface temperature of the coupled ocean model.
-!-- Dirichlet
-    IF ( .NOT. neutral )  THEN
-       IF ( ibc_pt_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)
-                j = bc_h(l)%j(m)
-                k = bc_h(l)%k(m)
-                pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i)
-             ENDDO
-          ENDDO
-!      
-!--    Neumann, zero-gradient
-       ELSEIF ( ibc_pt_b == 1 )  THEN
-          DO  l = 0, 1
-             !$OMP PARALLEL DO PRIVATE( i, j, k )
-             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
-             !$ACC PRESENT(bc_h, pt_p)
-             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)
-                pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i)
-             ENDDO
-          ENDDO
-       ENDIF
-       
-!      
-!--    Temperature at top boundary
-       IF ( ibc_pt_t == 0 )  THEN
-           pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
-!      
-!--        In case of nudging adjust top boundary to pt which is
-!--        read in from NUDGING-DATA
-           IF ( nudging )  THEN
-              pt_p(nzt+1,:,:) = pt_init(nzt+1)
-           ENDIF
-       ELSEIF ( ibc_pt_t == 1 )  THEN
-           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
-       ELSEIF ( ibc_pt_t == 2 )  THEN
-           !$ACC KERNELS PRESENT(pt_p, dzu)
-           pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
-           !$ACC END KERNELS
-       ENDIF
-    ENDIF
-!
-!-- Boundary conditions for total water content,
-!-- bottom and top boundary (see also temperature)
-    IF ( humidity )  THEN
-!
-!--    Surface conditions for constant_humidity_flux
-!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
-!--    the k coordinate belongs to the atmospheric grid point, therefore, set
-!--    q_p at k-1
-       IF ( ibc_q_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)
-                j = bc_h(l)%j(m)
-                k = bc_h(l)%k(m)
-                q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i)
-             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)
-                j = bc_h(l)%j(m)
-                k = bc_h(l)%k(m)
-                q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i)
-             ENDDO
-          ENDDO
-       ENDIF
-!
-!--    Top boundary
-       IF ( ibc_q_t == 0 ) THEN
-          q_p(nzt+1,:,:) = q(nzt+1,:,:)
-       ELSEIF ( ibc_q_t == 1 ) THEN
-          q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
-       ENDIF
-    ENDIF
-!
-!-- Boundary conditions for scalar,
-!-- bottom and top boundary (see also temperature)
-    IF ( passive_scalar )  THEN
-!
-!--    Surface conditions for constant_humidity_flux
-!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
-!--    the k coordinate belongs to the atmospheric grid point, therefore, set
-!--    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)
-                j = bc_h(l)%j(m)
-                k = bc_h(l)%k(m)
-                s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
-             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)
-                j = bc_h(l)%j(m)
-                k = bc_h(l)%k(m)
-                s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i)
-             ENDDO
-          ENDDO
-       ENDIF
-!
-!--    Top boundary condition
-       IF ( ibc_s_t == 0 )  THEN
-          s_p(nzt+1,:,:) = s(nzt+1,:,:)
-       ELSEIF ( ibc_s_t == 1 )  THEN
-          s_p(nzt+1,:,:) = s_p(nzt,:,:)
-       ELSEIF ( ibc_s_t == 2 )  THEN
-          s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
-       ENDIF
-
-    ENDIF  
-!
-!-- Set boundary conditions for subgrid TKE and dissipation (RANS mode)
-    CALL tcm_boundary_conds
-!
-!-- In case of inflow or nest boundary at the south boundary the boundary for v
-!-- is at nys and in case of inflow or nest boundary at the left boundary the
-!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
-!-- version) these levels are handled as a prognostic level, boundary values
-!-- have to be restored here.
-    IF ( bc_dirichlet_s )  THEN
-       v_p(:,nys,:) = v_p(:,nys-1,:)
-    ELSEIF ( bc_dirichlet_l ) THEN
-       u_p(:,:,nxl) = u_p(:,:,nxl-1)
-    ENDIF
-
-!
-!-- The same restoration for u at i=nxl and v at j=nys as above must be made
-!-- in case of nest boundaries. This must not be done in case of vertical nesting 
-!-- mode as in that case the lateral boundaries are actually cyclic. 
-!-- Lateral oundary conditions for TKE and dissipation are set 
-!-- in tcm_boundary_conds.
-    IF ( nesting_mode /= 'vertical'  .OR.  nesting_offline )  THEN
-       IF ( bc_dirichlet_s )  THEN
-          v_p(:,nys,:) = v_p(:,nys-1,:)
-       ENDIF
-       IF ( bc_dirichlet_l )  THEN
-          u_p(:,:,nxl) = u_p(:,:,nxl-1)
-       ENDIF
-    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
-       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
-       IF ( humidity )  THEN
-          q_p(:,nys-1,:) = q_p(:,nys,:)
-       ENDIF
-       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
-    ELSEIF ( bc_radiation_n )  THEN
-       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
-       IF ( humidity )  THEN
-          q_p(:,nyn+1,:) = q_p(:,nyn,:)
-       ENDIF
-       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
-    ELSEIF ( bc_radiation_l )  THEN
-       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
-       IF ( humidity )  THEN
-          q_p(:,:,nxl-1) = q_p(:,:,nxl)
-       ENDIF
-       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
-    ELSEIF ( bc_radiation_r )  THEN
-       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
-       IF ( humidity )  THEN
-          q_p(:,:,nxr+1) = q_p(:,:,nxr)
-       ENDIF
-       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
-    ENDIF
-
-!
-!-- Radiation boundary conditions for the velocities at the respective outflow.
-!-- The phase velocity is either assumed to the maximum phase velocity that
-!-- ensures numerical stability (CFL-condition) or calculated after
-!-- Orlanski(1976) and averaged along the outflow boundary.
-    IF ( bc_radiation_s )  THEN
-
-       IF ( use_cmax )  THEN
-          u_p(:,-1,:) = u(:,0,:)
-          v_p(:,0,:)  = v(:,1,:)
-          w_p(:,-1,:) = w(:,0,:)
-       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_s(k,0,i) - u_m_s(k,1,i)
-
-                IF ( denom /= 0.0_wp )  THEN
-                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,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_s(k,1,i) - v_m_s(k,2,i)
-
-                IF ( denom /= 0.0_wp )  THEN
-                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,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_s(k,0,i) - w_m_s(k,1,i)
-
-                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
-                   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,:)
-                w_m_n(:,:,:) = w(:,ny-1:ny,:)
-          ENDIF
-
-!
-!--       Calculate the new velocities
-          DO  k = nzb+1, nzt+1
-             DO  i = nxlg, nxrg
-                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *      &
-                                       ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
-
-                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *     &
-                                       ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
-
-                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *      &
-                                       ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
-             ENDDO
-          ENDDO
-
-!
-!--       Bottom boundary at the outflow
-          IF ( ibc_uv_b == 0 )  THEN
-             u_p(nzb,ny+1,:) = 0.0_wp
-             v_p(nzb,ny+1,:) = 0.0_wp   
-          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_wp
-
-!
-!--       Top boundary at the outflow
-          IF ( ibc_uv_t == 0 )  THEN
-             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
-             v_p(nzt+1,ny+1,:) = v_init(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_wp
-
-       ENDIF
-
-    ENDIF
-
-    IF ( bc_radiation_l )  THEN
-
-       IF ( use_cmax )  THEN
-          u_p(:,:,0)  = u(:,:,1)
-          v_p(:,:,-1) = v(:,:,0)
-          w_p(:,:,-1) = w(:,:,0)
-       ELSEIF ( .NOT. use_cmax )  THEN
-
-          c_max = dx / 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  j = nys, nyn
-
-                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
-
-                IF ( denom /= 0.0_wp )  THEN
-                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
-                   IF ( c_u(k,j) < 0.0_wp )  THEN
-                      c_u(k,j) = 0.0_wp
-                   ELSEIF ( c_u(k,j) > c_max )  THEN
-                      c_u(k,j) = c_max
-                   ENDIF
-                ELSE
-                   c_u(k,j) = c_max
-                ENDIF
-
-                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
-
-                IF ( denom /= 0.0_wp )  THEN
-                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
-                   IF ( c_v(k,j) < 0.0_wp )  THEN
-                      c_v(k,j) = 0.0_wp
-                   ELSEIF ( c_v(k,j) > c_max )  THEN
-                      c_v(k,j) = c_max
-                   ENDIF
-                ELSE
-                   c_v(k,j) = c_max
-                ENDIF
-
-                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
-
-                IF ( denom /= 0.0_wp )  THEN
-                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
-                   IF ( c_w(k,j) < 0.0_wp )  THEN
-                      c_w(k,j) = 0.0_wp
-                   ELSEIF ( c_w(k,j) > c_max )  THEN
-                      c_w(k,j) = c_max
-                   ENDIF
-                ELSE
-                   c_w(k,j) = c_max
-                ENDIF
-
-                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
-                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
-                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
-
-             ENDDO
-          ENDDO
-
-#if defined( __parallel )   
-          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
-          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
-                              MPI_SUM, comm1dy, ierr )   
-          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
-          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
-                              MPI_SUM, comm1dy, ierr )  
-          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
-          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
-                              MPI_SUM, comm1dy, 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 / (ny+1)
-          c_v_m = c_v_m / (ny+1)
-          c_w_m = c_w_m / (ny+1)
-
-!
-!--       Save old timelevels for the next timestep
-          IF ( intermediate_timestep_count == 1 )  THEN
-                u_m_l(:,:,:) = u(:,:,1:2)
-                v_m_l(:,:,:) = v(:,:,0:1)
-                w_m_l(:,:,:) = w(:,:,0:1)
-          ENDIF
-
-!
-!--       Calculate the new velocities
-          DO  k = nzb+1, nzt+1
-             DO  j = nysg, nyng
-                u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
-                                       ( u(k,j,0) - u(k,j,1) ) * ddx
-
-                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *          &
-                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
-
-                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *          &
-                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
-             ENDDO
-          ENDDO
-
-!
-!--       Bottom boundary at the outflow
-          IF ( ibc_uv_b == 0 )  THEN
-             u_p(nzb,:,0)  = 0.0_wp 
-             v_p(nzb,:,-1) = 0.0_wp
-          ELSE
-             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
-             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
-          ENDIF
-          w_p(nzb,:,-1) = 0.0_wp
-
-!
-!--       Top boundary at the outflow
-          IF ( ibc_uv_t == 0 )  THEN
-             u_p(nzt+1,:,0)  = u_init(nzt+1)
-             v_p(nzt+1,:,-1) = v_init(nzt+1)
-          ELSE
-             u_p(nzt+1,:,0)  = u_p(nzt,:,0)
-             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
-          ENDIF
-          w_p(nzt:nzt+1,:,-1) = 0.0_wp
-
-       ENDIF
-
-    ENDIF
-
-    IF ( bc_radiation_r )  THEN
-
-       IF ( use_cmax )  THEN
-          u_p(:,:,nx+1) = u(:,:,nx)
-          v_p(:,:,nx+1) = v(:,:,nx)
-          w_p(:,:,nx+1) = w(:,:,nx)
-       ELSEIF ( .NOT. use_cmax )  THEN
-
-          c_max = dx / 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  j = nys, nyn
-
-                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
-
-                IF ( denom /= 0.0_wp )  THEN
-                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
-                   IF ( c_u(k,j) < 0.0_wp )  THEN
-                      c_u(k,j) = 0.0_wp
-                   ELSEIF ( c_u(k,j) > c_max )  THEN
-                      c_u(k,j) = c_max
-                   ENDIF
-                ELSE
-                   c_u(k,j) = c_max
-                ENDIF
-
-                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
-
-                IF ( denom /= 0.0_wp )  THEN
-                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
-                   IF ( c_v(k,j) < 0.0_wp )  THEN
-                      c_v(k,j) = 0.0_wp
-                   ELSEIF ( c_v(k,j) > c_max )  THEN
-                      c_v(k,j) = c_max
-                   ENDIF
-                ELSE
-                   c_v(k,j) = c_max
-                ENDIF
-
-                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
-
-                IF ( denom /= 0.0_wp )  THEN
-                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
-                   IF ( c_w(k,j) < 0.0_wp )  THEN
-                      c_w(k,j) = 0.0_wp
-                   ELSEIF ( c_w(k,j) > c_max )  THEN
-                      c_w(k,j) = c_max
-                   ENDIF
-                ELSE
-                   c_w(k,j) = c_max
-                ENDIF
-
-                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
-                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
-                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
-
-             ENDDO
-          ENDDO
-
-#if defined( __parallel )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
-          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
-                              MPI_SUM, comm1dy, ierr )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
-          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
-                              MPI_SUM, comm1dy, ierr )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
-          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
-                              MPI_SUM, comm1dy, 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 / (ny+1)
-          c_v_m = c_v_m / (ny+1)
-          c_w_m = c_w_m / (ny+1)
-
-!
-!--       Save old timelevels for the next timestep
-          IF ( intermediate_timestep_count == 1 )  THEN
-                u_m_r(:,:,:) = u(:,:,nx-1:nx)
-                v_m_r(:,:,:) = v(:,:,nx-1:nx)
-                w_m_r(:,:,:) = w(:,:,nx-1:nx)
-          ENDIF
-
-!
-!--       Calculate the new velocities
-          DO  k = nzb+1, nzt+1
-             DO  j = nysg, nyng
-                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
-                                       ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
-
-                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *      &
-                                       ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
-
-                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *      &
-                                       ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
-             ENDDO
-          ENDDO
-
-!
-!--       Bottom boundary at the outflow
-          IF ( ibc_uv_b == 0 )  THEN
-             u_p(nzb,:,nx+1) = 0.0_wp
-             v_p(nzb,:,nx+1) = 0.0_wp 
-          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_wp
-
-!
-!--       Top boundary at the outflow
-          IF ( ibc_uv_t == 0 )  THEN
-             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
-             v_p(nzt+1,:,nx+1) = v_init(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_p(nzt:nzt+1,:,nx+1) = 0.0_wp
-
-       ENDIF
-
-    ENDIF
-
- END SUBROUTINE boundary_conds
Index: palm/trunk/SOURCE/dynamics_mod.f90
===================================================================
--- palm/trunk/SOURCE/dynamics_mod.f90	(revision 4280)
+++ palm/trunk/SOURCE/dynamics_mod.f90	(revision 4281)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Moved boundary conditions in dynamics module
+! 
+! 4097 2019-07-15 11:59:11Z suehring
 ! Avoid overlong lines - limit is 132 characters per line
 !
@@ -39,28 +42,76 @@
 
     USE arrays_3d, &
-        ONLY:  pt, pt_1, pt_2, pt_p, &
+        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, &
+               pt, pt_1, pt_2, pt_init, pt_p, &
                q, q_1, q_2, q_p, &
                s, s_1, s_2, s_p, &
-               u, u_1, u_2, u_p, &
-               v, v_1, v_2, v_p, &
-               w, w_1, w_2, w_p
+               u, u_1, u_2, u_init, u_p, u_m_l, u_m_n, u_m_r, u_m_s, &
+               v, v_1, v_2, v_p, v_init, v_m_l, v_m_n, v_m_r, v_m_s, &
+               w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s
 
     USE control_parameters, &
-        ONLY:  length, &
+        ONLY:  bc_dirichlet_l, &
+               bc_dirichlet_s, &
+               bc_radiation_l, &
+               bc_radiation_n, &
+               bc_radiation_r, &
+               bc_radiation_s, &
+               bc_pt_t_val, &
+               bc_q_t_val, &
+               bc_s_t_val, &
+               child_domain, &
+               coupling_mode, &
+               dt_3d, &
+               ibc_pt_b, &
+               ibc_pt_t, &
+               ibc_q_b, &
+               ibc_q_t, &
+               ibc_s_b, &
+               ibc_s_t, &
+               ibc_uv_b, &
+               ibc_uv_t, &
+               intermediate_timestep_count, &
+               length, &
+               nesting_offline, &
+               nudging, &
                restart_string, &
                humidity, &
                neutral, &
-               passive_scalar
+               passive_scalar, &
+               tsc, &
+               use_cmax
+
+    USE grid_variables, &
+        ONLY:  ddx, &
+               ddy, &
+               dx, &
+               dy
 
     USE indices, &
         ONLY:  nbgp, &
+               nx, &
                nxl, &
+               nxlg, &
                nxr, &
+               nxrg, &
+               ny, &
                nys, &
+               nysg, &
                nyn, &
+               nyng, &
                nzb, &
                nzt
 
     USE kinds
+
+    USE pegrid
+
+    USE pmc_interface, &
+        ONLY : nesting_mode
+
+    USE surface_mod, &
+        ONLY :  bc_h
+
 
     IMPLICIT NONE
@@ -90,4 +141,5 @@
        dynamics_exchange_horiz, &
        dynamics_prognostic_equations, &
+       dynamics_boundary_conditions, &
        dynamics_swap_timelevel, &
        dynamics_3d_data_averaging, &
@@ -169,4 +221,8 @@
     END INTERFACE dynamics_prognostic_equations
 
+    INTERFACE dynamics_boundary_conditions
+       MODULE PROCEDURE dynamics_boundary_conditions
+    END INTERFACE dynamics_boundary_conditions
+
     INTERFACE dynamics_swap_timelevel
        MODULE PROCEDURE dynamics_swap_timelevel
@@ -657,4 +713,818 @@
 
 
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute boundary conditions of dynamics model
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE dynamics_boundary_conditions
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i  !< grid index x direction
+    INTEGER(iwp) ::  j  !< grid index y direction
+    INTEGER(iwp) ::  k  !< grid index z direction
+    INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
+    INTEGER(iwp) ::  m  !< running index surface elements
+
+    REAL(wp)    ::  c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition
+    REAL(wp)    ::  denom !< horizontal gradient of velocity component normal to the outflow boundary
+
+!
+!-- Bottom boundary
+    IF ( ibc_uv_b == 1 )  THEN
+       u_p(nzb,:,:) = u_p(nzb+1,:,:)
+       v_p(nzb,:,:) = v_p(nzb+1,:,:)
+    ENDIF
+!
+!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
+!-- of downward-facing surfaces.
+    DO  l = 0, 1
+       !$OMP PARALLEL DO PRIVATE( i, j, k )
+       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+       !$ACC PRESENT(bc_h, w_p)
+       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)
+          w_p(k+bc_h(l)%koff,j,i) = 0.0_wp
+       ENDDO
+    ENDDO
+
+!
+!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
+    IF ( ibc_uv_t == 0 )  THEN
+        !$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init)
+        u_p(nzt+1,:,:) = u_init(nzt+1)
+        v_p(nzt+1,:,:) = v_init(nzt+1)
+        !$ACC END KERNELS
+    ELSEIF ( ibc_uv_t == 1 )  THEN
+        u_p(nzt+1,:,:) = u_p(nzt,:,:)
+        v_p(nzt+1,:,:) = v_p(nzt,:,:)
+    ENDIF
+
+!
+!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
+    IF (  .NOT.  child_domain  .AND.  .NOT.  nesting_offline  .AND.            &
+                 TRIM(coupling_mode) /= 'vnested_fine' )  THEN
+       !$ACC KERNELS PRESENT(w_p)
+       w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
+       !$ACC END KERNELS
+    ENDIF
+
+!
+!-- Temperature at bottom and top boundary.
+!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
+!-- the sea surface temperature of the coupled ocean model.
+!-- Dirichlet
+    IF ( .NOT. neutral )  THEN
+       IF ( ibc_pt_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)
+                j = bc_h(l)%j(m)
+                k = bc_h(l)%k(m)
+                pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i)
+             ENDDO
+          ENDDO
+!
+!--    Neumann, zero-gradient
+       ELSEIF ( ibc_pt_b == 1 )  THEN
+          DO  l = 0, 1
+             !$OMP PARALLEL DO PRIVATE( i, j, k )
+             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+             !$ACC PRESENT(bc_h, pt_p)
+             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)
+                pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i)
+             ENDDO
+          ENDDO
+       ENDIF
+
+!
+!--    Temperature at top boundary
+       IF ( ibc_pt_t == 0 )  THEN
+           pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
+!
+!--        In case of nudging adjust top boundary to pt which is
+!--        read in from NUDGING-DATA
+           IF ( nudging )  THEN
+              pt_p(nzt+1,:,:) = pt_init(nzt+1)
+           ENDIF
+       ELSEIF ( ibc_pt_t == 1 )  THEN
+           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
+       ELSEIF ( ibc_pt_t == 2 )  THEN
+           !$ACC KERNELS PRESENT(pt_p, dzu)
+           pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
+           !$ACC END KERNELS
+       ENDIF
+    ENDIF
+!
+!-- Boundary conditions for total water content,
+!-- bottom and top boundary (see also temperature)
+    IF ( humidity )  THEN
+!
+!--    Surface conditions for constant_humidity_flux
+!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
+!--    the k coordinate belongs to the atmospheric grid point, therefore, set
+!--    q_p at k-1
+       IF ( ibc_q_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)
+                j = bc_h(l)%j(m)
+                k = bc_h(l)%k(m)
+                q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i)
+             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)
+                j = bc_h(l)%j(m)
+                k = bc_h(l)%k(m)
+                q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i)
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    Top boundary
+       IF ( ibc_q_t == 0 ) THEN
+          q_p(nzt+1,:,:) = q(nzt+1,:,:)
+       ELSEIF ( ibc_q_t == 1 ) THEN
+          q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
+       ENDIF
+    ENDIF
+!
+!-- Boundary conditions for scalar,
+!-- bottom and top boundary (see also temperature)
+    IF ( passive_scalar )  THEN
+!
+!--    Surface conditions for constant_humidity_flux
+!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
+!--    the k coordinate belongs to the atmospheric grid point, therefore, set
+!--    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)
+                j = bc_h(l)%j(m)
+                k = bc_h(l)%k(m)
+                s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
+             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)
+                j = bc_h(l)%j(m)
+                k = bc_h(l)%k(m)
+                s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i)
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    Top boundary condition
+       IF ( ibc_s_t == 0 )  THEN
+          s_p(nzt+1,:,:) = s(nzt+1,:,:)
+       ELSEIF ( ibc_s_t == 1 )  THEN
+          s_p(nzt+1,:,:) = s_p(nzt,:,:)
+       ELSEIF ( ibc_s_t == 2 )  THEN
+          s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
+       ENDIF
+
+    ENDIF
+!
+!-- In case of inflow or nest boundary at the south boundary the boundary for v
+!-- is at nys and in case of inflow or nest boundary at the left boundary the
+!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
+!-- version) these levels are handled as a prognostic level, boundary values
+!-- have to be restored here.
+    IF ( bc_dirichlet_s )  THEN
+       v_p(:,nys,:) = v_p(:,nys-1,:)
+    ELSEIF ( bc_dirichlet_l ) THEN
+       u_p(:,:,nxl) = u_p(:,:,nxl-1)
+    ENDIF
+
+!
+!-- The same restoration for u at i=nxl and v at j=nys as above must be made
+!-- in case of nest boundaries. This must not be done in case of vertical nesting
+!-- mode as in that case the lateral boundaries are actually cyclic.
+!-- Lateral oundary conditions for TKE and dissipation are set
+!-- in tcm_boundary_conds.
+    IF ( nesting_mode /= 'vertical'  .OR.  nesting_offline )  THEN
+       IF ( bc_dirichlet_s )  THEN
+          v_p(:,nys,:) = v_p(:,nys-1,:)
+       ENDIF
+       IF ( bc_dirichlet_l )  THEN
+          u_p(:,:,nxl) = u_p(:,:,nxl-1)
+       ENDIF
+    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
+       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
+       IF ( humidity )  THEN
+          q_p(:,nys-1,:) = q_p(:,nys,:)
+       ENDIF
+       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
+    ELSEIF ( bc_radiation_n )  THEN
+       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
+       IF ( humidity )  THEN
+          q_p(:,nyn+1,:) = q_p(:,nyn,:)
+       ENDIF
+       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
+    ELSEIF ( bc_radiation_l )  THEN
+       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
+       IF ( humidity )  THEN
+          q_p(:,:,nxl-1) = q_p(:,:,nxl)
+       ENDIF
+       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
+    ELSEIF ( bc_radiation_r )  THEN
+       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
+       IF ( humidity )  THEN
+          q_p(:,:,nxr+1) = q_p(:,:,nxr)
+       ENDIF
+       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
+    ENDIF
+
+!
+!-- Radiation boundary conditions for the velocities at the respective outflow.
+!-- The phase velocity is either assumed to the maximum phase velocity that
+!-- ensures numerical stability (CFL-condition) or calculated after
+!-- Orlanski(1976) and averaged along the outflow boundary.
+    IF ( bc_radiation_s )  THEN
+
+       IF ( use_cmax )  THEN
+          u_p(:,-1,:) = u(:,0,:)
+          v_p(:,0,:)  = v(:,1,:)
+          w_p(:,-1,:) = w(:,0,:)
+       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_s(k,0,i) - u_m_s(k,1,i)
+
+                IF ( denom /= 0.0_wp )  THEN
+                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,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_s(k,1,i) - v_m_s(k,2,i)
+
+                IF ( denom /= 0.0_wp )  THEN
+                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,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_s(k,0,i) - w_m_s(k,1,i)
+
+                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
+                   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,:)
+                w_m_n(:,:,:) = w(:,ny-1:ny,:)
+          ENDIF
+
+!
+!--       Calculate the new velocities
+          DO  k = nzb+1, nzt+1
+             DO  i = nxlg, nxrg
+                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *      &
+                                       ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
+
+                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *     &
+                                       ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
+
+                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *      &
+                                       ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
+             ENDDO
+          ENDDO
+
+!
+!--       Bottom boundary at the outflow
+          IF ( ibc_uv_b == 0 )  THEN
+             u_p(nzb,ny+1,:) = 0.0_wp
+             v_p(nzb,ny+1,:) = 0.0_wp
+          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_wp
+
+!
+!--       Top boundary at the outflow
+          IF ( ibc_uv_t == 0 )  THEN
+             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
+             v_p(nzt+1,ny+1,:) = v_init(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_wp
+
+       ENDIF
+
+    ENDIF
+
+    IF ( bc_radiation_l )  THEN
+
+       IF ( use_cmax )  THEN
+          u_p(:,:,0)  = u(:,:,1)
+          v_p(:,:,-1) = v(:,:,0)
+          w_p(:,:,-1) = w(:,:,0)
+       ELSEIF ( .NOT. use_cmax )  THEN
+
+          c_max = dx / 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  j = nys, nyn
+
+                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
+
+                IF ( denom /= 0.0_wp )  THEN
+                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
+                   IF ( c_u(k,j) < 0.0_wp )  THEN
+                      c_u(k,j) = 0.0_wp
+                   ELSEIF ( c_u(k,j) > c_max )  THEN
+                      c_u(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_u(k,j) = c_max
+                ENDIF
+
+                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
+
+                IF ( denom /= 0.0_wp )  THEN
+                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
+                   IF ( c_v(k,j) < 0.0_wp )  THEN
+                      c_v(k,j) = 0.0_wp
+                   ELSEIF ( c_v(k,j) > c_max )  THEN
+                      c_v(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_v(k,j) = c_max
+                ENDIF
+
+                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
+
+                IF ( denom /= 0.0_wp )  THEN
+                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
+                   IF ( c_w(k,j) < 0.0_wp )  THEN
+                      c_w(k,j) = 0.0_wp
+                   ELSEIF ( c_w(k,j) > c_max )  THEN
+                      c_w(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_w(k,j) = c_max
+                ENDIF
+
+                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
+                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
+                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
+
+             ENDDO
+          ENDDO
+
+#if defined( __parallel )
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, ierr )
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, ierr )
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, 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 / (ny+1)
+          c_v_m = c_v_m / (ny+1)
+          c_w_m = c_w_m / (ny+1)
+
+!
+!--       Save old timelevels for the next timestep
+          IF ( intermediate_timestep_count == 1 )  THEN
+                u_m_l(:,:,:) = u(:,:,1:2)
+                v_m_l(:,:,:) = v(:,:,0:1)
+                w_m_l(:,:,:) = w(:,:,0:1)
+          ENDIF
+
+!
+!--       Calculate the new velocities
+          DO  k = nzb+1, nzt+1
+             DO  j = nysg, nyng
+                u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
+                                       ( u(k,j,0) - u(k,j,1) ) * ddx
+
+                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *          &
+                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
+
+                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *          &
+                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
+             ENDDO
+          ENDDO
+
+!
+!--       Bottom boundary at the outflow
+          IF ( ibc_uv_b == 0 )  THEN
+             u_p(nzb,:,0)  = 0.0_wp
+             v_p(nzb,:,-1) = 0.0_wp
+          ELSE
+             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
+             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
+          ENDIF
+          w_p(nzb,:,-1) = 0.0_wp
+
+!
+!--       Top boundary at the outflow
+          IF ( ibc_uv_t == 0 )  THEN
+             u_p(nzt+1,:,0)  = u_init(nzt+1)
+             v_p(nzt+1,:,-1) = v_init(nzt+1)
+          ELSE
+             u_p(nzt+1,:,0)  = u_p(nzt,:,0)
+             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
+          ENDIF
+          w_p(nzt:nzt+1,:,-1) = 0.0_wp
+
+       ENDIF
+
+    ENDIF
+
+    IF ( bc_radiation_r )  THEN
+
+       IF ( use_cmax )  THEN
+          u_p(:,:,nx+1) = u(:,:,nx)
+          v_p(:,:,nx+1) = v(:,:,nx)
+          w_p(:,:,nx+1) = w(:,:,nx)
+       ELSEIF ( .NOT. use_cmax )  THEN
+
+          c_max = dx / 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  j = nys, nyn
+
+                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
+
+                IF ( denom /= 0.0_wp )  THEN
+                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
+                   IF ( c_u(k,j) < 0.0_wp )  THEN
+                      c_u(k,j) = 0.0_wp
+                   ELSEIF ( c_u(k,j) > c_max )  THEN
+                      c_u(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_u(k,j) = c_max
+                ENDIF
+
+                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
+
+                IF ( denom /= 0.0_wp )  THEN
+                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
+                   IF ( c_v(k,j) < 0.0_wp )  THEN
+                      c_v(k,j) = 0.0_wp
+                   ELSEIF ( c_v(k,j) > c_max )  THEN
+                      c_v(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_v(k,j) = c_max
+                ENDIF
+
+                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
+
+                IF ( denom /= 0.0_wp )  THEN
+                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
+                   IF ( c_w(k,j) < 0.0_wp )  THEN
+                      c_w(k,j) = 0.0_wp
+                   ELSEIF ( c_w(k,j) > c_max )  THEN
+                      c_w(k,j) = c_max
+                   ENDIF
+                ELSE
+                   c_w(k,j) = c_max
+                ENDIF
+
+                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
+                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
+                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
+
+             ENDDO
+          ENDDO
+
+#if defined( __parallel )
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, ierr )
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, ierr )
+          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
+          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
+                              MPI_SUM, comm1dy, 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 / (ny+1)
+          c_v_m = c_v_m / (ny+1)
+          c_w_m = c_w_m / (ny+1)
+
+!
+!--       Save old timelevels for the next timestep
+          IF ( intermediate_timestep_count == 1 )  THEN
+                u_m_r(:,:,:) = u(:,:,nx-1:nx)
+                v_m_r(:,:,:) = v(:,:,nx-1:nx)
+                w_m_r(:,:,:) = w(:,:,nx-1:nx)
+          ENDIF
+
+!
+!--       Calculate the new velocities
+          DO  k = nzb+1, nzt+1
+             DO  j = nysg, nyng
+                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
+                                       ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
+
+                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *      &
+                                       ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
+
+                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *      &
+                                       ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
+             ENDDO
+          ENDDO
+
+!
+!--       Bottom boundary at the outflow
+          IF ( ibc_uv_b == 0 )  THEN
+             u_p(nzb,:,nx+1) = 0.0_wp
+             v_p(nzb,:,nx+1) = 0.0_wp
+          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_wp
+
+!
+!--       Top boundary at the outflow
+          IF ( ibc_uv_t == 0 )  THEN
+             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
+             v_p(nzt+1,:,nx+1) = v_init(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_p(nzt:nzt+1,:,nx+1) = 0.0_wp
+
+       ENDIF
+
+    ENDIF
+
+ END SUBROUTINE dynamics_boundary_conditions
 !------------------------------------------------------------------------------!
 ! Description:
Index: palm/trunk/SOURCE/module_interface.f90
===================================================================
--- palm/trunk/SOURCE/module_interface.f90	(revision 4280)
+++ palm/trunk/SOURCE/module_interface.f90	(revision 4281)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Added dynamics boundary conditions
+! 
+! 4272 2019-10-23 15:18:57Z schwenkel
 ! Further modularization of boundary conditions: moved boundary conditions to
 ! respective modules
@@ -179,4 +182,5 @@
                dynamics_exchange_horiz, &
                dynamics_prognostic_equations, &
+               dynamics_boundary_conditions, &
                dynamics_swap_timelevel, &
                dynamics_3d_data_averaging, &
@@ -197,4 +201,5 @@
                tcm_actions, &
                tcm_prognostic_equations, &
+               tcm_boundary_conds, &
                tcm_swap_timelevel, &
                tcm_3d_data_averaging, &
@@ -1299,4 +1304,7 @@
     IF ( debug_output_timestep )  CALL debug_message( 'module-specific boundary_conditions', 'start' )
 
+    CALL dynamics_boundary_conditions
+    CALL tcm_boundary_conds
+
     IF ( bulk_cloud_model    )  CALL bcm_boundary_conditions
     IF ( air_chemistry       )  CALL chem_boundary_conditions
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 4280)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 4281)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Moved boundary conditions to module interface
+! 
+! 4276 2019-10-28 16:03:29Z schwenkel
 ! Further modularization of lpm code components
 ! 
@@ -730,8 +733,6 @@
 !
 !--       Boundary conditions for the prognostic quantities (except of the
-!--       velocities at the outflow in case of a non-cyclic lateral wall)
-          CALL boundary_conds
-!
-!--       Boundary conditions for module-specific variables
+!--       velocities at the outflow in case of a non-cyclic lateral wall) and
+!--       boundary conditions for module-specific variables
           CALL module_interface_boundary_conditions
 !
