Index: palm/trunk/SOURCE/init_1d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_1d_model.f90	(revision 2334)
+++ palm/trunk/SOURCE/init_1d_model.f90	(revision 2337)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! revised calculation of mixing length
+! removed rounding of time step
+! corrected calculation of virtual potential temperature
+! 
+! 2334 2017-08-04 11:57:04Z gronemeier
 ! set c_m = 0.4 according to Detering and Etling (1985)
 ! 
@@ -120,5 +125,5 @@
     
     USE model_1d,                                                              &
-        ONLY:  e1d, e1d_p, kh1d, km1d, l1d, l_black, qs1d, rif1d,              &
+        ONLY:  e1d, e1d_p, kh1d, km1d, l1d, l1d_diss, l_black, qs1d, rif1d,    &
                simulated_time_1d, te_e, te_em, te_u, te_um, te_v, te_vm, ts1d, &
                u1d, u1d_p, us1d, usws1d, v1d, v1d_p, vsws1d, z01d, z0h1d
@@ -141,5 +146,5 @@
     ALLOCATE( 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),                              &
+              l_black(nzb:nzt+1), l1d(nzb:nzt+1), l1d_diss(nzb:nzt+1),         &
               rif1d(nzb:nzt+1),   te_e(nzb:nzt+1),                             &
               te_em(nzb:nzt+1),  te_u(nzb:nzt+1),    te_um(nzb:nzt+1),         &
@@ -183,9 +188,10 @@
        ENDIF
     ENDIF
-    l1d   = l_black
-    u1d   = u_init
-    u1d_p = u_init
-    v1d   = v_init
-    v1d_p = v_init
+    l1d      = l_black
+    l1d_diss = l_black
+    u1d      = u_init
+    u1d_p    = u_init
+    v1d      = v_init
+    v1d_p    = v_init
 
 !
@@ -226,5 +232,5 @@
 
 !
-!-- Integrate the 1D-model equations using the leap-frog scheme
+!-- Integrate the 1D-model equations using the Runge-Kutta scheme
     CALL time_integration_1d
 
@@ -261,5 +267,5 @@
         ONLY:  current_timestep_number_1d, damp_level_ind_1d, dt_1d,           &
                dt_pr_1d, dt_run_control_1d, e1d, e1d_p, end_time_1d,           &
-               kh1d, km1d, l1d, l_black, qs1d, rif1d, simulated_time_1d,       &
+               kh1d, km1d, l1d, l1d_diss, l_black, qs1d, rif1d, simulated_time_1d, &
                stop_dt_1d, te_e, te_em, te_u, te_um, te_v, te_vm, time_pr_1d,  &
                ts1d, time_run_control_1d, u1d, u1d_p, us1d, usws1d, v1d,       &
@@ -340,15 +346,17 @@
                    pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
                    flux = ( ( pt_init(k+1) - pt_init(k-1) ) +                  &
-                            0.61_wp * pt_init(k) *                             &
-                            ( q_init(k+1) - q_init(k-1) ) ) * dd2zu(k)
+                            0.61_wp * ( pt_init(k+1) * q_init(k+1) -           &
+                                        pt_init(k-1) * q_init(k-1)   )         &
+                          ) * dd2zu(k)
                 ENDIF
 
                 IF ( dissipation_1d == 'detering' )  THEN
 !
-!--                According to Detering, c_e=0.064
-                   dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d(k)
+!--                According to Detering, c_e=c_m^3
+                   dissipation = 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(k) / l_grid(k) )    &
-                                 * e1d(k) * SQRT( e1d(k) ) / l1d(k)
+                   dissipation = ( 0.19_wp                                     &
+                                   + 0.74_wp * l1d_diss(k) /  l_grid(k)        &
+                                 ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
                 ENDIF
 
@@ -380,5 +388,6 @@
                 pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
                 flux = ( ( pt_init(k+1) - pt_init(k-1) ) +                     &
-                         0.61_wp * pt_init(k) * ( q_init(k+1) - q_init(k-1) )  &
+                         0.61_wp * ( pt_init(k+1) * q_init(k+1) -              &
+                                     pt_init(k-1) * q_init(k-1)   )            &
                        ) * dd2zu(k)
              ENDIF
@@ -387,8 +396,8 @@
 !
 !--             According to Detering, c_e=0.064
-                dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d(k)
+                dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
              ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
-                dissipation = ( 0.19_wp + 0.74_wp * l1d(k) / l_grid(k) )       &
-                              * e1d(k) * SQRT( e1d(k) ) / l1d(k)
+                dissipation = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l_grid(k) )  &
+                              * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
              ENDIF
 
