Index: palm/trunk/SOURCE/model_1d_mod.f90
===================================================================
--- palm/trunk/SOURCE/model_1d_mod.f90	(revision 2339)
+++ palm/trunk/SOURCE/model_1d_mod.f90	(revision 2696)
@@ -1,5 +1,5 @@
 !> @file model_1d_mod.f90
 !------------------------------------------------------------------------------!
-! This file is part of PALM.
+! This file is part of the PALM model system.
 !
 ! PALM is free software: you can redistribute it and/or modify it under the
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! implement TKE-e closure
+! modification of dissipation production according to Detering and Etling
+! reduced factor for timestep criterion to 0.125 and first dt to 1s (TG)
+! 
+! 2339 2017-08-07 13:55:26Z gronemeier
 ! corrected timestamp in header
 ! 
@@ -131,5 +136,6 @@
                intermediate_timestep_count_max, kappa, km_constant,            &
                message_string, mixing_length_1d, prandtl_number,               &
-               roughness_length, run_description_header, simulated_time_chr, timestep_scheme, tsc, z0h_factor
+               roughness_length, run_description_header, simulated_time_chr,   &
+               timestep_scheme, tsc, z0h_factor
 
     USE indices,                                                               &
@@ -143,38 +149,47 @@
     IMPLICIT NONE
 
-
     INTEGER(iwp) ::  current_timestep_number_1d = 0  !< current timestep number (1d-model)
-    INTEGER(iwp) ::  damp_level_ind_1d               !< lower grid index of damping layer (1d-model) !#
+    INTEGER(iwp) ::  damp_level_ind_1d               !< lower grid index of damping layer (1d-model)
 
     LOGICAL ::  run_control_header_1d = .FALSE.  !< flag for output of run control header (1d-model)
     LOGICAL ::  stop_dt_1d = .FALSE.             !< termination flag, used in case of too small timestep (1d-model)
 
+    REAL(wp) ::  c_1 = 1.44_wp                 !< model constant
+    REAL(wp) ::  c_2 = 1.92_wp                 !< model constant
+    REAL(wp) ::  c_3 = 1.44_wp                 !< model constant
+    REAL(wp) ::  c_h = 0.0015_wp               !< model constant according to Detering and Etling (1985)
     REAL(wp) ::  c_m = 0.4_wp                  !< model constant, 0.4 according to Detering and Etling (1985)
-    REAL(wp) ::  damp_level_1d = -1.0_wp       !< namelist parameter    !#
+    REAL(wp) ::  c_mu = 0.09_wp                !< model constant
+    REAL(wp) ::  damp_level_1d = -1.0_wp       !< namelist parameter
     REAL(wp) ::  dt_1d = 60.0_wp               !< dynamic timestep (1d-model)
     REAL(wp) ::  dt_max_1d = 300.0_wp          !< timestep limit (1d-model)
-    REAL(wp) ::  dt_pr_1d = 9999999.9_wp       !< namelist parameter     !#
-    REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter     !#
-    REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter     !#
+    REAL(wp) ::  dt_pr_1d = 9999999.9_wp       !< namelist parameter
+    REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter
+    REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter
     REAL(wp) ::  qs1d                          !< characteristic humidity scale (1d-model)
     REAL(wp) ::  simulated_time_1d = 0.0_wp    !< updated simulated time (1d-model)
+    REAL(wp) ::  sig_diss = 1.3_wp             !< model constant
     REAL(wp) ::  time_pr_1d = 0.0_wp           !< updated simulated time for profile output (1d-model)
     REAL(wp) ::  time_run_control_1d = 0.0_wp  !< updated simulated time for run-control output (1d-model)
     REAL(wp) ::  ts1d                          !< characteristic temperature scale (1d-model)
