Index: /palm/trunk/SCRIPTS/document_changes
===================================================================
--- /palm/trunk/SCRIPTS/document_changes	(revision 4585)
+++ /palm/trunk/SCRIPTS/document_changes	(revision 4586)
@@ -3,7 +3,7 @@
 # This file is part of the PALM model system.
 #
-# PALM is free software: you can redistribute it and/or modify it under the 
-# terms of the GNU General Public License as published by the Free Software 
-# Foundation, either version 3 of the License, or (at your option) any later 
+# PALM is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
 # version.
 #
@@ -20,17 +20,20 @@
 # Current revisions:
 # -----------------
-# 
-# 
+#
+#
 # Former revisions:
 # -----------------
 # $Id$
+# Do not add whitespaces at current-revision section
+#
+# 4481 2020-03-31 18:55:54Z maronga
 # Bugfix for copyright updates
-# 
+#
 # 4370 2020-01-10 14:00:44Z raasch
 # script made bash compatible
-# 
+#
 # 3665 2019-01-10 08:28:24Z raasch
 # Corrected "Former revisions" section
-# 
+#
 # 2696 2017-12-14 17:12:51Z kanani
 # Change in file header (GPL part)
@@ -39,22 +42,22 @@
 # Changed the submission procedure in order to reduce the number of required
 # commits to one per change.
-# 
+#
 # 2117 2017-01-16 16:28:44Z maronga
-# 
+#
 # 1827 2016-04-07 12:12:23Z maronga
 # Added note that the script does not work on MAC OS
-# 
+#
 # 1813 2016-04-06 09:38:56Z maronga
 # Added update of the copyright statements in the file headers.
-# 
+#
 # 1810 2016-04-05 20:25:29Z maronga
 # document_changes now checks all subdirectories. Modified printing to screen.
-# 
+#
 # 1804 2016-04-05 16:30:18Z maronga
 # Removed printing of "this is an alpha version"
-# 
+#
 # 1326 2014-03-21 10:44:31Z maronga
 # Initial revision
-# 
+#
 # Description:
 # ------------
@@ -105,5 +108,5 @@
        input_dir=`pwd`
     fi
-    
+
     printf "\n"
     printf "#------------------------------------------------------------------------# \n"
@@ -111,10 +114,10 @@
     printf "|                                                                        | \n"
     printf "| This tool moves the change comments in the all PALM file headers from  | \n"
-    printf "| 'Current revisions' to 'Former revisions' and saves the time stamp.    | \n"  
+    printf "| 'Current revisions' to 'Former revisions' and saves the time stamp.    | \n"
     printf "#------------------------------------------------------------------------# \n"
-    
+
     printf "\n  *** Checking files in $input_dir and all recursive subdirectories...\n"
 
-    
+
 #
 #-- scan all (sub-)directories for files.
@@ -134,5 +137,5 @@
           continue
        fi
-       
+
        (( count_files = count_files + 1 ))
 
@@ -179,8 +182,8 @@
              then
 
-             
-                printf "\r%$(tput cols)s" " "   
+
+                printf "\r%$(tput cols)s" " "
                 printf "\r  \e[1;92m*** Comments found in $fn\e[0m\n"
-             
+
                 found_comment=true
                 cp $fn $fn~
@@ -195,5 +198,5 @@
              fi
           fi
-          
+
 #
 #--       get the timestamp from the current revision
@@ -204,5 +207,5 @@
             timestamp_string="$comment_char $timestamp"
           fi
-          
+
        done <"$fn"
 
@@ -210,9 +213,9 @@
 #
 #--    check for updates of the copyright statement
-       while read line 
-       do 
+       while read line
+       do
 
           line_tmp=""
-          line_tmp2=""     
+          line_tmp2=""
           line_tmp=${line:22:29}
           line_tmp2=${line:2:10}
@@ -226,5 +229,5 @@
              if  [[ "$year_in_file2" != "$current_year" ]]
              then
-                printf "\r%$(tput cols)s" " "   
+                printf "\r%$(tput cols)s" " "
                 printf "\r  \e[1;33m*** Copyright update required in $fn\e[0m\n"
 
@@ -240,5 +243,5 @@
        done <"$fn"
 
