Index: palm/trunk/SOURCE/model_1d_mod.f90
===================================================================
--- palm/trunk/SOURCE/model_1d_mod.f90	(revision 3049)
+++ palm/trunk/SOURCE/model_1d_mod.f90	(revision 3083)
@@ -25,4 +25,21 @@
 ! -----------------
 ! $Id$
+! Bugfixes:
+!   - preset te_diss and te_e to avoid runtime errors
+!   - implementation of buoyancy term to dissipation
+!     according to Sogachev et al. (2012)
+!   - where diss_p < 0 set diss_p = 0.1 diss
+!   - calculate progn eq(diss) starting from nzb+1
+! Changes:
+!   - add sig_e to TKE equation
+!   - adjust prognostic equation of diss
+!   - set model constants according to Koblitz (2013)
+!   - renamed c_m to c_0
+!   - rename l_black into l1d_init
+!   - calculate l_grid within init_1d_model and save it as l1d_init
+!   - calculate l1d according to DE85 if dissipation is a prognostic value
+!   - made annotations doxygen-readable
+! 
+! 3049 2018-05-29 13:52:36Z Giersch
 ! Error messages revised
 ! 
@@ -164,5 +181,6 @@
     USE kinds
 
-    USE pegrid
+    USE pegrid,                                                                &
+        ONLY:  myid
         
 
@@ -175,10 +193,10 @@
     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) ::  c_mu = 0.09_wp                !< model constant
+    REAL(wp) ::  alpha_buoyancy                !< model constant according to Koblitz (2013)
+    REAL(wp) ::  c_0 = 0.03_wp**0.25_wp        !< model constant according to Koblitz (2013)
+    REAL(wp) ::  c_1 = 1.52_wp                 !< model constant according to Koblitz (2013)
+    REAL(wp) ::  c_2 = 1.83_wp                 !< model constant according to Koblitz (2013)
+    REAL(wp) ::  c_3                           !< model constant
+    REAL(wp) ::  c_mu                          !< model constant
     REAL(wp) ::  damp_level_1d = -1.0_wp       !< namelist parameter
     REAL(wp) ::  dt_1d = 60.0_wp               !< dynamic timestep (1d-model)
@@ -187,7 +205,9 @@
     REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter
     REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter
+    REAL(wp) ::  lambda                        !< maximum mixing length
     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) ::  sig_diss = 2.95_wp            !< model constant according to Koblitz (2013)
+    REAL(wp) ::  sig_e = 2.95_wp               !< model constant according to Koblitz (2013)
     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)
@@ -198,5 +218,4 @@
     REAL(wp) ::  z01d                          !< roughness length for momentum (1d-model)
     REAL(wp) ::  z0h1d                         !< roughness length for scalars (1d-model)
-
 
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  diss1d   !< tke dissipation rate (1d-model)
@@ -279,6 +298,4 @@
 
     INTEGER(iwp) ::  k  !< loop index
-    
-    REAL(wp) ::  lambda !< maximum mixing length
 
 !
@@ -324,11 +341,11 @@
 !
 !--       Use the same mixing length as in 3D model (LES-mode)
-          !@todo: rename (delete?) this option
-          ! As the mixing length is different between RANS and LES mode, it
-          ! must be distinguished here between these modes. For RANS mode,
-          ! the mixing length is calculated accoding to Blackadar, which is
-          ! the other option at this point.
-          ! Maybe delete this option entirely (not appropriate in LES case)
-          ! 2018-03-20, gronemeier
+          !> @todo rename (delete?) this option
+          !>  As the mixing length is different between RANS and LES mode, it
+          !>  must be distinguished here between these modes. For RANS mode,
+          !>  the mixing length is calculated accoding to Blackadar, which is
+          !>  the other option at this point.
+          !>  Maybe delete this option entirely (not appropriate in LES case)
+          !>  2018-03-20, gronemeier
           DO  k = nzb+1, nzt
              l1d_init(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
@@ -359,5 +376,5 @@
        us1d = 0.1_wp   ! without initial friction the flow would not change
     ELSE
-       diss1d(nzb+1) = 1.0_wp
+       diss1d(nzb+1) = 0.001_wp
        e1d(nzb+1)  = 1.0_wp
        km1d(nzb+1) = 1.0_wp
@@ -372,7 +389,8 @@
 
 !
-!-- Tendencies must be preset in order to avoid runtime errors within the
-!-- first Runge-Kutta step
+!-- Tendencies must be preset in order to avoid runtime errors
+    te_diss  = 0.0_wp
     te_dissm = 0.0_wp
+    te_e  = 0.0_wp
     te_em = 0.0_wp
     te_um = 0.0_wp
@@ -381,5 +399,6 @@
 !
 !-- Set model constant
-    IF ( dissipation_1d == 'as_in_3d_model' )  c_m = 0.1_wp
+    IF ( dissipation_1d == 'as_in_3d_model' )  c_0 = 0.1_wp
+    c_mu = c_0**4
 
 !
@@ -423,5 +442,5 @@
 !-- Determine the time step at the start of a 1D-simulation and
 !-- determine and printout quantities used for run control
-    dt_1d = 1.0_wp
+    dt_1d = 0.01_wp
     CALL run_control_1d
 
@@ -483,5 +502,5 @@
 !--             dissipation rate
                 IF ( dissipation_1d == 'detering' )  THEN
-                   diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
+                   diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
                 ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
                    diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) &
