Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 4101)
+++ palm/trunk/SOURCE/Makefile	(revision 4102)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Remove dependency on pmc_interface for boundary_conds
+# 
+# 4070 2019-07-03 13:51:40Z gronemeier
 # Add new data output modules
 #
@@ -802,10 +805,8 @@
 	radiation_model_mod.o
 boundary_conds.o: \
-	basic_constants_and_equations_mod.o \
 	bulk_cloud_model_mod.o \
 	chemistry_model_mod.o \
 	mod_kinds.o \
 	modules.o \
-	pmc_interface_mod.o \
 	salsa_mod.o \
 	surface_mod.o \
Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 4101)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 4102)
@@ -25,4 +25,11 @@
 ! -----------------
 ! $Id$
+! - 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
 ! 
@@ -224,12 +231,9 @@
     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,  &
-               diss, diss_p, dzu, e_p, nc_p, nr_p, pt, pt_init, pt_p, q,       &
+               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,     &
                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 basic_constants_and_equations_mod,                                     &
-        ONLY:  kappa
-
     USE bulk_cloud_model_mod,                                                  &
         ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
@@ -239,12 +243,12 @@
 
     USE control_parameters,                                                    &
-        ONLY:  air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,  &
+        ONLY:  air_chemistry, 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, constant_diffusion, coupling_mode, dt_3d,         &
+               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, ocean_mode, passive_scalar,  &
-               rans_mode, rans_tke_e, tsc, salsa, use_cmax
+               tsc, salsa, use_cmax
 
     USE grid_variables,                                                        &
@@ -262,5 +266,5 @@
 
     USE pmc_interface,                                                         &
-        ONLY : nesting_mode, rans_mode_parent
+        ONLY : nesting_mode
  
     USE salsa_mod,                                                             &
@@ -268,9 +272,8 @@
 
     USE surface_mod,                                                           &
-        ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
-                surf_usm_h, surf_usm_v
+        ONLY :  bc_h
 
     USE turbulence_closure_mod,                                                &
-        ONLY:  c_0
+        ONLY:  tcm_boundary_conds
 
     IMPLICIT NONE
@@ -279,5 +282,4 @@
     INTEGER(iwp) ::  j  !< grid index y direction
     INTEGER(iwp) ::  k  !< grid index z direction
-    INTEGER(iwp) ::  kb !< variable to set respective boundary value, depends on facing. 
     INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
     INTEGER(iwp) ::  m  !< running index surface elements
@@ -297,8 +299,4 @@
 !-- of downward-facing surfaces. 
     DO  l = 0, 1
-!
-!--    Set kb, for upward-facing surfaces value at topography top (k-1) is set,
-!--    for downward-facing surfaces at topography bottom (k+1). 
-       kb = MERGE( -1, 1, l == 0 )
        !$OMP PARALLEL DO PRIVATE( i, j, k )
        !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
@@ -308,5 +306,5 @@
           j = bc_h(l)%j(m)
           k = bc_h(l)%k(m)
-          w_p(k+kb,j,i) = 0.0_wp
+          w_p(k+bc_h(l)%koff,j,i) = 0.0_wp
        ENDDO
     ENDDO
@@ -341,8 +339,4 @@
        IF ( ibc_pt_b == 0 )  THEN
           DO  l = 0, 1
-!      
-!--          Set kb, for upward-facing surfaces value at topography top (k-1)
-!--          is set, for downward-facing surfaces at topography bottom (k+1). 
-             kb = MERGE( -1, 1, l == 0 )
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, bc_h(l)%ns
@@ -350,5 +344,5 @@
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
-                pt_p(k+kb,j,i) = pt(k+kb,j,i)
+                pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i)
              ENDDO
           ENDDO
@@ -357,8 +351,4 @@
        ELSEIF ( ibc_pt_b == 1 )  THEN
           DO  l = 0, 1
-!      
-!--          Set kb, for upward-facing surfaces value at topography top (k-1) 
-!--          is set, for downward-facing surfaces at topography bottom (k+1). 
-             kb = MERGE( -1, 1, l == 0 )
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
@@ -368,5 +358,5 @@
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
-                pt_p(k+kb,j,i) = pt_p(k,j,i)
+                pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i)
              ENDDO
           ENDDO
@@ -391,192 +381,4 @@
        ENDIF
     ENDIF
