Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 4087)
+++ 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)