-       printf "\r%$(tput cols)s" " "   
+       printf "\r%$(tput cols)s" " "
        printf "\r\e[1m  *** Searched files: $count_files. Comments found: $count_changes. Copyright updates found: $count_updates\e[0m"
 
@@ -257,5 +260,5 @@
 #--       fix old time stamp
           (( line_time = line_stop + 2 ))
-          sed -i "${line_time}d" $fn      
+          sed -i "${line_time}d" $fn
           sed -i "${line_time}i$timestamp_string" $fn
 
@@ -273,6 +276,6 @@
           list_delete=${list_delete#?}
           sed -i "$list_delete" $fn
-          sed -i "${line_start}i${comment_char} " $fn
-          sed -i "${line_start}i${comment_char} " $fn
+          sed -i "${line_start}i${comment_char}" $fn
+          sed -i "${line_start}i${comment_char}" $fn
 
        fi
@@ -288,7 +291,7 @@
     then
        printf "  *** You can now proceed with\n      \e[0;91msvn commit -m 'your commit message' trunk\e[0m\n"
-       printf "  *** Please do not forget to commit your changes in the changelog at\n      \e[0;91mhttps://palm.muk.uni-hannover.de/trac/wiki/doc/tec/changelog\e[0m!\n"    
+       printf "  *** Please do not forget to commit your changes in the changelog at\n      \e[0;91mhttps://palm.muk.uni-hannover.de/trac/wiki/doc/tec/changelog\e[0m!\n"
     else
-        printf "  *** No comments found in files!\n"   
-    
+        printf "  *** No comments found in files!\n"
+
     fi
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 4585)
+++ /palm/trunk/SOURCE/header.f90	(revision 4586)
@@ -20,56 +20,60 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Renamed rif to Ri (gradient Richardson number, 1D model)
+! and zeta (= z_mo / ol, stability parameter, 3D model)
+!
+! 4573 2020-06-24 13:08:47Z oliver.maas
 ! added statement for pt_surface_heating_rate
-! 
+!
 ! 4536 2020-05-17 17:24:13Z raasch
 ! output of restart data format added
-! 
+!
 ! 4473 2020-03-25 21:04:07Z gronemeier
 ! revised message if wall_adjustment is used
-! 
+!
 ! 4444 2020-03-05 15:59:50Z raasch
 ! bugfix: cpp-directives for serial mode added
-! 
+!
 ! 4360 2020-01-07 11:25:50Z suehring
 ! Bugfix, character length too short, caused crash on NEC.
-! 
+!
 ! 4309 2019-11-26 18:49:59Z suehring
 ! replaced recycling_yshift by y_shift
-! 
+!
 ! 4301 2019-11-22 12:09:09Z oliver.maas
 !
 ! 4297 2019-11-21 10:37:50Z oliver.maas
 ! Adjusted format for simulated time and related quantities
-! 
+!
 ! 4297 2019-11-21 10:37:50Z oliver.maas
 ! adjusted message to the changed parameter recycling_yshift
-! 
+!
 ! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
-! 
+!
 ! 4223 2019-09-10 09:20:47Z gronemeier
 ! Write information about rotation angle
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 4168 2019-08-16 13:50:17Z suehring
 ! Replace function get_topography_top_index by topo_top_ind
-! 
+!
 ! 4069 2019-07-01 14:05:51Z Giersch
-! Masked output running index mid has been introduced as a local variable to 
-! avoid runtime error (Loop variable has been modified) in time_integration 
-! 
+! Masked output running index mid has been introduced as a local variable to
+! avoid runtime error (Loop variable has been modified) in time_integration
+!
 ! 4023 2019-06-12 13:20:01Z maronga
 ! Renamed "coupling start time" to "spinup time"
-! 
+!
 ! 4017 2019-06-06 12:16:46Z schwenkel
 ! unused variable removed
-! 
+!
 ! 3655 2019-01-07 16:51:22Z knoop
 ! Implementation of the PALM module interface
@@ -82,5 +86,5 @@
 ! ------------
 !> Writing a header with all important information about the current run.
-!> This subroutine is called three times, two times at the beginning 
+!> This subroutine is called three times, two times at the beginning
 !> (writing information on files RUN_CONTROL and HEADER) and one time at the
 !> end of the run, then writing additional information about CPU-usage on file