-
-!
-!-- Boundary conditions for TKE.
-!-- Generally Neumann conditions with de/dz=0 are assumed.
-    IF ( .NOT. constant_diffusion )  THEN
-
-       IF ( .NOT. rans_mode )  THEN
-          DO  l = 0, 1
-!
-!--         Set kb, for upward-facing surfaces value at topography top (k-1) is set,
-!--         for downward-facing surfaces at topography bottom (k+1). 
-             kb = MERGE( -1, 1, l == 0 )
-             !$OMP PARALLEL DO PRIVATE( i, j, k )
-             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
-             !$ACC PRESENT(bc_h, e_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)
-                e_p(k+kb,j,i) = e_p(k,j,i)
-             ENDDO
-          ENDDO
-       ELSE
-!
-!--       Use wall function within constant-flux layer
-!--       Note, grid points listed in bc_h are not included in any calculations in RANS mode and
-!--       are therefore not set here.
-!
-!--       Upward-facing surfaces
-!--       Default surfaces
-          DO  m = 1, surf_def_h(0)%ns
-             i = surf_def_h(0)%i(m)
-             j = surf_def_h(0)%j(m)
-             k = surf_def_h(0)%k(m)
-             e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
-          ENDDO
-!
-!--       Natural surfaces
-          DO  m = 1, surf_lsm_h%ns
-             i = surf_lsm_h%i(m)
-             j = surf_lsm_h%j(m)
-             k = surf_lsm_h%k(m)
-             e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2
-          ENDDO
-!
-!--       Urban surfaces
-          DO  m = 1, surf_usm_h%ns
-             i = surf_usm_h%i(m)
-             j = surf_usm_h%j(m)
-             k = surf_usm_h%k(m)
-             e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2
-          ENDDO
-!
-!--       Vertical surfaces
-          DO  l = 0, 3
-!
-!--          Default surfaces
-             DO  m = 1, surf_def_v(l)%ns
-                i = surf_def_v(l)%i(m)
-                j = surf_def_v(l)%j(m)
-                k = surf_def_v(l)%k(m)
-                e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
-             ENDDO
-!
-!--          Natural surfaces
-             DO  m = 1, surf_lsm_v(l)%ns
-                i = surf_lsm_v(l)%i(m)
-                j = surf_lsm_v(l)%j(m)
-                k = surf_lsm_v(l)%k(m)
-                e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2
-             ENDDO
-!
-!--          Urban surfaces
-             DO  m = 1, surf_usm_v(l)%ns
-                i = surf_usm_v(l)%i(m)
-                j = surf_usm_v(l)%j(m)
-                k = surf_usm_v(l)%k(m)
-                e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2
-             ENDDO
-          ENDDO
-       ENDIF
-
-       IF ( .NOT. child_domain )  THEN
-          !$ACC KERNELS PRESENT(e_p)
-          e_p(nzt+1,:,:) = e_p(nzt,:,:)
-          !$ACC END KERNELS
-!
-!--    Nesting case: if parent operates in RANS mode and child in LES mode,
-!--    no TKE is transfered. This case, set Neumann conditions at lateral and 
-!--    top child boundaries. 
-!--    If not ( both either in RANS or in LES mode ), TKE boundary condition
-!--    is treated in the nesting. 
-       ELSE
-
-          IF ( rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
-
-             e_p(nzt+1,:,:) = e_p(nzt,:,:)
-             IF ( bc_dirichlet_l )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
-             IF ( bc_dirichlet_r )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
-             IF ( bc_dirichlet_s )  e_p(:,nys-1,:) = e_p(:,nys,:)
-             IF ( bc_dirichlet_n )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
-
-          ENDIF
-       ENDIF
-    ENDIF
-
-!
-!-- Boundary conditions for TKE dissipation rate. 
-    IF ( rans_tke_e )  THEN
-!
-!--    Use wall function within constant-flux layer
-!--    Upward-facing surfaces
-!--    Default surfaces
-       DO  m = 1, surf_def_h(0)%ns
-          i = surf_def_h(0)%i(m)
-          j = surf_def_h(0)%j(m)
-          k = surf_def_h(0)%k(m)
-          diss_p(k,j,i) = surf_def_h(0)%us(m)**3          &
-                        / ( kappa * surf_def_h(0)%z_mo(m) )
-       ENDDO
-!
-!--    Natural surfaces
-       DO  m = 1, surf_lsm_h%ns
-          i = surf_lsm_h%i(m)
-          j = surf_lsm_h%j(m)
-          k = surf_lsm_h%k(m)
-          diss_p(k,j,i) = surf_lsm_h%us(m)**3          &
-                        / ( kappa * surf_lsm_h%z_mo(m) )
-       ENDDO
-!
-!--    Urban surfaces
-       DO  m = 1, surf_usm_h%ns
-          i = surf_usm_h%i(m)
-          j = surf_usm_h%j(m)
-          k = surf_usm_h%k(m)
-          diss_p(k,j,i) = surf_usm_h%us(m)**3          &
-                        / ( kappa * surf_usm_h%z_mo(m) )
-       ENDDO
-!
-!--    Vertical surfaces
-       DO  l = 0, 3
-!
-!--       Default surfaces
-          DO  m = 1, surf_def_v(l)%ns
-             i = surf_def_v(l)%i(m)
-             j = surf_def_v(l)%j(m)
-             k = surf_def_v(l)%k(m)
-             diss_p(k,j,i) = surf_def_v(l)%us(m)**3          &
-                           / ( kappa * surf_def_v(l)%z_mo(m) )
-          ENDDO
-!
-!--       Natural surfaces
-          DO  m = 1, surf_lsm_v(l)%ns
-             i = surf_lsm_v(l)%i(m)
-             j = surf_lsm_v(l)%j(m)
-             k = surf_lsm_v(l)%k(m)
-             diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3          &
-                           / ( kappa * surf_lsm_v(l)%z_mo(m) )
-          ENDDO
-!
-!--       Urban surfaces
-          DO  m = 1, surf_usm_v(l)%ns
-             i = surf_usm_v(l)%i(m)
-             j = surf_usm_v(l)%j(m)
-             k = surf_usm_v(l)%k(m)
-             diss_p(k,j,i) = surf_usm_v(l)%us(m)**3          &
-                           / ( kappa * surf_usm_v(l)%z_mo(m) )
-          ENDDO
-       ENDDO
-!
-!--    Limit change of diss to be between -90% and +100%. Also, set an absolute
-!--    minimum value
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb, nzt+1
-                diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i),          &
-                                          2.0_wp * diss(k,j,i) ), &
-                                     0.1_wp * diss(k,j,i),        &
-                                     0.0001_wp )
-             ENDDO
-          ENDDO
-       ENDDO
-
-       IF ( .NOT. child_domain )  THEN
-          diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
-       ENDIF
-    ENDIF
-
 !
 !-- Boundary conditions for salinity
@@ -586,8 +388,4 @@
 !--    given.
        DO  l = 0, 1
-!
-!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
-!--       for downward-facing surfaces at topography bottom (k+1). 
-          kb = MERGE( -1, 1, l == 0 )
           !$OMP PARALLEL DO PRIVATE( i, j, k )
           DO  m = 1, bc_h(l)%ns
@@ -595,5 +393,5 @@
              j = bc_h(l)%j(m)
              k = bc_h(l)%k(m)
-             sa_p(k+kb,j,i) = sa_p(k,j,i)
+             sa_p(k+bc_h(l)%koff,j,i) = sa_p(k,j,i)
           ENDDO
        ENDDO
@@ -620,8 +418,4 @@
 
           DO  l = 0, 1
-!
-!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
-!--          for downward-facing surfaces at topography bottom (k+1). 
-             kb = MERGE( -1, 1, l == 0 )
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, bc_h(l)%ns
@@ -629,5 +423,5 @@
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
-                q_p(k+kb,j,i) = q(k+kb,j,i)
+                q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i)
              ENDDO
           ENDDO
@@ -636,8 +430,4 @@
           
           DO  l = 0, 1
-!
-!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
-!--          for downward-facing surfaces at topography bottom (k+1). 
-             kb = MERGE( -1, 1, l == 0 )
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, bc_h(l)%ns
@@ -645,5 +435,5 @@
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
-                q_p(k+kb,j,i) = q_p(k,j,i)
+                q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i)
              ENDDO
           ENDDO
@@ -662,12 +452,14 @@
 !--       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 k-1
+!--       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(0)%ns
-             i = bc_h(0)%i(m)            
-             j = bc_h(0)%j(m)
-             k = bc_h(0)%k(m)
-             qc_p(k-1,j,i) = 0.0_wp
-             nc_p(k-1,j,i) = 0.0_wp
+             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
 !
@@ -683,12 +475,14 @@
 !--       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 k-1
+!--       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(0)%ns
-             i = bc_h(0)%i(m)            
-             j = bc_h(0)%j(m)
-             k = bc_h(0)%k(m)
-             qr_p(k-1,j,i) = 0.0_wp
-             nr_p(k-1,j,i) = 0.0_wp
+             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
 !
@@ -711,8 +505,4 @@
           
           DO  l = 0, 1
-!
-!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
-!--          for downward-facing surfaces at topography bottom (k+1). 
-             kb = MERGE( -1, 1, l == 0 )
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, bc_h(l)%ns
@@ -720,5 +510,5 @@
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
-                s_p(k+kb,j,i) = s(k+kb,j,i)
+                s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
              ENDDO
           ENDDO
@@ -727,8 +517,4 @@
           
           DO  l = 0, 1
-!
-!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
-!--          for downward-facing surfaces at topography bottom (k+1). 
-             kb = MERGE( -1, 1, l == 0 )
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, bc_h(l)%ns
@@ -736,5 +522,5 @@
                 j = bc_h(l)%j(m)
                 k = bc_h(l)%k(m)