@@ -513,19 +522,8 @@
                    b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /                 &
                        zu(nzb+1) * z0h1d )
-!
-!--                In the borderline case the formula for stable stratification 
-!--                must be applied, because otherwise a zero division would
-!--                occur in the argument of the logarithm.
-                   IF ( a == 0.0_wp  .OR.  b == 0.0_wp )  THEN
-                      ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /       &
-                             ( LOG( zu(nzb+1) / z0h1d ) +                      &
-                               5.0_wp * rif1d(nzb+1) *                         &
-                               ( zu(nzb+1) - z0h1d ) / zu(nzb+1)               &
-                             )
-                   ELSE
-                      ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /       &
-                             LOG( (a-1.0_wp) / (a+1.0_wp) *                    &
-                                  (b+1.0_wp) / (b-1.0_wp) )
-                   ENDIF
+
+                   ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /          &
+                          LOG( (a-1.0_wp) / (a+1.0_wp) *                       &
+                               (b+1.0_wp) / (b-1.0_wp) )
                 ENDIF
 
@@ -557,6 +555,7 @@
                    pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) )
                    flux = ( ( pt_init(k+1) - pt_init(k-1) )                    &
-                            + 0.61_wp * pt_init(k)                             &
-                            * ( q_init(k+1) - q_init(k-1) )                    &
+                            + 0.61_wp                                          &
+                            * (   pt_init(k+1) * q_init(k+1)                   &
+                                - pt_init(k-1) * q_init(k-1) )                 &
                           ) * dd2zu(k)
                 ENDIF
@@ -600,20 +599,9 @@
                    b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /  &
                                                      zu(nzb+1) * z01d ) )
-!
-!--                In the borderline case the formula for stable stratification 
-!--                must be applied, because otherwise a zero division would
-!--                occur in the argument of the logarithm.
-                   IF ( a == 1.0_wp  .OR.  b == 1.0_wp )  THEN
-                      us1d = kappa * uv_total / (                              &
-                             LOG( zu(nzb+1) / z01d ) +                         &
-                             5.0_wp * rif1d(nzb+1) * ( zu(nzb+1) - z01d ) /    &
-                                                  zu(nzb+1) )
-                   ELSE
-                      us1d = kappa * uv_total / (                              &
-                                 LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) /   &
-                                      (1.0_wp+a) ) +                           &
-                                 2.0_wp * ( ATAN( b ) - ATAN( a ) )            &
-                                                )
-                   ENDIF
+                   us1d = kappa * uv_total / (                                 &
+                              LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) /      &
+                                   (1.0_wp+a) ) +                              &
+                              2.0_wp * ( ATAN( b ) - ATAN( a ) )               &
+                                             )
                 ENDIF
 
@@ -649,39 +637,30 @@
                       b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /              &
                                          zu(nzb+1) * z0h1d )
-!
-!--                   In the borderline case the formula for stable stratification 
-!--                   must be applied, because otherwise a zero division would
-!--                   occur in the argument of the logarithm.
-                      IF ( a == 1.0_wp  .OR.  b == 1.0_wp )  THEN
-                         qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /      &
-                                ( LOG( zu(nzb+1) / z0h1d ) +                   &
-                                  5.0_wp * rif1d(nzb+1) *                      &
-                                  ( zu(nzb+1) - z0h1d ) / zu(nzb+1)            &
-                                )
-                      ELSE
-                         qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /      &
-                                LOG( (a-1.0_wp) / (a+1.0_wp) *                 &
-                                     (b+1.0_wp) / (b-1.0_wp) )
-                      ENDIF
-                   ENDIF                
+                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
+                             LOG( (a-1.0_wp) / (a+1.0_wp) *                    &
+                                  (b+1.0_wp) / (b-1.0_wp) )
+                   ENDIF
                 ELSE
                    qs1d = 0.0_wp