-    REAL(wp) ::  us1d                          !< friction velocity (1d-model)                 !#
-    REAL(wp) ::  usws1d                        !< u-component of the momentum flux (1d-model)  !#
-    REAL(wp) ::  vsws1d                        !< v-component of the momentum flux (1d-model)  !#
+    REAL(wp) ::  us1d                          !< friction velocity (1d-model)
+    REAL(wp) ::  usws1d                        !< u-component of the momentum flux (1d-model)
+    REAL(wp) ::  vsws1d                        !< v-component of the momentum flux (1d-model)
     REAL(wp) ::  z01d                          !< roughness length for momentum (1d-model)
     REAL(wp) ::  z0h1d                         !< roughness length for scalars (1d-model)
 
 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  e1d      !< tke (1d-model)                         !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  diss1d   !< tke dissipation rate (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  diss1d_p !< prognostic value of tke dissipation rate (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  e1d      !< tke (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  e1d_p    !< prognostic value of tke (1d-model)
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  kh1d     !< turbulent diffusion coefficient for heat (1d-model)    !#
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  km1d     !< turbulent diffusion coefficient for momentum (1d-model)!#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  kh1d     !< turbulent diffusion coefficient for heat (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  km1d     !< turbulent diffusion coefficient for momentum (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black  !< mixing length Blackadar (1d-model)
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d      !< mixing length for turbulent diffusion coefficients (1d-model) !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d      !< mixing length for turbulent diffusion coefficients (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d_diss !< mixing length for dissipation (1d-model)
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rif1d    !< Richardson flux number (1d-model)      !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rif1d    !< Richardson flux number (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_diss  !< tendency of diss (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_dissm !< weighted tendency of diss for previous sub-timestep (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_e     !< tendency of e (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_em    !< weighted tendency of e for previous sub-timestep (1d-model)
@@ -183,7 +198,7 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_v     !< tendency of v (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_vm    !< weighted tendency of v for previous sub-timestep (1d-model)
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u1d      !< u-velocity component (1d-model)        !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u1d      !< u-velocity component (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  u1d_p    !< prognostic value of u-velocity component (1d-model)
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v1d      !< v-velocity component (1d-model)        !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v1d      !< v-velocity component (1d-model)
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  v1d_p    !< prognostic value of v-velocity component (1d-model)
 
@@ -227,7 +242,7 @@
 !
 !-- Public variables
-    PUBLIC  damp_level_1d, damp_level_ind_1d, dt_pr_1d, dt_run_control_1d,     &
-            e1d, end_time_1d, kh1d, km1d, l1d, rif1d, u1d, us1d, usws1d, v1d,  &
-            vsws1d
+    PUBLIC  damp_level_1d, damp_level_ind_1d, diss1d, dt_pr_1d,                &
+            dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, rif1d, u1d,  &
+            us1d, usws1d, v1d, vsws1d
 
 
@@ -238,15 +253,17 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=9) ::  time_to_string  !<
+    CHARACTER (LEN=9) ::  time_to_string  !< function to transform time from real to character string
+
+    INTEGER(iwp) ::  k  !< loop index
     
-    INTEGER(iwp) ::  k  !< 
-    
-    REAL(wp) ::  lambda !<
+    REAL(wp) ::  lambda !< maximum mixing length
 
 !
 !-- Allocate required 1D-arrays
-    ALLOCATE( e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1),               &
+    ALLOCATE( diss1d(nzb:nzt+1), diss1d_p(nzb:nzt+1),                          &
+              e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1),               &
               km1d(nzb:nzt+1), l_black(nzb:nzt+1), l1d(nzb:nzt+1),             &
-              l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_e(nzb:nzt+1),          &
+              l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_diss(nzb:nzt+1),       &
+              te_dissm(nzb:nzt+1), te_e(nzb:nzt+1),                            &
               te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1),             &
               te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1),               &
@@ -259,4 +276,5 @@
        kh1d = km_constant / prandtl_number
     ELSE
+       diss1d = 0.0_wp; diss1d_p = 0.0_wp
        e1d = 0.0_wp; e1d_p = 0.0_wp
        kh1d = 0.0_wp; km1d = 0.0_wp
@@ -309,4 +327,5 @@
        us1d = 0.1_wp   ! without initial friction the flow would not change
     ELSE
+       diss1d(nzb+1) = 1.0_wp
        e1d(nzb+1)  = 1.0_wp
        km1d(nzb+1) = 1.0_wp
@@ -323,4 +342,5 @@
 !-- Tendencies must be preset in order to avoid runtime errors within the
 !-- first Runge-Kutta step
+    te_dissm = 0.0_wp
     te_em = 0.0_wp
     te_um = 0.0_wp
@@ -354,11 +374,10 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=9) ::  time_to_string  !<
-    
+    CHARACTER (LEN=9) ::  time_to_string  !< function to transform time from real to character string
+
     INTEGER(iwp) ::  k  !< loop index
-    
+
     REAL(wp) ::  a            !< auxiliary variable
     REAL(wp) ::  b            !< auxiliary variable
-    REAL(wp) ::  dissipation  !< dissipation of TKE
     REAL(wp) ::  dpt_dz       !< vertical temperature gradient
     REAL(wp) ::  flux         !< vertical temperature gradient
@@ -372,5 +391,5 @@
 !-- Determine the time step at the start of a 1D-simulation and
 !-- determine and printout quantities used for run control
-    dt_1d = 10.0_wp
+    dt_1d = 1.0_wp
     CALL run_control_1d
 
@@ -392,6 +411,6 @@
 
 !
-!--       Compute all tendency terms. If a Prandtl-layer is simulated, k starts 
-!--       at nzb+2.
+!--       Compute all tendency terms. If a constant-flux layer is simulated, 
+!--       k starts at nzb+2.
           DO  k = nzb_diff, nzt
 
@@ -414,5 +433,5 @@
              DO  k = nzb_diff, nzt
 !
-!--             TKE
+!--             TKE and dissipation rate
                 kmzm = 0.5_wp * ( km1d(k-1) + km1d(k) )
                 kmzp = 0.5_wp * ( km1d(k) + km1d(k+1) )
@@ -428,14 +447,15 @@
                 ENDIF
 
+!
+!--             Calculate dissipation rate if no prognostic equation is used for
+!--             dissipation rate
                 IF ( dissipation_1d == 'detering' )  THEN
-!
-!--                According to Detering, c_e=c_m^3
-                   dissipation = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
+                   diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
                 ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
-                   dissipation = ( 0.19_wp                                     &
-                                   + 0.74_wp * l1d_diss(k) /  l_grid(k)        &
-                                 ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
-                ENDIF
-
+                   diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) /  l_grid(k)  &
+                               ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
+                ENDIF
+!
+!--             TKE
                 te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
                                     + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
@@ -446,10 +466,30 @@
                                    - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
                                                  ) * ddzw(k)                   &
-                                   - dissipation
+                                   - diss1d(k)
+
+                IF ( dissipation_1d == 'prognostic' )  THEN
+!
+!--                dissipation rate
+                   te_diss(k) = km1d(k) *                                      &
+                                  ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2  &
+                                  + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2  &
+                                  ) * c_1 * c_mu**0.75 / c_h * f / us1d        &
+                                    * SQRT(e1d(k))                             &
+                                  - g / pt_0 * kh1d(k) * flux * c_3            &
+                                    * diss1d(k) / ( e1d(k) + 1.0E-20_wp )      &
+                                  + ( kmzp * ( diss1d(k+1) - diss1d(k) )       &
+                                           * ddzu(k+1)                         &
+                                    - kmzm * ( diss1d(k) - diss1d(k-1) )       &
+                                           * ddzu(k)                           &
+                                    ) * ddzw(k) / sig_diss                     &
+                                  - c_2 * diss1d(k)**2 / ( e1d(k) + 1.0E-20_wp )
+
+                ENDIF
+
              ENDDO
           ENDIF
 
 !
-!--       Tendency terms at the top of the Prandtl-layer.
+!--       Tendency terms at the top of the constant-flux layer.
 !--       Finite differences of the momentum fluxes are computed using half the 
 !--       normal grid length (2.0*ddzw(k)) for the sake of enhanced accuracy 
@@ -470,11 +510,12 @@
              ENDIF
 
+!
+!--          Calculate dissipation rate if no prognostic equation is used for
+!--          dissipation rate
              IF ( dissipation_1d == 'detering' )  THEN
-!
-!--             According to Detering, c_e=0.064
-                dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
+                diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
              ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
-                dissipation = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l_grid(k) )  &
-                              * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
+                diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l_grid(k) )    &
+                            * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
              ENDIF
 
@@ -491,13 +532,16 @@
 !
 !--          TKE
-             te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2   &
-                                 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2   &
-                                 )                                             &
-                                 - g / pt_0 * kh1d(k) * flux                   &
-                                 +           (                                 &
-                                  kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)     &
-                                - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)       &
-                                              ) * ddzw(k)                      &
-                                - dissipation
+             IF ( .NOT. dissipation_1d == 'prognostic' )  THEN
+                te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
+                                    + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
+                                    )                                          &
+                                    - g / pt_0 * kh1d(k) * flux                &
+                                    +           (                              &
+                                     kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
+                                   - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
+                                                 ) * ddzw(k)                   &
+                                   - diss1d(k)
+             ENDIF
+
           ENDIF
 