-                s_p(k+kb,j,i) = s_p(k,j,i)
+                s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i)
              ENDDO
           ENDDO
@@ -750,5 +536,8 @@
        ENDIF
 
-    ENDIF    
+    ENDIF  
+!
+!-- Set boundary conditions for subgrid TKE and dissipation (RANS mode)
+    CALL tcm_boundary_conds
 !
 !-- Top/bottom boundary conditions for chemical species
@@ -760,15 +549,8 @@
 !-- version) these levels are handled as a prognostic level, boundary values
 !-- have to be restored here.
-!-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
     IF ( bc_dirichlet_s )  THEN
        v_p(:,nys,:) = v_p(:,nys-1,:)
-       IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
-    ELSEIF ( bc_dirichlet_n )  THEN
-       IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
     ELSEIF ( bc_dirichlet_l ) THEN
        u_p(:,:,nxl) = u_p(:,:,nxl-1)
-       IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
-    ELSEIF ( bc_dirichlet_r )  THEN
-       IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
     ENDIF
 
@@ -777,4 +559,6 @@
 !-- 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
@@ -787,9 +571,9 @@
 
 !
-!-- Lateral boundary conditions for scalar quantities at the outflow
+!-- 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 ( .NOT. constant_diffusion )  e_p(:,nys-1,:) = e_p(:,nys,:)
-       IF ( rans_tke_e )  diss_p(:,nys-1,:) = diss_p(:,nys,:)
        IF ( humidity )  THEN
           q_p(:,nys-1,:) = q_p(:,nys,:)
@@ -806,6 +590,4 @@
     ELSEIF ( bc_radiation_n )  THEN
        pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
-       IF ( .NOT. constant_diffusion )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
-       IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:)
        IF ( humidity )  THEN
           q_p(:,nyn+1,:) = q_p(:,nyn,:)
@@ -822,6 +604,4 @@
     ELSEIF ( bc_radiation_l )  THEN
        pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
-       IF ( .NOT. constant_diffusion )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
-       IF ( rans_tke_e )  diss_p(:,:,nxl-1) = diss_p(:,:,nxl)
        IF ( humidity )  THEN
           q_p(:,:,nxl-1) = q_p(:,:,nxl)
@@ -838,6 +618,4 @@
     ELSEIF ( bc_radiation_r )  THEN
        pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
-       IF ( .NOT. constant_diffusion )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
-       IF ( rans_tke_e )  diss_p(:,:,nxr+1) = diss_p(:,:,nxr)
        IF ( humidity )  THEN
           q_p(:,:,nxr+1) = q_p(:,:,nxr)
Index: palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4101)
+++ palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4102)
@@ -27,4 +27,8 @@
 ! -----------------
 ! $Id$
+! Slightly revise setting of boundary conditions at horizontal walls, use 
+! data-structure offset index instead of pre-calculate it for each facing
+! 
+! 4080 2019-07-09 18:17:37Z suehring
 ! Restore accidantly removed limitation to positive values
 ! 
@@ -804,5 +808,4 @@
     INTEGER(iwp) ::  j                            !< grid index y direction.
     INTEGER(iwp) ::  k                            !< grid index z direction.
-    INTEGER(iwp) ::  kb                           !< variable to set respective boundary value, depends on facing. 
     INTEGER(iwp) ::  l                            !< running index boundary type, for up- and downward-facing walls.
     INTEGER(iwp) ::  m                            !< running index surface elements.
@@ -839,15 +842,8 @@
        CASE ( 'set_bc_bottomtop' )                   
 !
-!--          Bottom boundary condtions for chemical species      
+!--       Boundary condtions for chemical species at horizontal walls      
           DO  lsp = 1, nspec                                                      
-             IF ( ibc_cs_b == 0 )  THEN                   
+             IF ( ibc_cs_b == 0 )  THEN
                 DO  l = 0, 1 
-!
-!--                Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e.
-!--                the chem_species(nspec)%conc_p value at the topography top (k-1)
-!--                is set; for downward-facing surfaces (l=1), kb=1, i.e. the 
-!--                value at the topography bottom (k+1) is set. 
-
-                   kb = MERGE( -1, 1, l == 0 )
                    !$OMP PARALLEL DO PRIVATE( i, j, k )
                    DO  m = 1, bc_h(l)%ns
@@ -855,5 +851,6 @@
                        j = bc_h(l)%j(m)
                        k = bc_h(l)%k(m)
-                      chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i) 
+                      chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
+                                      chem_species(lsp)%conc(k+bc_h(l)%koff,j,i) 
                    ENDDO                                        
                 ENDDO                                       
@@ -863,5 +860,4 @@
 !--             in boundary_conds there is som extra loop over m here for passive tracer
                 DO  l = 0, 1
-                   kb = MERGE( -1, 1, l == 0 )
                    !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
                    DO m = 1, bc_h(l)%ns
@@ -869,5 +865,6 @@
                       j = bc_h(l)%j(m)
                       k = bc_h(l)%k(m)
-                      chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i)
+                      chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
+                                         chem_species(lsp)%conc_p(k,j,i)
 
                    ENDDO