@@ -497,5 +516,5 @@
                                      kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
                                    - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
-                                                 ) * ddzw(k)                   &
+                                                 ) * ddzw(k) / sig_e           &
                                    - diss1d(k)
 
@@ -503,17 +522,24 @@
 !
 !--                dissipation rate
-                   te_diss(k) = km1d(k) *                                      &
+                   IF ( rif1d(k) >= 0.0_wp )  THEN
+                      alpha_buoyancy = 1.0_wp - l1d(k) / lambda
+                   ELSE
+                      alpha_buoyancy = 1.0_wp - ( 1.0_wp + ( c_2 - 1.0_wp )    &
+                                                         / ( c_2 - c_1    ) )  &
+                                              * l1d(k) / lambda
+                   ENDIF
+                   c_3 = ( c_1 - c_2 ) * alpha_buoyancy
+                   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))                             &
+                                  ) * ( c_1 + (c_2 - c_1) * l1d(k) / lambda )  &
                                   - g / pt_0 * kh1d(k) * flux * c_3            &
-                                    * diss1d(k) / ( e1d(k) + 1.0E-20_wp )      &
-                                  + ( kmzp * ( diss1d(k+1) - diss1d(k) )       &
+                                  - c_2 * diss1d(k)                            &
+                                ) * 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 )
+                                  ) * ddzw(k) / sig_diss
 
                 ENDIF
@@ -546,5 +572,5 @@
 !--          dissipation rate
              IF ( dissipation_1d == 'detering' )  THEN
-                diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
+                diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
              ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
                 diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) )  &
@@ -565,4 +591,8 @@
 !--          TKE
              IF ( .NOT. dissipation_1d == 'prognostic' )  THEN
+                !> @query why integrate over 2dz
+                !>   Why is it allowed to integrate over two delta-z for e
+                !>   while for u and v it is not?
+                !>   2018-04-23, gronemeier
                 te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
                                     + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
@@ -572,5 +602,5 @@
                                      kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
                                    - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
-                                                 ) * ddzw(k)                   &
+                                                 ) * ddzw(k) / sig_e           &
                                    - diss1d(k)
              ENDIF
@@ -595,10 +625,4 @@
              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
@@ -606,4 +630,12 @@
 !--          value is reduced to 10 percent of its old value.
              WHERE ( e1d_p < 0.0_wp )  e1d_p = 0.1_wp * e1d
+
+             IF ( dissipation_1d == 'prognostic' )  THEN
+                DO  k = nzb+1, nzt
+                   diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
+                                                       tsc(3) * te_dissm(k) )
+                ENDDO
+                WHERE ( diss1d_p < 0.0_wp )  diss1d_p = 0.1_wp * diss1d
+             ENDIF
           ENDIF
 
@@ -643,5 +675,6 @@
                    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)
+                         te_dissm(k) = -9.5625_wp * te_diss(k)  &
+                                     +  5.3125_wp * te_dissm(k)
                       ENDDO
                    ENDIF
@@ -650,5 +683,4 @@
              ENDIF
           ENDIF
-
 
 !
@@ -701,5 +733,8 @@
 
              ENDIF    ! constant_flux_layer
-
+             !> @todo combine if clauses
+             !>   The previous and following if clauses can be combined into a
+             !>   single clause
+             !>   2018-04-23, gronemeier
 !
 !--          Compute the Richardson-flux numbers,
