Index: palm/trunk/SOURCE/model_1d_mod.f90
===================================================================
--- palm/trunk/SOURCE/model_1d_mod.f90	(revision 4449)
+++ palm/trunk/SOURCE/model_1d_mod.f90	(revision 4586)
@@ -20,9 +20,12 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! renamed Richardson flux number into gradient Richardson number
+!
+! 4449 2020-03-09 14:43:16Z suehring
 ! Set intermediate_timestep_count back to zero after 1D-model integration.
 ! This is required e.g. for initial calls of calc_mean_profile.
@@ -120,5 +123,5 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d_init !< initial mixing length (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 ::  ri1d    !< gradient Richardson 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)
@@ -174,5 +177,5 @@
 !-- Public variables
     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,  &
+            dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, ri1d, u1d,   &
             us1d, usws1d, v1d, vsws1d
 
@@ -196,5 +199,5 @@
               e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1),               &
               km1d(nzb:nzt+1), l1d(nzb:nzt+1), l1d_init(nzb:nzt+1),            &
-              l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_diss(nzb:nzt+1),       &
+              l1d_diss(nzb:nzt+1), ri1d(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),             &
@@ -211,5 +214,5 @@
        e1d = 0.0_wp; e1d_p = 0.0_wp
        kh1d = 0.0_wp; km1d = 0.0_wp
-       rif1d = 0.0_wp
+       ri1d = 0.0_wp
 !
 !--    Compute the mixing length
@@ -414,5 +417,5 @@
 !
 !--                dissipation rate
-                   IF ( rif1d(k) >= 0.0_wp )  THEN
+                   IF ( ri1d(k) >= 0.0_wp )  THEN
                       alpha_buoyancy = 1.0_wp - l1d(k) / lambda
                    ELSE
@@ -604,10 +607,10 @@
              IF ( constant_flux_layer )  THEN
 !
-!--             Compute theta* using Rif numbers of the previous time step
-                IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
+!--             Compute theta* using Ri numbers of the previous time step
+                IF ( ri1d(nzb+1) >= 0.0_wp )  THEN
 !
 !--                Stable stratification
                    ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /          &
-                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
+                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *  &
                                           ( zu(nzb+1) - z0h1d ) / zu(nzb+1)    &
                           )
@@ -615,6 +618,6 @@
 !
 !--                Unstable stratification
-                   a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
-                   b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /                 &
+                   a = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) )
+                   b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) /                  &
                        zu(nzb+1) * z0h1d )
 
@@ -630,8 +633,8 @@
              !>   2018-04-23, gronemeier
 !
-!--          Compute the Richardson-flux numbers,
+!--          Compute the gradient Richardson numbers,
 !--          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.
+!--          There, the Ri numbers of the previous time step are used.
 
              IF ( constant_flux_layer )  THEN
@@ -643,5 +646,5 @@
                    flux = ts1d + 0.61_wp * pt_init(k) * qs1d
                 ENDIF
-                rif1d(nzb+1) = zu(nzb+1) * kappa * g * flux / &
+                ri1d(nzb+1) = zu(nzb+1) * kappa * g * flux /                   &
                                ( pt_0 * ( us1d**2 + 1E-30_wp ) )
              ENDIF
@@ -659,6 +662,6 @@
                           ) * dd2zu(k)
                 ENDIF
-                IF ( rif1d(k) >= 0.0_wp )  THEN
-                   rif1d(k) = g / pt_0 * flux /                                &
+                IF ( ri1d(k) >= 0.0_wp )  THEN
+                   ri1d(k) = g / pt_0 * flux /                                 &
                               (  ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2     &
                                + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2     &
@@ -666,17 +669,17 @@
                               )
                 ELSE
-                   rif1d(k) = g / pt_0 * flux /                                &
+                   ri1d(k) = g / pt_0 * flux /                                 &
                               (  ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2     &
                                + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2     &
                                + 1E-30_wp                                      &
-                              ) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
+                              ) * ( 1.0_wp - 16.0_wp * ri1d(k) )**0.25_wp
                 ENDIF
              ENDDO
 !
-!--          Richardson-numbers must remain restricted to a realistic value
+!--          Richardson numbers must remain restricted to a realistic value
 !--          range. It is exceeded excessively for very small velocities
 !--          (u,v --> 0).