Index: palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- palm/trunk/SOURCE/salsa_mod.f90	(revision 4101)
+++ palm/trunk/SOURCE/salsa_mod.f90	(revision 4102)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Slightly revise setting of boundary conditions at horizontal walls, use 
+! data-structure offset index instead of pre-calculate it for each facing
+! 
+! 4079 2019-07-09 18:04:41Z suehring
 ! Application of monotonic flux limiter for the vertical scalar advection 
 ! up to the topography top (only for the cache-optimized version at the 
@@ -8017,5 +8021,4 @@
     INTEGER(iwp) ::  j    !< grid index y direction
     INTEGER(iwp) ::  k    !< grid index y direction
-    INTEGER(iwp) ::  kb   !< variable to set respective boundary value, depends on facing.
     INTEGER(iwp) ::  l    !< running index boundary type, for up- and downward-facing walls
     INTEGER(iwp) ::  m    !< running index surface elements
@@ -8030,8 +8033,4 @@
 !--       belongs to the atmospheric grid point, therefore, set s_p at k-1
           DO  l = 0, 1
-!
-!--          Set kb, for upward-facing surfaces value at topography top (k-1) is 
-!--          set, for downward-facing surfaces at topography bottom (k+1)
-             kb = MERGE ( -1, 1, l == 0 )
              !$OMP PARALLEL PRIVATE( ib, ic, icc, ig, i, j, k )
              !$OMP DO
@@ -8043,13 +8042,16 @@
 
                 DO  ib = 1, nbins_aerosol
-                   aerosol_number(ib)%conc_p(k+kb,j,i) = aerosol_number(ib)%conc(k+kb,j,i)
+                   aerosol_number(ib)%conc_p(k+bc_h(l)%koff,j,i) =             &
+                                    aerosol_number(ib)%conc(k+bc_h(l)%koff,j,i)
                    DO  ic = 1, ncomponents_mass
                       icc = ( ic - 1 ) * nbins_aerosol + ib
-                      aerosol_mass(icc)%conc_p(k+kb,j,i) = aerosol_mass(icc)%conc(k+kb,j,i)
+                      aerosol_mass(icc)%conc_p(k+bc_h(l)%koff,j,i) =           &
+                                    aerosol_mass(icc)%conc(k+bc_h(l)%koff,j,i)
                    ENDDO
                 ENDDO
                 IF ( .NOT. salsa_gases_from_chem )  THEN
                    DO  ig = 1, ngases_salsa
-                      salsa_gas(ig)%conc_p(k+kb,j,i) = salsa_gas(ig)%conc(k+kb,j,i)
+                      salsa_gas(ig)%conc_p(k+bc_h(l)%koff,j,i) =               &
+                                    salsa_gas(ig)%conc(k+bc_h(l)%koff,j,i)
                    ENDDO
                 ENDIF
@@ -8063,8 +8065,4 @@
 
           DO l = 0, 1
-!
-!--          Set kb, for upward-facing surfaces value at topography top (k-1) is 
-!--          set, for downward-facing surfaces at topography bottom (k+1)
-             kb = MERGE( -1, 1, l == 0 )
              !$OMP PARALLEL PRIVATE( ib, ic, icc, ig, i, j, k )
              !$OMP DO
@@ -8076,13 +8074,16 @@
 
                 DO  ib = 1, nbins_aerosol
-                   aerosol_number(ib)%conc_p(k+kb,j,i) = aerosol_number(ib)%conc_p(k,j,i)
+                   aerosol_number(ib)%conc_p(k+bc_h(l)%koff,j,i) =             &
+                                               aerosol_number(ib)%conc_p(k,j,i)
                    DO  ic = 1, ncomponents_mass
                       icc = ( ic - 1 ) * nbins_aerosol + ib
-                      aerosol_mass(icc)%conc_p(k+kb,j,i) = aerosol_mass(icc)%conc_p(k,j,i)
+                      aerosol_mass(icc)%conc_p(k+bc_h(l)%koff,j,i) =           &
+                                               aerosol_mass(icc)%conc_p(k,j,i)
                    ENDDO
                 ENDDO
                 IF ( .NOT. salsa_gases_from_chem ) THEN
                    DO  ig = 1, ngases_salsa
-                      salsa_gas(ig)%conc_p(k+kb,j,i) = salsa_gas(ig)%conc_p(k,j,i)
+                      salsa_gas(ig)%conc_p(k+bc_h(l)%koff,j,i) =               &
+                                               salsa_gas(ig)%conc_p(k,j,i)
                    ENDDO
                 ENDIF
Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 4101)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 4102)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! - Revise initialization of the boundary data structure
+! - Add new data structure to set boundary conditions at vertical walls
+! 
+! 3943 2019-05-02 09:50:41Z maronga
 ! Removed qsws_eb as it is no longer needed.
 ! 
@@ -306,10 +310,12 @@
 !-- are applied 
     TYPE bc_type
-
-       INTEGER(iwp) :: ns                                  !< number of surface elements on the PE
-
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< x-index linking to the PALM 3D-grid  
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< y-index linking to the PALM 3D-grid    
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< z-index linking to the PALM 3D-grid   
+       INTEGER(iwp) :: ioff !< offset value in x-direction, used to determine index of surface element
+       INTEGER(iwp) :: joff !< offset value in y-direction, used to determine index of surface element
+       INTEGER(iwp) :: koff !< offset value in z-direction, used to determine index of surface element
+       INTEGER(iwp) :: ns   !< number of surface elements on the PE
+
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i !< x-index linking to the PALM 3D-grid  
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j !< y-index linking to the PALM 3D-grid    
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k !< z-index linking to the PALM 3D-grid   
 
        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i) 
@@ -600,4 +606,5 @@
 
     TYPE (bc_type), DIMENSION(0:1)           ::  bc_h        !< boundary condition data type, horizontal upward- and downward facing surfaces
+    TYPE (bc_type), DIMENSION(0:3)           ::  bc_v        !< boundary condition data type, vertical surfaces
 
     TYPE (surf_type), DIMENSION(0:2), TARGET ::  surf_def_h  !< horizontal default surfaces (Up, Down, and Top)
@@ -668,6 +675,6 @@
 !
 !-- Public variables
-    PUBLIC bc_h, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file, ns_v_on_file, surf_def_h, &
-           surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type,                  &
+    PUBLIC bc_h, bc_v, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file, ns_v_on_file,       &
+           surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type,      &
            vertical_surfaces_exist, surf_bulk_cloud_model, surf_microphysics_morrison,             &
            surf_microphysics_seifert
@@ -687,5 +694,6 @@
 ! Description:
 ! ------------
-!> Initialize data type for setting boundary conditions at horizontal surfaces. 
+!> Initialize data type for setting boundary conditions at horizontal and  
+!> vertical surfaces. 
 !------------------------------------------------------------------------------!
     SUBROUTINE init_bc
@@ -696,92 +704,175 @@
        INTEGER(iwp) ::  j         !< loop index along y-direction
        INTEGER(iwp) ::  k         !< loop index along y-direction
