Index: /palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/boundary_conds.f90	(revision 3128)
+++ /palm/trunk/SOURCE/boundary_conds.f90	(revision 3129)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! - Use wall function for e_p and diss_p in case of rans_tke_e
+! - move limitation of diss_p from tcm_prognostic
+! 
+! 2938 2018-03-27 15:52:42Z suehring
 ! Set boundary condition for TKE and TKE dissipation rate in case of nesting
 ! and if parent model operates in RANS mode but child model in LES mode.
@@ -182,8 +186,8 @@
     USE arrays_3d,                                                             &
         ONLY:  c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l,  &
-               diss_p, dzu, e_p, nc_p, nr_p, pt, pt_p, q, q_p, qc_p, qr_p, s,  & 
+               diss, diss_p, dzu, e_p, nc_p, nr_p, pt, pt_p, q, q_p, qc_p, qr_p, s,  & 
                s_p, sa, sa_p, u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,  &
                v, vg, 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, pt_init
+               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init, ddzu
 
     USE chemistry_model_mod,                                                   &
@@ -196,5 +200,5 @@
                ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t,         &
                ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r,     &
-               inflow_s, intermediate_timestep_count,                          &
+               inflow_s, intermediate_timestep_count, kappa,                   &
                microphysics_morrison, microphysics_seifert, nest_domain,       &
                nest_bound_l, nest_bound_n, nest_bound_r, nest_bound_s, nudging,&
@@ -217,5 +221,9 @@
 
     USE surface_mod,                                                           &
-        ONLY :  bc_h
+        ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
+                surf_usm_h, surf_usm_v
+
+    USE turbulence_closure_mod,                                                &
+        ONLY:  c_0
 
     IMPLICIT NONE
@@ -344,4 +352,59 @@
              ENDDO
           ENDDO
+       ELSE
+!
+!--       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)
+             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
 
@@ -372,6 +435,82 @@
 !
 !-- Boundary conditions for TKE dissipation rate. 
-    IF ( rans_tke_e .AND. .NOT. nest_domain )  THEN
-       diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
+    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. nest_domain )  THEN
+          diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
+       ENDIF
     ENDIF
 
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 3128)
+++ /palm/trunk/SOURCE/modules.f90	(revision 3129)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! add target attribute to km and kh, necessary for output in tcm_data_output_3d
+! 
+! 3120 2018-07-11 18:30:57Z gronemeier
 ! +les_dynamic
 ! 
@@ -747,6 +750,6 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_v    !< 6th-order advective flux at south face of grid box - v-component
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_w    !< 6th-order advective flux at south face of grid box - w-component
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  kh          !< eddy diffusivity for heat 
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  km          !< eddy diffusivity for momentum
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  kh  !< eddy diffusivity for heat 
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  km  !< eddy diffusivity for momentum
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  prr         !< rain rate
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p_loc       !< local array in multigrid/sor solver containing the pressure which is iteratively advanced in each iteration step
Index: /palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3128)
+++ /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3129)
@@ -25,4 +25,12 @@
 ! -----------------
 ! $Id$
+! - move limitation of diss to boundary_conds
+! - move boundary conditions for e and diss to boundary_conds
+! - consider non-default surfaces in tcm_diffusivities
+! - use z_mo within surface layer instead of calculating it
+! - resort output after case select -> reduced code duplication
+! - when using rans_tke_e and 1d-model, do not use e1d, km1d and diss1d
+! 
+! 3121 2018-07-11 18:46:49Z gronemeier
 ! - created the function phi_m
 ! - implemented km = u* * kappa * zp / phi_m in production_e_init for all
@@ -117,5 +125,4 @@
 !> @todo Check for random disturbances
 !> @note <Enter notes on the module>
-!> @bug  TKE-e closure still crashes due to too small dt
 !------------------------------------------------------------------------------!
  MODULE turbulence_closure_mod
@@ -202,8 +209,9 @@
 
     !> @todo remove debug variables
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: diss_prod1, diss_adve1, diss_diff1, &
-                                               diss_prod2, diss_adve2, diss_diff2, &
-                                               diss_prod3, diss_adve3, diss_diff3, &
-                                               dummy1, dummy2, dummy3
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
+       diss_prod1, diss_adve1, diss_diff1, &
+       diss_prod2, diss_adve2, diss_diff2, &
+       diss_prod3, diss_adve3, diss_diff3, &
+       dummy1, dummy2, dummy3
 
 
@@ -689,13 +697,15 @@
     CHARACTER (LEN=*) ::  variable   !< name of variable
 