-             WHERE ( rif1d < -5.0_wp )  rif1d = -5.0_wp
-             WHERE ( rif1d > 1.0_wp )  rif1d = 1.0_wp
+             WHERE ( ri1d < -5.0_wp )  ri1d = -5.0_wp
+             WHERE ( ri1d > 1.0_wp )  ri1d = 1.0_wp
 
 !
@@ -685,9 +688,9 @@
                 uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 )
 
-                IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
+                IF ( ri1d(nzb+1) >= 0.0_wp )  THEN
 !
 !--                Stable stratification
                    us1d = kappa * uv_total / (                                 &
-                             LOG( zu(nzb+1) / z01d ) + 5.0_wp * rif1d(nzb+1) * &
+                             LOG( zu(nzb+1) / z01d ) + 5.0_wp * ri1d(nzb+1) *  &
                                               ( zu(nzb+1) - z01d ) / zu(nzb+1) &
                                              )
@@ -695,6 +698,6 @@
 !
 !--                Unstable stratification
-                   a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) ) )
-                   b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /  &
+                   a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) ) )
+                   b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) /   &
                                                      zu(nzb+1) * z01d ) )
                    us1d = kappa * uv_total / (                                 &
@@ -728,9 +731,9 @@
 !
 !--                Compute q*
-                   IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
+                   IF ( ri1d(nzb+1) >= 0.0_wp )  THEN
 !
 !--                   Stable stratification
                       qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
-                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * rif1d(nzb+1) * &
+                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *  &
                                           ( zu(nzb+1) - z0h1d ) / zu(nzb+1)    &
                           )
@@ -738,6 +741,6 @@
 !
 !--                   Unstable stratification
-                      a = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) )
-                      b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) /              &
+                      a = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) )
+                      b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) /               &
                                          zu(nzb+1) * z0h1d )
                       qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
@@ -759,11 +762,11 @@
                 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) )
+                      IF ( ri1d(k) >= 0.0_wp )  THEN
+                         l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * ri1d(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) )
+                                       SQRT( 1.0_wp - 16.0_wp * ri1d(k) )
                       ENDIF
                    ENDDO
@@ -793,10 +796,10 @@
 !--          Prandtl-Kolmogorov, respectively
              IF ( constant_flux_layer )  THEN
-                IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
+                IF ( ri1d(nzb+1) >= 0.0_wp )  THEN
                    km1d(nzb+1) = us1d * kappa * zu(nzb+1) /                    &
-                                 ( 1.0_wp + 5.0_wp * rif1d(nzb+1) )
+                                 ( 1.0_wp + 5.0_wp * ri1d(nzb+1) )
                 ELSE
                    km1d(nzb+1) = us1d * kappa * zu(nzb+1) *                    &
-                                 ( 1.0_wp - 16.0_wp * rif1d(nzb+1) )**0.25_wp
+                                 ( 1.0_wp - 16.0_wp * ri1d(nzb+1) )**0.25_wp
                 ENDIF
              ENDIF
@@ -823,8 +826,8 @@
 !--          kh = phim / phih * km
              DO  k = nzb+1, nzt
-                IF ( rif1d(k) >= 0.0_wp )  THEN
+                IF ( ri1d(k) >= 0.0_wp )  THEN
                    kh1d(k) = km1d(k)
                 ELSE
-                   kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp
+                   kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * ri1d(k) )**0.25_wp
                 ENDIF
              ENDDO
@@ -1026,6 +1029,6 @@
        WRITE ( 17, 101 )
        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), diss1d(k)
+          WRITE ( 17, 103)  k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k),      &
+                            ri1d(k), km1d(k), kh1d(k), l1d(k), diss1d(k)
        ENDDO
        WRITE ( 17, 101 )
@@ -1045,5 +1048,5 @@
 101 FORMAT ('#',111('-'))
 102 FORMAT ('#  k     zu      u          v          pt         e          ',   &
-            'rif        Km         Kh         l          diss   ')
+            'Ri         Km         Kh         l          diss   ')
 103 FORMAT (1X,I4,1X,F7.1,9(1X,E10.3))
 