@@ -513,10 +557,16 @@
           ENDDO
           IF ( .NOT. constant_diffusion )  THEN
+
              DO  k = nzb+1, nzt
-
                 e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + &
                                               tsc(3) * te_em(k) )
-
              ENDDO
+
+             IF ( dissipation_1d == 'prognostic' )  THEN
+                DO  k = nzb_diff, nzt
+                   diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
+                                                 tsc(3) * te_dissm(k) )
+                ENDDO
+             ENDIF
 !
 !--          Eliminate negative TKE values, which can result from the
@@ -540,4 +590,9 @@
                       te_em(k) = te_e(k)
                    ENDDO
+                   IF ( dissipation_1d == 'prognostic' )  THEN
+                      DO k = nzb+1, nzt
+                         te_dissm(k) = te_diss(k)
+                      ENDDO
+                   ENDIF
                 ENDIF
 
@@ -554,4 +609,9 @@
                       te_em(k) = -9.5625_wp * te_e(k) + 5.3125_wp * te_em(k)
                    ENDDO
+                   IF ( dissipation_1d == 'prognostic' )  THEN
+                      DO k = nzb+1, nzt
+                         te_dissm(k) = -9.5625_wp * te_diss(k) + 5.3125_wp * te_dissm(k)
+                      ENDDO
+                   ENDIF
                 ENDIF
 