-    INTEGER(iwp) ::  av   !< flag for (non-)average output
-    INTEGER(iwp) ::  i    !< loop index
-    INTEGER(iwp) ::  j    !< loop index
-    INTEGER(iwp) ::  k    !< loop index
-    INTEGER(iwp) ::  nzb_do   !< vertical output index (bottom)
-    INTEGER(iwp) ::  nzt_do   !< vertical output index (top)
+    INTEGER(iwp) ::  av        !< flag for (non-)average output
+    INTEGER(iwp) ::  flag_nr   !< number of masking flag
+    INTEGER(iwp) ::  i         !< loop index
+    INTEGER(iwp) ::  j         !< loop index
+    INTEGER(iwp) ::  k         !< loop index
+    INTEGER(iwp) ::  nzb_do    !< vertical output index (bottom)
+    INTEGER(iwp) ::  nzt_do    !< vertical output index (top)
 
     LOGICAL ::  found   !< flag if output variable is found
     LOGICAL ::  two_d   !< flag parameter that indicates 2D variables (horizontal cross sections)
+    LOGICAL ::  resorted  !< flag if output is already resorted
 
     REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
@@ -704,18 +714,17 @@
        !< array to which output data is resorted to
 
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
+    
     found = .TRUE.
+    resorted = .FALSE.
+!
+!-- Set masking flag for topography for not resorted arrays
+    flag_nr = 0
 
     SELECT CASE ( TRIM( variable ) )
-
 
        CASE ( 'diss_xy', 'diss_xz', 'diss_yz' )
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss
           ELSE
              IF ( .NOT. ALLOCATED( diss_av ) ) THEN