@@ -88,5 +92,5 @@
 !-----------------------------------------------------------------------------!
  SUBROUTINE header
- 
+
 
     USE arrays_3d,                                                             &
@@ -148,22 +152,22 @@
     IMPLICIT NONE
 
-    
+
     CHARACTER (LEN=2)  ::  do2d_mode           !< mode of 2D data output (xy, xz, yz)
-    
+
     CHARACTER (LEN=5)  ::  section_chr         !< string indicating grid information where to output 2D slices
-    
+
     CHARACTER (LEN=10) ::  host_chr            !< string for hostname
-    
+
     CHARACTER (LEN=16) ::  begin_chr           !< string indication start time for the data output
     CHARACTER (LEN=16) ::  coor_chr            !< dummy string
-    
+
     CHARACTER (LEN=26) ::  ver_rev             !< string for run identification
 
 #if defined( __parallel )
-    CHARACTER (LEN=32) ::  cpl_name            !< name of child domain (nesting mode only) 
+    CHARACTER (LEN=32) ::  cpl_name            !< name of child domain (nesting mode only)
 #endif
-    
+
     CHARACTER (LEN=40) ::  output_format       !< netcdf format
-        
+
     CHARACTER (LEN=70) ::  char1               !< dummy varialbe used for various strings
     CHARACTER (LEN=70) ::  char2               !< string containing informating about the advected distance in case of Galilei transformation
@@ -175,9 +179,9 @@
     CHARACTER (LEN=70) ::  do3d_chr            !< string indicating 3D output variables
     CHARACTER (LEN=70) ::  domask_chr          !< string indicating masked output variables
-    CHARACTER (LEN=70) ::  run_classification  !< string classifying type of run, e.g. nested, coupled, etc. 
-    
+    CHARACTER (LEN=70) ::  run_classification  !< string classifying type of run, e.g. nested, coupled, etc.
+
     CHARACTER (LEN=85) ::  r_upper             !< string indicating model top boundary condition for various quantities
     CHARACTER (LEN=85) ::  r_lower             !< string indicating bottom boundary condition for various quantities
-    
+
     CHARACTER (LEN=86) ::  coordinates         !< string indicating height coordinates for profile-prescribed variables
     CHARACTER (LEN=86) ::  gradients           !< string indicating gradients of profile-prescribed variables between the prescribed height coordinates
@@ -218,5 +222,5 @@
     INTEGER(iwp) ::  npe_total      !< number of total PEs in a coupler (parent + child)
 #endif
-    
+
 
     REAL(wp) ::  cpuseconds_per_simulated_second  !< CPU time (in s) per simulated second
@@ -506,11 +510,11 @@
     IF ( .NOT. ocean_mode )  THEN
        WRITE ( io, 250 )  dx, dy
-       
+
        DO i = 1, number_stretch_level_start+1
           WRITE ( io, 253 )  i, dz(i)
        ENDDO
-       
+
        WRITE( io, 251 ) (nx+1)*dx, (ny+1)*dy, zu(nzt+1)
-       
+
        IF ( ANY( dz_stretch_level_start_index < nzt+1 ) )  THEN
           WRITE( io, '(A)', advance='no') ' Vertical stretching starts at height:'
@@ -535,5 +539,5 @@
           ENDDO
        ENDIF
-       
+
     ELSE
        WRITE ( io, 250 )  dx, dy
@@ -541,7 +545,7 @@
           WRITE ( io, 253 )  i, dz(i)
        ENDDO
-       
+
        WRITE ( io, 251 ) (nx+1)*dx, (ny+1)*dy, zu(0)
-       
+
        IF ( ANY( dz_stretch_level_start_index > 0 ) )  THEN
           WRITE( io, '(A)', advance='no') ' Vertical stretching starts at height:'
@@ -603,5 +607,5 @@
        ENDDO
 
- 
+
        IF ( .NOT. large_scale_forcing )  THEN
           WRITE ( io, 426 )  TRIM( coordinates ), TRIM( temperatures ), &
@@ -620,5 +624,5 @@
     i = 1
     DO  WHILE ( ug_vertical_gradient_level_ind(i) /= -9999 )
-      
+
        WRITE (coor_chr,'(F6.2,1X)')  ug(ug_vertical_gradient_level_ind(i))
        ugcomponent = TRIM( ugcomponent ) // '  ' // TRIM( coor_chr )