@@ -562,8 +622,8 @@
 !
 !--       Boundary conditions for the prognostic variables.
-!--       At the top boundary (nzt+1) u,v and e keep their initial values 
-!--       (ug(nzt+1), vg(nzt+1), 0).
+!--       At the top boundary (nzt+1) u, v, e, and diss keep their initial 
+!--       values (ug(nzt+1), vg(nzt+1), 0, 0).
 !--       At the bottom boundary, Dirichlet condition is used for u and v (0)
-!--       and Neumann condition for e (e(nzb)=e(nzb+1)).
+!--       and Neumann condition for e and diss (e(nzb)=e(nzb+1)).
           u1d_p(nzb) = 0.0_wp
           v1d_p(nzb) = 0.0_wp
@@ -575,4 +635,7 @@
           IF ( .NOT. constant_diffusion )  THEN
              e1d  = e1d_p
+             IF ( dissipation_1d == 'prognostic' )  THEN
+                diss1d = diss1d_p
+             ENDIF
           ENDIF
 
@@ -582,5 +645,5 @@
 
 !
-!--          First compute the vertical fluxes in the Prandtl-layer
+!--          First compute the vertical fluxes in the constant-flux layer
              IF ( constant_flux_layer )  THEN
 !
@@ -609,7 +672,7 @@
 !
 !--          Compute the Richardson-flux numbers,