-
-       INTEGER(iwp), DIMENSION(0:1) ::  num_h         !< number of surfaces
-       INTEGER(iwp), DIMENSION(0:1) ::  num_h_kji     !< number of surfaces for (j,i)-grid point
-       INTEGER(iwp), DIMENSION(0:1) ::  start_index_h !< local start index
-
-!
-!--    First of all, count the number of upward- and downward-facing surfaces
-       num_h = 0
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             DO  k = nzb+1, nzt
-!
-!--             Check if current gridpoint belongs to the atmosphere
-                IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
-!
-!--                Upward-facing
-                   IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )              &
-                      num_h(0) = num_h(0) + 1
-!
-!--                Downward-facing
-                   IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )              &
-                      num_h(1) = num_h(1) + 1
-                ENDIF
+       INTEGER(iwp) ::  l         !< running index for differently aligned surfaces
+
+       INTEGER(iwp), DIMENSION(0:1) ::  num_h         !< number of horizontal surfaces on subdomain
+       INTEGER(iwp), DIMENSION(0:1) ::  num_h_kji     !< number of horizontal surfaces at (j,i)-grid point
+       INTEGER(iwp), DIMENSION(0:1) ::  start_index_h !< local start index of horizontal surface elements
+       
+       INTEGER(iwp), DIMENSION(0:3) ::  num_v         !< number of vertical surfaces on subdomain
+       INTEGER(iwp), DIMENSION(0:3) ::  num_v_kji     !< number of vertical surfaces at (j,i)-grid point
+       INTEGER(iwp), DIMENSION(0:3) ::  start_index_v !< local start index of vertical surface elements
+!
+!--    Set offset indices, i.e. index difference between surface element and 
+!--    surface-bounded grid point.
+!--    Horizontal surfaces - no horizontal offsets
+       bc_h(:)%ioff = 0
+       bc_h(:)%joff = 0
+!
+!--    Horizontal surfaces, upward facing (0) and downward facing (1)
+       bc_h(0)%koff   = -1
+       bc_h(1)%koff   = 1
+!
+!--    Vertical surfaces - no vertical offset
+       bc_v(0:3)%koff = 0
+!
+!--    North- and southward facing - no offset in x
+       bc_v(0:1)%ioff = 0
+!
+!--    Northward facing offset in y
+       bc_v(0)%joff = -1
+!
+!--    Southward facing offset in y
+       bc_v(1)%joff = 1
+!
+!--    East- and westward facing - no offset in y
+       bc_v(2:3)%joff = 0
+!
+!--    Eastward facing offset in x
+       bc_v(2)%ioff = -1
+!
+!--    Westward facing offset in y
+       bc_v(3)%ioff = 1
+!
+!--    Initialize data structure for horizontal surfaces, i.e. count the number
+!--    of surface elements, allocate and initialize the respective index arrays, 
+!--    and set the respective start and end indices at each (j,i)-location.
+       DO  l = 0, 1
+!
+!--       Count the number of upward- and downward-facing surfaces on subdomain
+          num_h(l) = 0
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                DO  k = nzb+1, nzt
+!         
+!--                Check if current gridpoint belongs to the atmosphere
+                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
+                      IF ( .NOT. BTEST( wall_flags_0(k+bc_h(l)%koff,           &
+                                                     j+bc_h(l)%joff,           &
+                                                     i+bc_h(l)%ioff), 0 ) )    &
+                         num_h(l) = num_h(l) + 1
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
+!         
+!--       Save the number of horizontal surface elements
+          bc_h(l)%ns = num_h(l)
+!
+!--       ALLOCATE arrays for horizontal surfaces
+          ALLOCATE( bc_h(l)%i(1:bc_h(l)%ns) )
+          ALLOCATE( bc_h(l)%j(1:bc_h(l)%ns) )
+          ALLOCATE( bc_h(l)%k(1:bc_h(l)%ns) )
+          ALLOCATE( bc_h(l)%start_index(nysg:nyng,nxlg:nxrg) )
+          ALLOCATE( bc_h(l)%end_index(nysg:nyng,nxlg:nxrg)   )
+          bc_h(l)%start_index = 1
+          bc_h(l)%end_index   = 0
+          
+          num_h(l)         = 1
+          start_index_h(l) = 1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+          
+                num_h_kji(l) = 0
+                DO  k = nzb+1, nzt
+!         
+!--                Check if current gridpoint belongs to the atmosphere
+                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
+!         
+!--                   Upward-facing
+                      IF ( .NOT. BTEST( wall_flags_0(k+bc_h(l)%koff,           &
+                                                     j+bc_h(l)%joff,           &
+                                                     i+bc_h(l)%ioff), 0 )      &
+                         )  THEN
+                         bc_h(l)%i(num_h(l)) = i
+                         bc_h(l)%j(num_h(l)) = j
+                         bc_h(l)%k(num_h(l)) = k
+                         num_h_kji(l)        = num_h_kji(l) + 1
+                         num_h(l)            = num_h(l) + 1
+                      ENDIF
+                   ENDIF
+                ENDDO
+                bc_h(l)%start_index(j,i) = start_index_h(l)
+                bc_h(l)%end_index(j,i)   = bc_h(l)%start_index(j,i) +          &
+                                           num_h_kji(l) - 1
+                start_index_h(l)         = bc_h(l)%end_index(j,i) + 1
              ENDDO
           ENDDO
        ENDDO
-!
-!--    Save the number of surface elements
-       bc_h(0)%ns = num_h(0)
-       bc_h(1)%ns = num_h(1)
-!
-!--    ALLOCATE data type variables
-!--    Upward facing
-       ALLOCATE( bc_h(0)%i(1:bc_h(0)%ns) )
-       ALLOCATE( bc_h(0)%j(1:bc_h(0)%ns) )
-       ALLOCATE( bc_h(0)%k(1:bc_h(0)%ns) )
-       ALLOCATE( bc_h(0)%start_index(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE( bc_h(0)%end_index(nysg:nyng,nxlg:nxrg)   )
-       bc_h(0)%start_index = 1
-       bc_h(0)%end_index   = 0
-!
-!--    Downward facing
-       ALLOCATE( bc_h(1)%i(1:bc_h(1)%ns) )
-       ALLOCATE( bc_h(1)%j(1:bc_h(1)%ns) )
-       ALLOCATE( bc_h(1)%k(1:bc_h(1)%ns) )
-       ALLOCATE( bc_h(1)%start_index(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE( bc_h(1)%end_index(nysg:nyng,nxlg:nxrg)   )
-       bc_h(1)%start_index = 1
-       bc_h(1)%end_index   = 0
-!
-!--    Store the respective indices on data type
-       num_h(0:1)         = 1
-       start_index_h(0:1) = 1
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-
-             num_h_kji(0:1) = 0
-             DO  k = nzb+1, nzt
-!
-!--             Check if current gridpoint belongs to the atmosphere
-                IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
-!
-!--                Upward-facing
-                   IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )  THEN
-                      bc_h(0)%i(num_h(0)) = i
-                      bc_h(0)%j(num_h(0)) = j
-                      bc_h(0)%k(num_h(0)) = k
-                      num_h_kji(0)        = num_h_kji(0) + 1
-                      num_h(0)            = num_h(0) + 1
+
+!
+!--    Initialize data structure for vertical surfaces, i.e. count the number
+!--    of surface elements, allocate and initialize the respective index arrays, 
+!--    and set the respective start and end indices at each (j,i)-location.
+       DO  l = 0, 3
+!
+!--       Count the number of upward- and downward-facing surfaces on subdomain
+          num_v(l) = 0
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                DO  k = nzb+1, nzt
+!         
+!--                Check if current gridpoint belongs to the atmosphere
+                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
+                      IF ( .NOT. BTEST( wall_flags_0(k+bc_v(l)%koff,           &
+                                                     j+bc_v(l)%joff,           &
+                                                     i+bc_v(l)%ioff), 0 ) )    &
+                         num_v(l) = num_v(l) + 1
                    ENDIF
-!
-!--                Downward-facing
-                   IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )  THEN
-                      bc_h(1)%i(num_h(1)) = i
-                      bc_h(1)%j(num_h(1)) = j
-                      bc_h(1)%k(num_h(1)) = k
-                      num_h_kji(1)        = num_h_kji(1) + 1
-                      num_h(1)            = num_h(1) + 1
+                ENDDO
+             ENDDO
+          ENDDO
+!         
+!--       Save the number of horizontal surface elements
+          bc_v(l)%ns = num_v(l)
+!
+!--       ALLOCATE arrays for horizontal surfaces
+          ALLOCATE( bc_v(l)%i(1:bc_v(l)%ns) )
+          ALLOCATE( bc_v(l)%j(1:bc_v(l)%ns) )
+          ALLOCATE( bc_v(l)%k(1:bc_v(l)%ns) )
+          ALLOCATE( bc_v(l)%start_index(nysg:nyng,nxlg:nxrg) )
+          ALLOCATE( bc_v(l)%end_index(nysg:nyng,nxlg:nxrg)   )
+          bc_v(l)%start_index = 1
+          bc_v(l)%end_index   = 0
+          
+          num_v(l)         = 1
+          start_index_v(l) = 1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+          
+                num_v_kji(l) = 0
+                DO  k = nzb+1, nzt
+!         
+!--                Check if current gridpoint belongs to the atmosphere
+                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
+!         
+!--                   Upward-facing
+                      IF ( .NOT. BTEST( wall_flags_0(k+bc_v(l)%koff,           &
+                                                     j+bc_v(l)%joff,           &
+                                                     i+bc_v(l)%ioff), 0 )      &
+                         )  THEN
+                         bc_v(l)%i(num_v(l)) = i
+                         bc_v(l)%j(num_v(l)) = j
+                         bc_v(l)%k(num_v(l)) = k
+                         num_v_kji(l)        = num_v_kji(l) + 1
+                         num_v(l)            = num_v(l) + 1
+                      ENDIF
                    ENDIF
-                ENDIF
+                ENDDO
+                bc_v(l)%start_index(j,i) = start_index_v(l)
+                bc_v(l)%end_index(j,i)   = bc_v(l)%start_index(j,i) +          &
+                                           num_v_kji(l) - 1
+                start_index_v(l)         = bc_v(l)%end_index(j,i) + 1
              ENDDO
-             bc_h(0)%start_index(j,i) = start_index_h(0)
-             bc_h(0)%end_index(j,i)   = bc_h(0)%start_index(j,i) + num_h_kji(0) - 1
-             start_index_h(0)         = bc_h(0)%end_index(j,i) + 1
-
-             bc_h(1)%start_index(j,i) = start_index_h(1)
-             bc_h(1)%end_index(j,i)   = bc_h(1)%start_index(j,i) + num_h_kji(1) - 1
-             start_index_h(1)         = bc_h(1)%end_index(j,i) + 1
           ENDDO
        ENDDO
Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4101)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4102)
@@ -25,4 +25,13 @@
 ! -----------------
 ! $Id$