@@ -672,5 +676,5 @@
           i = i + 1
        ENDIF
- 
+
     ENDDO
 
@@ -782,5 +786,5 @@
 !-- Complex terrain
     IF ( complex_terrain )  THEN
-       WRITE( io, 280 ) 
+       WRITE( io, 280 )
        IF ( turbulent_inflow )  THEN
           WRITE( io, 281 )  zu(topo_top_ind(0,0,0))
@@ -846,5 +850,5 @@
        r_upper = 'e(nzt+1) = e(nzt) = e(nzt-1)'
 
-       WRITE ( io, 301 )  'e', r_lower, r_upper       
+       WRITE ( io, 301 )  'e', r_lower, r_upper
 
     ENDIF
@@ -953,5 +957,5 @@
     WRITE ( io, 317 )  bc_lr, bc_ns
     IF ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )  THEN
-       WRITE ( io, 318 )  use_cmax, pt_damping_width, pt_damping_factor       
+       WRITE ( io, 318 )  use_cmax, pt_damping_width, pt_damping_factor
        IF ( turbulent_inflow )  THEN
           IF ( y_shift == 0 ) THEN
@@ -1009,5 +1013,5 @@
                           TRIM( gradients ), TRIM( slices )
     ELSE
-       WRITE ( io, 428 ) 
+       WRITE ( io, 428 )
     ENDIF
 
@@ -1022,5 +1026,5 @@
        i = 1
        DO  WHILE ( q_vertical_gradient_level_ind(i) /= -9999 )
-          
+
           WRITE (coor_chr,'(E8.1,4X)')  q_init(q_vertical_gradient_level_ind(i))
           temperatures = TRIM( temperatures ) // '  ' // TRIM( coor_chr )
@@ -1028,8 +1032,8 @@
           WRITE (coor_chr,'(E8.1,4X)')  q_vertical_gradient(i)
           gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
-          
+
           WRITE (coor_chr,'(I8,4X)')  q_vertical_gradient_level_ind(i)
           slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
-          
+
           WRITE (coor_chr,'(F8.1,4X)')  q_vertical_gradient_level(i)
           coordinates = TRIM( coordinates ) // '  '  // TRIM( coor_chr )
@@ -1058,5 +1062,5 @@
        i = 1
        DO  WHILE ( s_vertical_gradient_level_ind(i) /= -9999 )
-          
+
           WRITE (coor_chr,'(E8.1,4X)')  s_init(s_vertical_gradient_level_ind(i))
           temperatures = TRIM( temperatures ) // '  ' // TRIM( coor_chr )
@@ -1064,8 +1068,8 @@
           WRITE (coor_chr,'(E8.1,4X)')  s_vertical_gradient(i)
           gradients = TRIM( gradients ) // '  ' // TRIM( coor_chr )
-          
+
           WRITE (coor_chr,'(I8,4X)')  s_vertical_gradient_level_ind(i)
           slices = TRIM( slices ) // '  ' // TRIM( coor_chr )
-          
+
           WRITE (coor_chr,'(F8.1,4X)')  s_vertical_gradient_level(i)
           coordinates = TRIM( coordinates ) // '  '  // TRIM( coor_chr )
@@ -1081,5 +1085,5 @@
        WRITE ( io, 422 )  TRIM( coordinates ), TRIM( temperatures ),           &
                           TRIM( gradients ), TRIM( slices )
-    ENDIF    
+    ENDIF
 
 !
@@ -1586,8 +1590,8 @@
     ENDIF
     IF ( humidity  .AND.  q_surface_initial_change /= 0.0_wp )  THEN
-       WRITE ( io, 477 )  q_surface_initial_change       
+       WRITE ( io, 477 )  q_surface_initial_change
     ENDIF
     IF ( passive_scalar  .AND.  q_surface_initial_change /= 0.0_wp )  THEN
-       WRITE ( io, 478 )  q_surface_initial_change       
+       WRITE ( io, 478 )  q_surface_initial_change
     ENDIF
 