@@ -789,8 +824,8 @@
 !--             compatibility with the 3D model.
                 IF ( ibc_e_b == 2 )  THEN
-                   e1d(nzb+1) = ( us1d / c_m )**2
+                   e1d(nzb+1) = ( us1d / c_0 )**2
                 ENDIF
                 IF ( dissipation_1d == 'prognostic' )  THEN
-                   e1d(nzb+1) = us1d**2 / SQRT( c_mu )
+                   e1d(nzb+1) = ( us1d / c_0 )**2
                    diss1d(nzb+1) = us1d**3 / ( kappa * zu(nzb+1) )
                    diss1d(nzb) = diss1d(nzb+1)
@@ -829,26 +864,33 @@
 !--          in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
 !--          too large.
-             IF ( mixing_length_1d == 'blackadar' )  THEN
+             IF ( dissipation_1d /= 'prognostic' )  THEN
+                IF ( mixing_length_1d == 'blackadar' )  THEN
+                   DO  k = nzb+1, nzt
+                      IF ( rif1d(k) >= 0.0_wp )  THEN
+                         l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
+                         l1d_diss(k) = l1d(k)
+                      ELSE
+                         l1d(k) = l1d_init(k)
+                         l1d_diss(k) = l1d_init(k) *                           &
+                                       SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
+                      ENDIF
+                   ENDDO
+                ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
+                   DO  k = nzb+1, nzt
+                      dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
+                      IF ( dpt_dz > 0.0_wp )  THEN
+                         l_stable = 0.76_wp * SQRT( e1d(k) )                   &
+                                  / SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
+                      ELSE
+                         l_stable = l1d_init(k)
+                      ENDIF
+                      l1d(k) = MIN( l1d_init(k), l_stable )
+                      l1d_diss(k) = l1d(k)
+                   ENDDO
+                ENDIF
+             ELSE
                 DO  k = nzb+1, nzt
-                   IF ( rif1d(k) >= 0.0_wp )  THEN
-                      l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
-                      l1d_diss(k) = l1d(k)
-                   ELSE
-                      l1d(k) = l1d_init(k)
-                      l1d_diss(k) = l1d_init(k) *                              &
-                                    SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
-                   ENDIF
-                ENDDO
-             ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
-                DO  k = nzb+1, nzt
-                   dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
-                   IF ( dpt_dz > 0.0_wp )  THEN
-                      l_stable = 0.76_wp * SQRT( e1d(k) ) /                    &
-                                     SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
-                   ELSE
-                      l_stable = l1d_init(k)
-                   ENDIF
-                   l1d(k) = MIN( l1d_init(k), l_stable )
-                   l1d_diss(k) = l1d(k)
+                   l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) )                   &
+                          / ( diss1d(k) + 1.0E-30_wp )
                 ENDDO
              ENDIF
@@ -867,4 +909,5 @@
                 ENDIF
              ENDIF
+
              IF ( dissipation_1d == 'prognostic' )  THEN
                 DO  k = nzb_diff, nzt
@@ -873,5 +916,5 @@
              ELSE
                 DO  k = nzb_diff, nzt
-                   km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
+                   km1d(k) = c_0 * SQRT( e1d(k) ) * l1d(k)
                 ENDDO
              ENDIF
@@ -1018,11 +1061,16 @@
 
     REAL(wp) ::  dt_diff  !< time step accorind to diffusion criterion
+    REAL(wp) ::  dt_old   !< previous time step
     REAL(wp) ::  fac      !< factor of criterion
     REAL(wp) ::  value    !< auxiliary variable
 
 !
+!-- Save previous time step
+    dt_old = dt_1d
+
+!
 !-- Compute the currently feasible time step according to the diffusion
 !-- criterion. At nzb+1 the half grid length is used.
-    fac = 0.125  !0.35_wp                                                       !### changed from 0.35
+    fac = 0.125
     dt_diff = dt_max_1d
     DO  k = nzb+2, nzt
@@ -1034,6 +1082,10 @@
 
 !
+!-- Limit the new time step to a maximum of 10 times the previous time step
+    dt_1d = MIN( dt_old * 10.0_wp, dt_1d )
+
+!
 !-- Set flag when the time step becomes too small
-    IF ( dt_1d < ( 0.00001_wp * dt_max_1d ) )  THEN
+    IF ( dt_1d < ( 1.0E-15_wp * dt_max_1d ) )  THEN
        stop_dt_1d = .TRUE.
 