+! - Modularize setting of boundary conditions for TKE and dissipation
+! - Neumann boundary condition for TKE at model top is set also in child domain
+! - Revise setting of Neumann boundary conditions at non-cyclic lateral 
+!   boundaries
+! - Bugfix, set Neumann boundary condition for TKE at vertical wall instead of 
+!   an implicit Dirichlet boundary condition which implied a sink of TKE
+!   at vertical walls
+! 
+! 4048 2019-06-21 21:00:21Z knoop
 ! write out preprocessor directives; remove tailing whitespaces
 !
@@ -192,58 +201,77 @@
 !> @todo Check for random disturbances
 !> @note <Enter notes on the module>
-!------------------------------------------------------------------------------!
+!-----------------------------------------------------------------------------!
  MODULE turbulence_closure_mod
 
 
-    USE arrays_3d,                                                             &
-        ONLY:  diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3,    &
-               e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m,           &
+    USE arrays_3d,                                                            &
+        ONLY:  diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3,   &
+               e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m,          &
                te_m, tend, u, v, vpt, w
 
-    USE basic_constants_and_equations_mod,                                     &
+    USE basic_constants_and_equations_mod,                                    &
         ONLY:  g, kappa, lv_d_cp, lv_d_rd, rd_d_rv
 
-    USE control_parameters,                                                    &
-        ONLY:  constant_diffusion, dt_3d, e_init, humidity,                    &
-               initializing_actions, intermediate_timestep_count,              &
-               intermediate_timestep_count_max, km_constant,                   &
-               les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number,  &
-               pt_reference, rans_mode, rans_tke_e, rans_tke_l,                &
-               timestep_scheme, turbulence_closure,                            &
-               turbulent_inflow, use_upstream_for_tke, vpt_reference,          &
+    USE control_parameters,                                                   &
+        ONLY:  bc_dirichlet_l,                                                &
+               bc_dirichlet_n,                                                &
+               bc_dirichlet_r,                                                &
+               bc_dirichlet_s,                                                &
+               bc_radiation_l,                                                &
+               bc_radiation_n,                                                &
+               bc_radiation_r,                                                &
+               bc_radiation_s,                                                &
+               child_domain,                                                  &
+               constant_diffusion, dt_3d, e_init, humidity,                   &
+               initializing_actions, intermediate_timestep_count,             &
+               intermediate_timestep_count_max, km_constant,                  &
+               les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, &
+               pt_reference, rans_mode, rans_tke_e, rans_tke_l,               &
+               timestep_scheme, turbulence_closure,                           &
+               turbulent_inflow, use_upstream_for_tke, vpt_reference,         &
                ws_scheme_sca, current_timestep_number
 
-    USE advec_ws,                                                              &
+    USE advec_ws,                                                             &
         ONLY:  advec_s_ws
 
-    USE advec_s_bc_mod,                                                        &
+    USE advec_s_bc_mod,                                                       &
         ONLY:  advec_s_bc
 
-    USE advec_s_pw_mod,                                                        &
+    USE advec_s_pw_mod,                                                       &
         ONLY:  advec_s_pw
 
-    USE advec_s_up_mod,                                                        &
+    USE advec_s_up_mod,                                                       &
         ONLY:  advec_s_up
 
-    USE cpulog,                                                                &
+    USE cpulog,                                                               &
         ONLY:  cpu_log, log_point_s
 
-    USE indices,                                                               &
-        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,     &
+    USE indices,                                                              &
+        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,    &
                wall_flags_0
 
     USE kinds
 
-    USE ocean_mod,                                                             &
+    USE ocean_mod,                                                            &
         ONLY:  prho_reference
 
     USE pegrid
 
-    USE plant_canopy_model_mod,                                                &
+    USE plant_canopy_model_mod,                                               &
         ONLY:  pcm_tendency
 
-    USE statistics,                                                            &
+    USE statistics,                                                           &
         ONLY:  hom, hom_sum, statistic_regions
-
+        
+    USE surface_mod,                                                          &
+        ONLY:  bc_h,                                                          &
+               bc_v,                                                          &
+               get_topography_top_index_ji,                                   &
+               surf_def_h,                                                    &
+               surf_def_v,                                                    &
+               surf_lsm_h,                                                    &
+               surf_lsm_v,                                                    &
+               surf_usm_h,                                                    &
+               surf_usm_v
 
     IMPLICIT NONE
@@ -271,5 +299,6 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  l_wall !< near-wall mixing length
 
-
+!
+!-- Public variables
     PUBLIC c_0, rans_const_c, rans_const_sigma
 
@@ -277,22 +306,28 @@
 
     PRIVATE