@@ -1727,5 +1731,5 @@
               ' Tunnel width: ', F6.2 )
 274 FORMAT (  ' Tunnel in ', A, ' direction.' / &
-              ' Tunnel height: ', F6.2, / &    
+              ' Tunnel height: ', F6.2, / &
               ' Tunnel-wall depth: ', F6.2      / &
               ' Tunnel width: ', F6.2, / &
@@ -1756,5 +1760,5 @@
              '       z_mo = ',F6.2,' m   z0 =',F7.4,' m   z0h =',F8.5,&
              ' m   kappa =',F5.2/ &
-             '       Rif value range:   ',F8.2,' <= rif <=',F6.2)
+             '       zeta value range:   ',F8.2,' <= zeta <=',F6.2)
 306 FORMAT ('       Predefined constant heatflux:   ',F9.6,' K m/s')
 307 FORMAT ('       Heatflux has a random normal distribution')
@@ -1762,5 +1766,5 @@
 309 FORMAT ('       Predefined constant salinityflux:   ',F9.6,' psu m/s')
 310 FORMAT (//'    1D-Model:'// &
-             '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
+             '       Ri value range:   ',F6.2,' <= Ri <=',F6.2)
 311 FORMAT ('       Predefined constant humidity flux: ',E10.3,' kg/kg m/s')
 312 FORMAT ('       Predefined surface humidity')
Index: /palm/trunk/SOURCE/model_1d_mod.f90
===================================================================
--- /palm/trunk/SOURCE/model_1d_mod.f90	(revision 4585)
+++ /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))
 
Index: /palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_mod.f90	(revision 4585)
+++ /palm/trunk/SOURCE/surface_mod.f90	(revision 4586)
@@ -20,9 +20,12 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! renamed Richardson flux number into gradient Richardson number (1D model)
+!
+! 4559 2020-06-11 08:51:48Z raasch
 ! File re-formatted to follow the PALM coding standard
 !
@@ -166,5 +169,5 @@
 
     USE model_1d_mod,                                                                              &
-        ONLY:  rif1d,                                                                              &
+        ONLY:  ri1d,                                                                               &
                us1d,                                                                               &
                usws1d,                                                                             &
@@ -2535,5 +2538,5 @@
        IF ( .NOT. constant_diffusion )  THEN
           IF ( constant_flux_layer )  THEN
-             surf%ol(num_h)   = surf%z_mo(num_h) / ( rif1d(nzb+1) + 1.0E-20_wp )
+             surf%ol(num_h)   = surf%z_mo(num_h) / ( ri1d(nzb+1) + 1.0E-20_wp )
              surf%us(num_h)   = us1d
              surf%usws(num_h) = usws1d
Index: /palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4585)
+++ /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4586)
@@ -20,9 +20,12 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Renamed rif to richardson_number (Richardson Flux Number to Gradient Richardson Number)
+!
+! 4510 2020-04-29 14:19:18Z raasch
 ! file re-formatted to follow the PALM coding standard
 !
@@ -4010,5 +4013,5 @@
 
     REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_stratification  !< mixing length according to stratification
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif                !< Richardson flux number
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  richardson_number  !< gradient Richardson number
 
     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
@@ -4111,5 +4114,5 @@
 
       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
-      !$ACC PRIVATE(ml_stratification, duv2_dz2, rif, dvar_dz) &
+      !$ACC PRIVATE(ml_stratification, duv2_dz2, richardson_number, dvar_dz) &
       !$ACC PRESENT(diss, e, u, v, var, wall_flags_total_0) &
       !$ACC PRESENT(dd2zu, ml_blackadar, delta)
@@ -4117,5 +4120,5 @@
           DO  j = nys, nyn
 !
-!--          Calculate Richardson-flux number
+!--          Calculate Richardson number
              IF ( use_single_reference_value )  THEN
 
@@ -4128,5 +4131,7 @@
                             + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp
 
-                   rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
+                   richardson_number(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2,        &
+                                                    -5.0_wp ),                                     &
+                                               1.0_wp )
                 ENDDO
              ELSE
@@ -4143,5 +4148,7 @@
                             + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp
 
-                   rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
+                   richardson_number(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2,           &
+                                                    -5.0_wp ),                                     &
+                                               1.0_wp )
                 ENDDO
              ENDIF
@@ -4150,8 +4157,10 @@
              !$ACC LOOP PRIVATE(k)
              DO  k = nzb+1, nzt