-                ENDIF             
+                ENDIF
 
              ENDIF   !  constant_flux_layer
 
 !
-!--          Compute the diabatic mixing length
+!--          Compute the diabatic mixing length. The unstable stratification
+!--          must not be considered for l1d (km1d) as it is already considered
+!--          in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
+!--          too large.
              IF ( mixing_length_1d == 'blackadar' )  THEN
                 DO  k = nzb+1, nzt
                    IF ( rif1d(k) >= 0.0_wp )  THEN
                       l1d(k) = l_black(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
+                      l1d_diss(k) = l1d(k)
                    ELSE
-                      l1d(k) = l_black(k) *                                    &
-                               ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
+                      l1d(k) = l_black(k)
+                      l1d_diss(k) = l_black(k) *                               &
+                                    SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
                    ENDIF
-                   l1d(k) = l_black(k)
                 ENDDO
-
              ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
                 DO  k = nzb+1, nzt
@@ -694,4 +673,5 @@
                    ENDIF
                    l1d(k) = MIN( l_grid(k), l_stable )
+                   l1d_diss(k) = l1d(k)
                 ENDDO
              ENDIF
@@ -700,7 +680,5 @@
 !--          Compute the diffusion coefficients for momentum via the
 !--          corresponding Prandtl-layer relationship and according to
-!--          Prandtl-Kolmogorov, respectively. The unstable stratification is
-!--          computed via the adiabatic mixing length, for the unstability has
-!--          already been taken account of via the TKE (cf. also Diss.).
+!--          Prandtl-Kolmogorov, respectively
              IF ( constant_flux_layer )  THEN
                 IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
@@ -713,11 +691,5 @@
              ENDIF
              DO  k = nzb_diff, nzt
-                km1d(k) = c_m * SQRT( e1d(k) )
-
-                IF ( rif1d(k) >= 0.0_wp )  THEN
-                   km1d(k) = km1d(k) * l1d(k)
-                ELSE
-                   km1d(k) = km1d(k) * l_black(k)
-                ENDIF
+                km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
              ENDDO
 
@@ -882,5 +854,5 @@
     
     USE model_1d,                                                              &
-        ONLY:  dt_1d, dt_max_1d, km1d, old_dt_1d, stop_dt_1d
+        ONLY:  dt_1d, dt_max_1d, km1d, stop_dt_1d
     
     USE pegrid
@@ -922,16 +894,4 @@
     ENDIF
 
-!
-!-- A more or less simple new time step value is obtained taking only the
-!-- first two significant digits
-    div = 1000.0_wp
-    DO  WHILE ( dt_1d < div )
-       div = div / 10.0_wp
-    ENDDO
-    dt_1d = NINT( dt_1d * 100.0_wp / div ) * div / 100.0_wp
-
-    old_dt_1d = dt_1d
-
-
  END SUBROUTINE timestep_1d
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 2334)
+++ palm/trunk/SOURCE/modules.f90	(revision 2337)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! -old_dt_1d
+! +l1d_diss
+! 
+! 2326 2017-08-01 07:23:24Z gronemeier
 ! Updated variable descriptions
 ! 
@@ -1757,5 +1761,4 @@
     REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter
     REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter
-    REAL(wp) ::  old_dt_1d = 1.0E-10_wp        !< previous timestep (1d-model)
     REAL(wp) ::  qs1d                          !< characteristic humidity scale (1d-model)
     REAL(wp) ::  simulated_time_1d = 0.0_wp    !< updated simulated time (1d-model)
@@ -1775,5 +1778,6 @@
     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 (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 ::  te_e     !< tendency of e (1d-model)