@@ -723,24 +732,97 @@
                 diss_av = REAL( fill_value, KIND = wp )
              ENDIF
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_av(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_av
           ENDIF
-
           IF ( mode == 'xy' ) grid = 'zu'
 
        CASE ( 'kh_xy', 'kh_xz', 'kh_yz' )
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO k = nzb_do, nzt_do
-                      local_pf(i,j,k) = kh(k,j,i)
-                   ENDDO
-                ENDDO
+             to_be_resorted => kh
+          ELSE
+             IF ( .NOT. ALLOCATED( kh_av ) ) THEN
+                ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                kh_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => kh_av
+          ENDIF
+          IF ( mode == 'xy' ) grid = 'zu'
+
+       CASE ( 'km_xy', 'km_xz', 'km_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => km
+          ELSE
+             IF ( .NOT. ALLOCATED( km_av ) ) THEN
+                ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                km_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => km_av
+          ENDIF
+          IF ( mode == 'xy' ) grid = 'zu'
+
+       CASE DEFAULT
+          found = .FALSE.
+          grid  = 'none'
+
+    END SELECT
+
+    IF ( found .AND. .NOT. resorted )  THEN
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_do, nzt_do
+                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
+                                         REAL( fill_value, KIND = wp ),        &
+                                         BTEST( wall_flags_0(k,j,i), flag_nr ) ) 
              ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+ 
+ END SUBROUTINE tcm_data_output_2d
+
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Define 3D output variables.
+!------------------------------------------------------------------------------!
+ SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
+ 
+
+    USE averaging,                                                             &
+        ONLY:  diss_av, kh_av, km_av
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  variable   !< name of variable
+
+    INTEGER(iwp) ::  av        !< flag for (non-)average output
+    INTEGER(iwp) ::  flag_nr   !< number of masking flag
+    INTEGER(iwp) ::  i         !< loop index
+    INTEGER(iwp) ::  j         !< loop index
+    INTEGER(iwp) ::  k         !< loop index
+    INTEGER(iwp) ::  nzb_do    !< lower limit of the data output (usually 0)
+    INTEGER(iwp) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
+
+    LOGICAL ::  found     !< flag if output variable is found
+    LOGICAL ::  resorted  !< flag if output is already resorted
+
+    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
+
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< local
+       !< array to which output data is resorted to
+
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
+
+    found = .TRUE.
+    resorted = .FALSE.
+!
+!-- Set masking flag for topography for not resorted arrays
+    flag_nr = 0
+
+    SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'diss' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => diss
           ELSE
              IF ( .NOT. ALLOCATED( diss_av ) ) THEN
@@ -748,117 +830,10 @@
                 diss_av = REAL( fill_value, KIND = wp )
              ENDIF
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO k = nzb_do, nzt_do
-                      local_pf(i,j,k) = kh_av(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-          IF ( mode == 'xy' ) grid = 'zu'
-
-       CASE ( 'km_xy', 'km_xz', 'km_yz' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO k = nzb_do, nzt_do
-                      local_pf(i,j,k) = km(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             IF ( .NOT. ALLOCATED( diss_av ) ) THEN
-                ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                diss_av = REAL( fill_value, KIND = wp )
-             ENDIF
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO k = nzb_do, nzt_do
-                      local_pf(i,j,k) = km_av(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ENDIF
-
-          IF ( mode == 'xy' ) grid = 'zu'
-
-       CASE DEFAULT
-          found = .FALSE.
-          grid  = 'none'
-
-    END SELECT
- 
- END SUBROUTINE tcm_data_output_2d
-
- 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Define 3D output variables.
-!------------------------------------------------------------------------------!
- SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
- 
-
-    USE averaging,                                                             &
-        ONLY:  diss_av, kh_av, km_av
-
-    IMPLICIT NONE
-
-    CHARACTER (LEN=*) ::  variable   !< name of variable
-
-    INTEGER(iwp) ::  av     !< flag for (non-)average output
-    INTEGER(iwp) ::  i      !< loop index
-    INTEGER(iwp) ::  j      !< loop index
-    INTEGER(iwp) ::  k      !< loop index
-    INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
-    INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
-
-    LOGICAL ::  found   !< flag if output variable is found
-
-    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
-
-    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< local
-       !< array to which output data is resorted to
-
-
-    found = .TRUE.
-
-
-    SELECT CASE ( TRIM( variable ) )
-
-
-       CASE ( 'diss' )
-          IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
-          ELSE
-             IF ( .NOT. ALLOCATED( diss_av ) ) THEN
-                ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                diss_av = REAL( fill_value, KIND = wp )
-             ENDIF
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_av(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_av
           ENDIF
 
        CASE ( 'kh' )
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = kh(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => kh
           ELSE
              IF ( .NOT. ALLOCATED( kh_av ) ) THEN
@@ -866,22 +841,10 @@
                 kh_av = REAL( fill_value, KIND = wp )
              ENDIF
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = kh_av(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => kh_av
           ENDIF
 
        CASE ( 'km' )
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = km(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => km
           ELSE
              IF ( .NOT. ALLOCATED( km_av ) ) THEN
@@ -889,143 +852,65 @@
                 km_av = REAL( fill_value, KIND = wp )
              ENDIF
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = km_av(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => km_av
           ENDIF
 
        CASE ( 'dummy3' )                                                        !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = dummy3(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => dummy3
           ENDIF
 
        CASE ( 'diss1' )                                                         !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = dummy1(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => dummy1
           ENDIF
 
        CASE ( 'diss2' )                                                         !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = dummy2(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => dummy2
           ENDIF
 
        CASE ( 'diss_prod1' )                                                    !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_prod1(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_prod1
           ENDIF
 
        CASE ( 'diss_adve1' )                                                    !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_adve1(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_adve1
           ENDIF
 
        CASE ( 'diss_diff1' )                                                    !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_diff1(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_diff1
           ENDIF
 
        CASE ( 'diss_prod2' )                                                    !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_prod2(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_prod2
           ENDIF
 
        CASE ( 'diss_adve2' )                                                    !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_adve2(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_adve2
           ENDIF
 
        CASE ( 'diss_diff2' )                                                    !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_diff2(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_diff2
           ENDIF
 
        CASE ( 'diss_prod3' )                                                    !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_prod3(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_prod3
           ENDIF
 
        CASE ( 'diss_adve3' )                                                    !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_adve3(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_adve3
           ENDIF
 
        CASE ( 'diss_diff3' )                                                    !> @todo remove later
           IF ( av == 0 )  THEN
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = diss_diff3(k,j,i)
-                   ENDDO
-                ENDDO
-             ENDDO
+             to_be_resorted => diss_diff3
           ENDIF
           
@@ -1034,4 +919,19 @@
 
     END SELECT
+
+
+    IF ( found .AND. .NOT. resorted )  THEN
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_do, nzt_do
+                local_pf(i,j,k) = MERGE(                                 &
+                                   to_be_resorted(k,j,i),                &
+                                   REAL( fill_value, KIND = wp ),        &
+                                   BTEST( wall_flags_0(k,j,i), flag_nr ) )
+             ENDDO
+          ENDDO
+       ENDDO
+       resorted = .TRUE.
+    ENDIF
 
  END SUBROUTINE tcm_data_output_3d
@@ -1168,34 +1068,63 @@
 
        IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
-!
-!--       Transfer initial profiles to the arrays of the 3D model
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                e(:,j,i)  = e1d
-                kh(:,j,i) = kh1d
-                km(:,j,i) = km1d
+
+          IF ( .NOT. rans_tke_e ) THEN
+!
+!--          Transfer initial profiles to the arrays of the 3D model
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   e(:,j,i)  = e1d
+                   kh(:,j,i) = kh1d
+                   km(:,j,i) = km1d
+                ENDDO
              ENDDO
-          ENDDO
-
-          IF ( constant_diffusion )  THEN
-             e = 0.0_wp
-          ENDIF
-
-          IF ( rans_tke_e )  THEN
-             IF ( dissipation_1d == 'prognostic' )  THEN    !> @query Why must this be checked?
-                DO  i = nxlg, nxrg                          !>   Should 'diss' not always
-                   DO  j = nysg, nyng                       !>   be prognostic in case rans_tke_e?
-                      diss(:,j,i) = diss1d
-                   ENDDO
-                ENDDO
-             ELSE
+
+             IF ( constant_diffusion )  THEN
+                e = 0.0_wp
+             ENDIF
+
+          ELSE
+!
+!--          In case of TKE-e closure in RANS mode, do not use e, diss, and km
+!--          profiles from 1D model. Instead, initialize with constant profiles
+             IF ( constant_diffusion )  THEN
+                km = km_constant
+                kh = km / prandtl_number
+                e  = 0.0_wp
+             ELSEIF ( e_init > 0.0_wp )  THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
                       DO  k = nzb+1, nzt
-                         diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km1d(k)
+                         km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
                       ENDDO
                    ENDDO
                 ENDDO
+                km(nzb,:,:)   = km(nzb+1,:,:)
+                km(nzt+1,:,:) = km(nzt,:,:)
+                kh = km / prandtl_number
+                e  = e_init
+             ELSE
+                IF ( .NOT. ocean )  THEN
+                   kh   = 0.01_wp   ! there must exist an initial diffusion, because 
+                   km   = 0.01_wp   ! otherwise no TKE would be produced by the
+                                     ! production terms, as long as not yet
+                                     ! e = (u*/cm)**2 at k=nzb+1
+                ELSE
+                   kh   = 0.00001_wp
+                   km   = 0.00001_wp
+                ENDIF
+                e    = 0.0_wp
              ENDIF
+
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb+1, nzt
+                      diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+             diss(nzb,:,:) = diss(nzb+1,:,:)
+             diss(nzt+1,:,:) = diss(nzt,:,:)
+
           ENDIF
 
@@ -1381,6 +1310,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  bc_lr_cyc, bc_ns_cyc, f, kappa, message_string,                 &
-               wall_adjustment_factor
+        ONLY:  bc_lr_cyc, bc_ns_cyc, f, message_string, wall_adjustment_factor
 
     USE grid_variables,                                                        &
@@ -1973,5 +1901,5 @@
 
     USE surface_mod,                                                           &
-        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_usm_h 
+        ONLY :  surf_def_h, surf_lsm_h, surf_usm_h 
 
     IMPLICIT NONE
@@ -2556,24 +2484,4 @@
 
 !
-!--    Use special boundary condition in case of TKE-e closure
-       IF ( rans_tke_e )  THEN
-          surf_s = surf_def_h(0)%start_index(j,i)
-          surf_e = surf_def_h(0)%end_index(j,i)
-          DO  m = surf_s, surf_e
-             k = surf_def_h(0)%k(m)
-             e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
-          ENDDO
-
-          DO  l = 0, 3
-             surf_s = surf_def_v(l)%start_index(j,i)
-             surf_e = surf_def_v(l)%end_index(j,i)
-             DO  m = surf_s, surf_e
-                k = surf_def_v(l)%k(m)
-                e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
-             ENDDO
-          ENDDO
-       ENDIF
-
-!
 !--    Calculate tendencies for the next Runge-Kutta step
        IF ( timestep_scheme(1:5) == 'runge' )  THEN
@@ -2675,30 +2583,4 @@
 
 !
-!--    Use special boundary condition in case of TKE-e closure
-       surf_s = surf_def_h(0)%start_index(j,i)
-       surf_e = surf_def_h(0)%end_index(j,i)
-       DO  m = surf_s, surf_e
-          k = surf_def_h(0)%k(m)
-          diss_p(k,j,i) = surf_def_h(0)%us(m)**3 / kappa * ddzu(k)
-       ENDDO
-
-       DO  l = 0, 1
-          surf_s = surf_def_v(l)%start_index(j,i)
-          surf_e = surf_def_v(l)%end_index(j,i)
-          DO  m = surf_s, surf_e
-             k = surf_def_v(l)%k(m)
-             diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dy )
-          ENDDO
-       ENDDO
-
-       DO  l = 2, 3
-          surf_s = surf_def_v(l)%start_index(j,i)
-          surf_e = surf_def_v(l)%end_index(j,i)
-          DO  m = surf_s, surf_e
-             k = surf_def_v(l)%k(m)
-             diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dx )
-          ENDDO
-       ENDDO
-!
 !--    Calculate tendencies for the next Runge-Kutta step
        IF ( timestep_scheme(1:5) == 'runge' )  THEN
@@ -2715,13 +2597,4 @@
           ENDIF
        ENDIF
-!
-!--    Limit change of diss to be between -90% and +100%. Also, set an absolute
-!--    minimum value
-       DO  k = nzb, nzt+1
-          diss_p(k,j,i) = MIN( MAX( diss_p(k,j,i),         &
-                                    0.1_wp * diss(k,j,i),  &
-                                    0.0001_wp ),           &
-                               2.0_wp * diss(k,j,i) )
-       ENDDO
 
        IF ( intermediate_timestep_count == 1 )  dummy1(:,j,i) = diss_p(:,j,i)   !> @todo remove later
@@ -4777,5 +4650,4 @@
 !> Computation of the turbulent diffusion coefficients for momentum and heat 
 !> according to Prandtl-Kolmogorov.
-!> @todo consider non-default surfaces
 !------------------------------------------------------------------------------!
  SUBROUTINE tcm_diffusivities_default( var, var_reference )
@@ -4792,5 +4664,6 @@
 
     USE surface_mod,                                                           &
-        ONLY :  bc_h, surf_def_h, surf_def_v
+        ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
+                surf_usm_h, surf_usm_v
 
     IMPLICIT NONE
@@ -4956,4 +4829,5 @@
 !
 !--    Up- and downward facing surfaces
+!--    Default surfaces
        DO  n = 0, 1
           DO  m = 1, surf_def_h(n)%ns
@@ -4961,27 +4835,56 @@
              j = surf_def_h(n)%j(m)
              k = surf_def_h(n)%k(m)
-             km(k,j,i) = kappa * surf_def_h(n)%us(m) * dzu(k)
+             km(k,j,i) = kappa * surf_def_h(n)%us(m) * surf_def_h(n)%z_mo(m)
              kh(k,j,i) = 1.35_wp * km(k,j,i)
           ENDDO
        ENDDO
 !
-!--    North- and southward facing surfaces
-       DO  n = 0, 1
+!--    Upward facing surfaces
+!--    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)
+          km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m)
+          kh(k,j,i) = 1.35_wp * km(k,j,i)
+       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)
+          km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m)
+          kh(k,j,i) = 1.35_wp * km(k,j,i)
+       ENDDO
+
+!
+!--    North-, south-, west and eastward facing surfaces
+       DO  n = 0, 3
+!
+!--       Default surfaces
           DO  m = 1, surf_def_v(n)%ns
              i = surf_def_v(n)%i(m)
              j = surf_def_v(n)%j(m)
              k = surf_def_v(n)%k(m)
-             km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dy
+             km(k,j,i) = kappa * surf_def_v(n)%us(m) * surf_def_v(n)%z_mo(m)
              kh(k,j,i) = 1.35_wp * km(k,j,i)
           ENDDO
-       ENDDO
-!
-!--    West- and eastward facing surfaces
-       DO  n = 2, 3
-          DO  m = 1, surf_def_v(n)%ns
-             i = surf_def_v(n)%i(m)
-             j = surf_def_v(n)%j(m)
-             k = surf_def_v(n)%k(m)
-             km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dx
+!
+!--       Natural surfaces
+          DO  m = 1, surf_lsm_v(n)%ns
+             i = surf_lsm_v(n)%i(m)
+             j = surf_lsm_v(n)%j(m)
+             k = surf_lsm_v(n)%k(m)
+             km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m)
+             kh(k,j,i) = 1.35_wp * km(k,j,i)
+          ENDDO
+!
+!--       Urban surfaces
+          DO  m = 1, surf_usm_v(n)%ns
+             i = surf_usm_v(n)%i(m)
+             j = surf_usm_v(n)%j(m)
+             k = surf_usm_v(n)%k(m)
+             km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m)
              kh(k,j,i) = 1.35_wp * km(k,j,i)
           ENDDO