-                IF ( rif(k) >= 0.0_wp )  THEN
-                   ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) )
+                IF ( richardson_number(k) >= 0.0_wp )  THEN
+                   ml_stratification(k) = ml_blackadar(k)                                          &
+                                        / ( 1.0_wp + 5.0_wp * richardson_number(k) )
                 ELSE
-                   ml_stratification(k) = ml_blackadar(k) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
+                   ml_stratification(k) = ml_blackadar(k)                                          &
+                                        * SQRT( 1.0_wp - 16.0_wp * richardson_number(k) )
                 ENDIF
              ENDDO
@@ -4274,5 +4283,5 @@
 
     REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_stratification  !< mixing length according to stratification
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif                !< Richardson flux number
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  richardson_number  !< gradient Richardson number
 
     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
@@ -4348,5 +4357,5 @@
 
 !
-!--    Calculate Richardson-flux number
+!--    Calculate Richardson number
        IF ( use_single_reference_value )  THEN
           DO  k = nzb+1, nzt
@@ -4356,5 +4365,6 @@
                       + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp
 
-             rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
+             richardson_number(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),   &
+                                         1.0_wp )
           ENDDO
        ELSE
@@ -4365,5 +4375,6 @@
                       + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp
 
-             rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ), 1.0_wp )
+             richardson_number(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ),      &
+                                         1.0_wp )
           ENDDO
        ENDIF
@@ -4371,8 +4382,9 @@
 !--    Calculate diabatic mixing length using Dyer-profile functions
        DO  k = nzb+1, nzt
-          IF ( rif(k) >= 0.0_wp )  THEN
-             ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) )
+          IF ( richardson_number(k) >= 0.0_wp )  THEN
+             ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * richardson_number(k) )
           ELSE
-             ml_stratification(k) = ml_blackadar(k) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
+             ml_stratification(k) = ml_blackadar(k)                                                &
+                                  * SQRT( 1.0_wp - 16.0_wp * richardson_number(k) )
           ENDIF
 
@@ -4772,8 +4784,8 @@
     REAL(wp)     ::  var_reference       !< reference temperature
 
-    !DIR$ ATTRIBUTES ALIGN:64:: ml_local_profile, ml_stratification, rif
+    !DIR$ ATTRIBUTES ALIGN:64:: ml_local_profile, ml_stratification, richardson_number
     REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_local_profile   !< mixing length (all heights)
     REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_stratification  !< mixing length according to stratification
-    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif                !< Richardson flux number
+    REAL(wp), DIMENSION(nzb+1:nzt) ::  richardson_number  !< gradient Richardson number
 
     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
@@ -4937,5 +4949,5 @@
        !$OMP DO
        !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
-       !$ACC PRIVATE(dvar_dz, duv2_dz2, ml_stratification, ml_local_profile, rif) &
+       !$ACC PRIVATE(dvar_dz, duv2_dz2, ml_stratification, ml_local_profile, richardson_number) &
        !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, delta, ml_blackadar) &
        !$ACC PRESENT(kh, km, sums_l_l, rmask)
@@ -4943,5 +4955,5 @@
           DO  j = nysg, nyng
 !
-!--          Calculate Richardson-flux number
+!--          Calculate Richardson number
              IF ( use_single_reference_value )  THEN
                 !$ACC LOOP PRIVATE(k)
@@ -4953,5 +4965,7 @@
                             + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp
 
-                   rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
+                   richardson_number(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2,        &
+                                                    -5.0_wp ),                                     &
+                                               1.0_wp )
                 ENDDO
              ELSE
@@ -4964,5 +4978,7 @@
                             + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 + 1E-30_wp
 
-                   rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
+                   richardson_number(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2,           &
+                                                    -5.0_wp ),                                     &
+                                               1.0_wp )
                 ENDDO
              ENDIF
@@ -4972,6 +4988,7 @@
              !$ACC LOOP PRIVATE(k)
              DO  k = nzb+1, nzt
-                IF ( rif(k) >= 0.0_wp )  THEN
-                   ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) )
+                IF ( richardson_number(k) >= 0.0_wp )  THEN
+                   ml_stratification(k) = ml_blackadar(k)                                          &
+                                        / ( 1.0_wp + 5.0_wp * richardson_number(k) )
                 ELSE
                    ml_stratification(k) = ml_blackadar(k)
