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