-
-
-    PUBLIC &
-       tcm_check_parameters, &
-       tcm_check_data_output, &
-       tcm_define_netcdf_grid, &
-       tcm_init_arrays, &
-       tcm_init, &
-       tcm_actions, &
-       tcm_prognostic_equations, &
-       tcm_swap_timelevel, &
-       tcm_3d_data_averaging, &
-       tcm_data_output_2d, &
-       tcm_data_output_3d, &
+!
+!-- Public subroutines
+    PUBLIC                                                                     &
+       tcm_boundary_conds,                                                     &
+       tcm_check_parameters,                                                   &
+       tcm_check_data_output,                                                  &
+       tcm_define_netcdf_grid,                                                 &
+       tcm_init_arrays,                                                        &
+       tcm_init,                                                               &
+       tcm_actions,                                                            &
+       tcm_prognostic_equations,                                               &
+       tcm_swap_timelevel,                                                     &
+       tcm_3d_data_averaging,                                                  &
+       tcm_data_output_2d,                                                     &
+       tcm_data_output_3d,                                                     &
        tcm_diffusivities
 
 !
 !-- PALM interfaces:
+!-- Boundary conditions for subgrid TKE and dissipation
+    INTERFACE tcm_boundary_conds
+       MODULE PROCEDURE tcm_boundary_conds
+    END INTERFACE tcm_boundary_conds
+!
 !-- Input parameter checks to be done in check_parameters
     INTERFACE tcm_check_parameters
@@ -371,4 +406,249 @@
  CONTAINS
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check parameters routine for turbulence closure module.
+!------------------------------------------------------------------------------!
+ SUBROUTINE tcm_boundary_conds
+
+    USE pmc_interface,                                                         &
+        ONLY : rans_mode_parent
+ 
+    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
+!
+!-- Boundary conditions for TKE.
+    IF ( .NOT. constant_diffusion )  THEN
+!
+!--    In LES mode, Neumann conditions with de/x_i=0 are assumed at solid walls.
+!--    Note, only TKE is prognostic in this case and dissipation is only 
+!--    a diagnostic quantity. 
+       IF ( .NOT. rans_mode )  THEN
+!
+!--       Horizontal walls, upward- and downward-facing
+          DO  l = 0, 1
+             !$OMP PARALLEL DO PRIVATE( i, j, k )
+             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+             !$ACC PRESENT(bc_h, e_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)
+                e_p(k+bc_h(l)%koff,j,i) = e_p(k,j,i)
+             ENDDO
+          ENDDO
+!
+!--       Vertical walls
+          DO  l = 0, 3
+             !$OMP PARALLEL DO PRIVATE( i, j, k )
+             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
+             !$ACC PRESENT(bc_v, e_p)
+             DO  m = 1, bc_v(l)%ns
+                i = bc_v(l)%i(m)            
+                j = bc_v(l)%j(m)
+                k = bc_v(l)%k(m)
+                e_p(k,j+bc_v(l)%joff,i+bc_v(l)%ioff) = e_p(k,j,i)
+             ENDDO
+          ENDDO
+!
+!--    In RANS mode, wall function is used as boundary condition for TKE
+       ELSE
+!
+!--       Use wall function within constant-flux layer
+!--       Note, grid points listed in bc_h are not included in any calculations in RANS mode and
+!--       are therefore not set here.
+!
+!--       Upward-facing surfaces
+!--       Default surfaces
+          DO  m = 1, surf_def_h(0)%ns
+             i = surf_def_h(0)%i(m)
+             j = surf_def_h(0)%j(m)
+             k = surf_def_h(0)%k(m)
+             e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
+          ENDDO
+!
+!--       Natural surfaces
+          DO  m = 1, surf_lsm_h%ns
+             i = surf_lsm_h%i(m)
+             j = surf_lsm_h%j(m)
+             k = surf_lsm_h%k(m)
+             e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2
+          ENDDO
+!
+!--       Urban surfaces
+          DO  m = 1, surf_usm_h%ns
+             i = surf_usm_h%i(m)
+             j = surf_usm_h%j(m)
+             k = surf_usm_h%k(m)
+             e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2
+          ENDDO
+!
+!--       Vertical surfaces
+          DO  l = 0, 3
+!
+!--          Default surfaces
+             DO  m = 1, surf_def_v(l)%ns
+                i = surf_def_v(l)%i(m)
+                j = surf_def_v(l)%j(m)
+                k = surf_def_v(l)%k(m)
+                e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
+             ENDDO
+!
+!--          Natural surfaces
+             DO  m = 1, surf_lsm_v(l)%ns
+                i = surf_lsm_v(l)%i(m)
+                j = surf_lsm_v(l)%j(m)
+                k = surf_lsm_v(l)%k(m)
+                e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2
+             ENDDO
+!
+!--          Urban surfaces
+             DO  m = 1, surf_usm_v(l)%ns
+                i = surf_usm_v(l)%i(m)
+                j = surf_usm_v(l)%j(m)
+                k = surf_usm_v(l)%k(m)
+                e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    Set Neumann boundary condition for TKE at model top. Do this also
+!--    in case of a nested run.
+       !$ACC KERNELS PRESENT(e_p)
+       e_p(nzt+1,:,:) = e_p(nzt,:,:)
+       !$ACC END KERNELS
+!
+!--    Nesting case: if parent operates in RANS mode and child in LES mode,
+!--    no TKE is transfered. This case, set Neumann conditions at lateral and 
+!--    top child boundaries. 
+!--    If not ( both either in RANS or in LES mode ), TKE boundary condition
+!--    is treated in the nesting. 
+       If ( child_domain )  THEN
+          IF ( rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
+
+             e_p(nzt+1,:,:) = e_p(nzt,:,:)
+             IF ( bc_dirichlet_l )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
+             IF ( bc_dirichlet_r )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
+             IF ( bc_dirichlet_s )  e_p(:,nys-1,:) = e_p(:,nys,:)
+             IF ( bc_dirichlet_n )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
+
+          ENDIF
+       ENDIF
+!
+!--    At in- and outflow boundaries also set Neumann boundary conditions 
+!--    for the SGS-TKE. An exception is made for the child domain if 
+!--    both parent and child operate in RANS mode. This case no 
+!--    lateral Neumann boundary conditions will be set but Dirichlet
+!--    conditions will be set in the nesting.
+       IF ( .NOT. child_domain  .AND.  .NOT. rans_mode_parent  .AND.           &
+            .NOT. rans_mode )  THEN
+          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
+             e_p(:,nys-1,:) = e_p(:,nys,:)
+             IF ( rans_tke_e )  diss_p(:,nys-1,:) = diss_p(:,nys,:)
+          ENDIF
+          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
+             e_p(:,nyn+1,:) = e_p(:,nyn,:)
+             IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:)  
+          ENDIF
+          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
+             e_p(:,:,nxl-1) = e_p(:,:,nxl)
+             IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:)  
+          ENDIF
+          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
+             e_p(:,:,nxr+1) = e_p(:,:,nxr)
+             IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 
+          ENDIF
+       ENDIF
+    ENDIF
+
+!
+!-- Boundary conditions for TKE dissipation rate in RANS mode. 
+    IF ( rans_tke_e )  THEN
+!
+!--    Use wall function within constant-flux layer
+!--    Upward-facing surfaces
+!--    Default surfaces
+       DO  m = 1, surf_def_h(0)%ns
+          i = surf_def_h(0)%i(m)
+          j = surf_def_h(0)%j(m)
+          k = surf_def_h(0)%k(m)
+          diss_p(k,j,i) = surf_def_h(0)%us(m)**3          &
+                        / ( kappa * surf_def_h(0)%z_mo(m) )
+       ENDDO
+!
+!--    Natural surfaces
+       DO  m = 1, surf_lsm_h%ns
+          i = surf_lsm_h%i(m)
+          j = surf_lsm_h%j(m)
+          k = surf_lsm_h%k(m)
+          diss_p(k,j,i) = surf_lsm_h%us(m)**3          &
+                        / ( kappa * surf_lsm_h%z_mo(m) )
+       ENDDO
+!
+!--    Urban surfaces
+       DO  m = 1, surf_usm_h%ns
+          i = surf_usm_h%i(m)
+          j = surf_usm_h%j(m)
+          k = surf_usm_h%k(m)
+          diss_p(k,j,i) = surf_usm_h%us(m)**3          &
+                        / ( kappa * surf_usm_h%z_mo(m) )
+       ENDDO
+!
+!--    Vertical surfaces
+       DO  l = 0, 3
+!
+!--       Default surfaces
+          DO  m = 1, surf_def_v(l)%ns
+             i = surf_def_v(l)%i(m)
+             j = surf_def_v(l)%j(m)
+             k = surf_def_v(l)%k(m)
+             diss_p(k,j,i) = surf_def_v(l)%us(m)**3          &
+                           / ( kappa * surf_def_v(l)%z_mo(m) )
+          ENDDO
+!
+!--       Natural surfaces
+          DO  m = 1, surf_lsm_v(l)%ns
+             i = surf_lsm_v(l)%i(m)
+             j = surf_lsm_v(l)%j(m)
+             k = surf_lsm_v(l)%k(m)
+             diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3          &
+                           / ( kappa * surf_lsm_v(l)%z_mo(m) )
+          ENDDO
+!
+!--       Urban surfaces
+          DO  m = 1, surf_usm_v(l)%ns
+             i = surf_usm_v(l)%i(m)
+             j = surf_usm_v(l)%j(m)
+             k = surf_usm_v(l)%k(m)
+             diss_p(k,j,i) = surf_usm_v(l)%us(m)**3          &
+                           / ( kappa * surf_usm_v(l)%z_mo(m) )
+          ENDDO
+       ENDDO
+!
+!--    Limit change of diss to be between -90% and +100%. Also, set an absolute
+!--    minimum value
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb, nzt+1
+                diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i),          &
+                                          2.0_wp * diss(k,j,i) ), &
+                                     0.1_wp * diss(k,j,i),        &
+                                     0.0001_wp )
+             ENDDO
+          ENDDO
+       ENDDO
+
+       
+       diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
+     
+    ENDIF
+ END SUBROUTINE tcm_boundary_conds 
+ 
 !------------------------------------------------------------------------------!
 ! Description:
@@ -926,7 +1206,4 @@
     USE model_1d_mod,                                                          &
         ONLY:  e1d, kh1d, km1d
-
-    USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index_ji
 
     IMPLICIT NONE
@@ -1713,7 +1990,4 @@
         ONLY:  phi_m
 
-    USE surface_mod,                                                           &
-        ONLY :  surf_def_h, surf_lsm_h, surf_usm_h
-
     IMPLICIT NONE
 
@@ -2418,8 +2692,4 @@
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  bulk_cloud_model
-
-    USE surface_mod,                                                           &
-        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
-                surf_usm_v
 
     IMPLICIT NONE
@@ -3122,8 +3392,4 @@
         ONLY:  bulk_cloud_model
 
-    USE surface_mod,                                                           &
-        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
-                surf_usm_v
-
     IMPLICIT NONE
 
@@ -3767,7 +4033,4 @@
         ONLY:  use_sgs_for_particles, wang_kernel
 
-    USE surface_mod,                                                           &
-       ONLY :  bc_h
-
     IMPLICIT NONE
 
@@ -3929,5 +4192,7 @@
 
 !
-!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
+!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:).
+!-- Note, bc cannot be set in tcm_boundary conditions as the dissipation 
+!-- in LES mode is only a diagnostic quantity. 
     IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles  .OR.                  &
          wang_kernel  .OR.  collision_turbulence  ) )  THEN
@@ -3973,7 +4238,4 @@
     USE particle_attributes,                                                   &
         ONLY:  use_sgs_for_particles, wang_kernel
-
-    USE surface_mod,                                                           &
-       ONLY :  bc_h
 
     IMPLICIT NONE
@@ -4062,4 +4324,6 @@
 !--    For each surface type determine start and end index (in case of elevated
 !--    topography several up/downward facing surfaces may exist.
+!--    Note, bc cannot be set in tcm_boundary conditions as the dissipation 
+!--    in LES mode is only a diagnostic quantity. 
        surf_s = bc_h(0)%start_index(j,i)
        surf_e = bc_h(0)%end_index(j,i)
@@ -4322,8 +4586,4 @@
         ONLY:  phi_m
 
-    USE surface_mod,                                                           &
-        ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
-                surf_usm_h, surf_usm_v
-
     INTEGER(iwp) ::  i          !< loop index
     INTEGER(iwp) ::  j          !< loop index
Index: palm/trunk/SOURCE/vertical_nesting_mod.f90
===================================================================
--- palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 4101)
+++ palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 4102)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! - Slightly revise setting of boundary conditions at horizontal walls, use 
+!   data-structure offset index instead of pre-calculate it for each facing
+! 
+! 4101 2019-07-17 15:14:26Z gronemeier
 ! remove old_dt
 ! 
@@ -2711,6 +2715,5 @@
         INTEGER(iwp)                          ::  i
         INTEGER(iwp)                          ::  j
-        INTEGER(iwp)                          ::  k
-        INTEGER(iwp)                          ::  kb !< variable to set respective boundary value, depends on facing. 
+        INTEGER(iwp)                          ::  k 
         INTEGER(iwp)                          ::  l  !< running index boundary type, for up- and downward-facing walls
         INTEGER(iwp)                          ::  m  !< running index surface elements
@@ -2860,15 +2863,11 @@
             IF ( ibc_pt_b == 1 )  THEN 
                DO  l = 0, 1 
-!
-!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
-!--       for downward-facing surfaces at topography bottom (k+1). 
-                    kb = MERGE( -1, 1, l == 0 )
-                    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(k+kb,j,i) = pt(k,j,i)
-                    ENDDO
-                 ENDDO
+                  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(k+bc_h(l)%koff,j,i) = pt(k,j,i)
+                  ENDDO
+               ENDDO
             ENDIF
    