-!--          first at the top of the Prandtl-layer using u* of the previous
-!--          time step (+1E-30, if u* = 0), then in the remaining area. There
-!--          the rif-numbers of the previous time step are used.
+!--          first at the top of the constant-flux layer using u* of the 
+!--          previous time step (+1E-30, if u* = 0), then in the remaining area.
+!--          There the rif-numbers of the previous time step are used.
 
              IF ( constant_flux_layer )  THEN
@@ -689,10 +752,15 @@
 
 !
-!--             Boundary condition for the turbulent kinetic energy at the top
-!--             of the Prandtl-layer. c_m = 0.4 according to Detering.
+!--             Boundary condition for the turbulent kinetic energy and 
+!--             dissipation rate at the top of the constant-flux layer.
 !--             Additional Neumann condition de/dz = 0 at nzb is set to ensure
 !--             compatibility with the 3D model.
                 IF ( ibc_e_b == 2 )  THEN
                    e1d(nzb+1) = ( us1d / c_m )**2
+                ENDIF
+                IF ( dissipation_1d == 'prognostic' )  THEN
+                   e1d(nzb+1) = us1d**2 / SQRT( c_mu )
+                   diss1d(nzb+1) = us1d**3 / ( kappa * zu(nzb+1) )
+                   diss1d(nzb) = diss1d(nzb+1)
                 ENDIF
                 e1d(nzb) = e1d(nzb+1)
@@ -767,7 +835,13 @@
                 ENDIF
              ENDIF
-             DO  k = nzb_diff, nzt
-                km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
-             ENDDO
+             IF ( dissipation_1d == 'prognostic' )  THEN
+                DO  k = nzb_diff, nzt
+                   km1d(k) = c_mu * e1d(k)**2 / ( diss1d(k) + 1.0E-30_wp )
+                ENDDO
+             ELSE
+                DO  k = nzb_diff, nzt
+                   km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
+                ENDDO
+             ENDIF
 
 !
@@ -918,5 +992,5 @@
 !-- Compute the currently feasible time step according to the diffusion
 !-- criterion. At nzb+1 the half grid length is used.
-    fac = 0.35_wp
+    fac = 0.125  !0.35_wp                                                       !### changed from 0.35
     dt_diff = dt_max_1d
     DO  k = nzb+2, nzt
@@ -977,5 +1051,5 @@
        DO  k = nzt+1, nzb, -1
           WRITE ( 17, 103)  k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &
-                            rif1d(k), km1d(k), kh1d(k), l1d(k), zu(k), k
+                            rif1d(k), km1d(k), kh1d(k), l1d(k), diss1d(k)
        ENDDO
        WRITE ( 17, 101 )
@@ -993,9 +1067,8 @@
 100 FORMAT ('# ',A/'#',10('-')/'# 1d-model profiles')
 104 FORMAT (//'# Time: ',A)
-101 FORMAT ('#',79('-'))
-102 FORMAT ('#  k     zu      u      v     pt      e    rif    Km    Kh     ', &
-            'l      zu      k')
-103 FORMAT (1X,I4,1X,F7.1,1X,F6.2,1X,F6.2,1X,F6.2,1X,F6.2,1X,F5.2,1X,F5.2, &
-            1X,F5.2,1X,F6.2,1X,F7.1,2X,I4)
+101 FORMAT ('#',111('-'))
+102 FORMAT ('#  k     zu      u          v          pt         e          ',   &
+            'rif        Km         Kh         l          diss   ')
+103 FORMAT (1X,I4,1X,F7.1,9(1X,E10.3))
 
 
