Index: palm/trunk/SOURCE/model_1d_mod.f90
===================================================================
--- palm/trunk/SOURCE/model_1d_mod.f90	(revision 4586)
+++ palm/trunk/SOURCE/model_1d_mod.f90	(revision 4677)
@@ -1,20 +1,19 @@
 !> @file model_1d_mod.f90
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! 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
-! version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
+! 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.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
 !
 ! Copyright 1997-2020 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 !
 ! Current revisions:
@@ -25,4 +24,7 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+! 4586 2020-07-01 16:16:43Z gronemeier
 ! renamed Richardson flux number into gradient Richardson number
 !
@@ -53,28 +55,25 @@
 !>        the 1D model uses different turbulence closure approaches at least if
 !>        the 3D model is set to LES-mode.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  MODULE model_1d_mod
 
-    USE arrays_3d,                                                             &
-        ONLY:  dd2zu, ddzu, ddzw, dzu, dzw, pt_init, q_init, ug, u_init,       &
-               vg, v_init, zu
-
-    USE basic_constants_and_equations_mod,                                     &
+    USE arrays_3d,                                                                                 &
+        ONLY:  dd2zu, ddzu, ddzw, dzu, dzw, pt_init, q_init, ug, u_init, vg, v_init, zu
+
+    USE basic_constants_and_equations_mod,                                                         &
         ONLY:  g, kappa, pi
 
-    USE control_parameters,                                                    &
-        ONLY:  constant_diffusion, constant_flux_layer, dissipation_1d, f,     &
-               humidity, ibc_e_b, intermediate_timestep_count,                 &
-               intermediate_timestep_count_max, km_constant,                   &
-               message_string, mixing_length_1d, prandtl_number,               &
-               roughness_length, run_description_header, simulated_time_chr,   &
-               timestep_scheme, tsc, z0h_factor
-
-    USE indices,                                                               &
+    USE control_parameters,                                                                        &
+        ONLY:  constant_diffusion, constant_flux_layer, dissipation_1d, f, humidity, ibc_e_b,      &
+               intermediate_timestep_count, intermediate_timestep_count_max, km_constant,          &
+               message_string, mixing_length_1d, prandtl_number, roughness_length,                 &
+               run_description_header, simulated_time_chr, timestep_scheme, tsc, z0h_factor
+
+    USE indices,                                                                                   &
         ONLY:  nzb, nzb_diff, nzt
 
     USE kinds
 
-    USE pegrid,                                                                &
+    USE pegrid,                                                                                    &
         ONLY:  myid
 
@@ -176,7 +175,6 @@
 !
 !-- 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, ri1d, 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, ri1d, u1d, us1d, usws1d, v1d, vsws1d
 
 
@@ -185,5 +183,5 @@
  SUBROUTINE init_1d_model
 
-    USE grid_variables,                                                        &
+    USE grid_variables,                                                                            &
         ONLY:  dx, dy
 
@@ -196,11 +194,11 @@
 !
 !-- Allocate required 1D-arrays
-    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), l1d(nzb:nzt+1), l1d_init(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),             &
-              te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(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), l1d(nzb:nzt+1), l1d_init(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),                                 &
+              te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1),                                   &
               u1d_p(nzb:nzt+1),  v1d(nzb:nzt+1), v1d_p(nzb:nzt+1) )
 
@@ -223,6 +221,5 @@
 !--       Blackadar mixing length
           IF ( f /= 0.0_wp )  THEN
-             lambda = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /        &
-                               ABS( f ) + 1E-10_wp
+             lambda = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / ABS( f ) + 1E-10_wp
           ELSE
              lambda = 30.0_wp
@@ -258,7 +255,7 @@
 
 !
-!-- Set initial horizontal velocities at the lowest grid levels to a very small
-!-- value in order to avoid too small time steps caused by the diffusion limit
-!-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!)
+!-- Set initial horizontal velocities at the lowest grid levels to a very small value in order to
+!-- avoid too small time steps caused by the diffusion limit in the initial phase of a run (at k=1,
+!-- dz/2 occurs in the limiting formula!)
     u1d(0:1)   = 0.1_wp
     u1d_p(0:1) = 0.1_wp
@@ -310,9 +307,9 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Runge-Kutta time differencing scheme for the 1D-model.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 
  SUBROUTINE time_integration_1d
@@ -335,6 +332,6 @@
 
 !
-!-- Determine the time step at the start of a 1D-simulation and
-!-- determine and printout quantities used for run control
+!-- Determine the time step at the start of a 1D-simulation and determine and printout quantities
+!-- used for run control
     dt_1d = 0.01_wp
     CALL run_control_1d
@@ -345,10 +342,8 @@
 
 !
-!--    Depending on the timestep scheme, carry out one or more intermediate
-!--    timesteps
+!--    Depending on the timestep scheme, carry out one or more intermediate timesteps
 
        intermediate_timestep_count = 0
-       DO  WHILE ( intermediate_timestep_count < &
-                   intermediate_timestep_count_max )
+       DO  WHILE ( intermediate_timestep_count < intermediate_timestep_count_max )
 
           intermediate_timestep_count = intermediate_timestep_count + 1
@@ -357,6 +352,5 @@
 
 !
-!--       Compute all tendency terms. If a constant-flux 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
 
@@ -365,13 +359,13 @@
 !
 !--          u-component
-             te_u(k) =  f * ( v1d(k) - vg(k) ) + ( &
-                              kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) &
-                            - kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k)   &
+             te_u(k) =  f * ( v1d(k) - vg(k) ) + (                                                 &
+                              kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1)                             &
+                            - kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k)                               &
                                                  ) * ddzw(k)
 !
 !--          v-component
-             te_v(k) = -f * ( u1d(k) - ug(k) ) + (                     &
-                              kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) &
-                            - kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k)   &
+             te_v(k) = -f * ( u1d(k) - ug(k) ) + (                                                 &
+                              kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1)                             &
+                            - kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k)                               &
                                                  ) * ddzw(k)
           ENDDO
@@ -387,29 +381,28 @@
                 ELSE
                    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+1) * q_init(k+1) -           &
-                                        pt_init(k-1) * q_init(k-1)   )         &
+                   flux = ( ( pt_init(k+1) - pt_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
 
 !
-!--             Calculate dissipation rate if no prognostic equation is used for
-!--             dissipation rate
+!--             Calculate dissipation rate if no prognostic equation is used for dissipation rate.
                 IF ( dissipation_1d == 'detering' )  THEN
                    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) &
-                               ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
+                   diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(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&
-                                    )                                          &
-                                    - 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) / sig_e           &
+                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) / sig_e                               &
                                    - diss1d(k)
 
@@ -420,20 +413,20 @@
                       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    ) )  &
+                      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 + 1.0_wp
-                   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_2 - c_1) * l1d(k) / lambda )  &
-                                  - g / pt_0 * kh1d(k) * flux * c_3            &
-                                  - 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)                           &
+                   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_2 - c_1) * l1d(k) / lambda )                      &
+                                  - g / pt_0 * kh1d(k) * flux * c_3                                &
+                                  - 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
 
@@ -445,6 +438,6 @@
 !
 !--       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
+!--       Finite differences of the momentum fluxes are computed using half the normal grid length
+!--       (2.0*ddzw(k)) for the sake of enhanced accuracy
           IF ( constant_flux_layer )  THEN
 
@@ -457,17 +450,16 @@
              ELSE
                 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+1) * q_init(k+1) -              &
-                                     pt_init(k-1) * q_init(k-1)   )            &
+                flux = ( ( pt_init(k+1) - pt_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
 
 !
-!--          Calculate dissipation rate if no prognostic equation is used for
-!--          dissipation rate
+!--          Calculate dissipation rate if no prognostic equation is used for dissipation rate.
              IF ( dissipation_1d == 'detering' )  THEN
                 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) )  &
+                diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) )                      &
                             * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
              ENDIF
@@ -475,11 +467,11 @@
 !
 !--          u-component
-             te_u(k) = f * ( v1d(k) - vg(k) ) + (                              &
-                       kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d       &
+             te_u(k) = f * ( v1d(k) - vg(k) ) + (                                                  &
+                       kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d                           &
                                                 ) * 2.0_wp * ddzw(k)
 !
 !--          v-component
-             te_v(k) = -f * ( u1d(k) - ug(k) ) + (                             &
-                       kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d       &
+             te_v(k) = -f * ( u1d(k) - ug(k) ) + (                                                 &
+                       kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d                           &
                                                  ) * 2.0_wp * ddzw(k)
 !
@@ -490,12 +482,12 @@
                 !>   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&
-                                    )                                          &
-                                    - 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) / sig_e           &
+                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) / sig_e                               &
                                    - diss1d(k)
              ENDIF
@@ -507,8 +499,6 @@
           DO  k = nzb+1, nzt
 
-             u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + &
-                                           tsc(3) * te_um(k) )
-             v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + &
-                                           tsc(3) * te_vm(k) )
+             u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + tsc(3) * te_um(k) )
+             v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + tsc(3) * te_vm(k) )
 
           ENDDO
@@ -516,18 +506,15 @@
 
              DO  k = nzb+1, nzt
-                e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + &
-                                              tsc(3) * te_em(k) )
+                e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + tsc(3) * te_em(k) )
              ENDDO
 
 !
-!--          Eliminate negative TKE values, which can result from the
-!--          integration due to numerical inaccuracies. In such cases the TKE
-!--          value is reduced to 10 percent of its old value.
+!--          Eliminate negative TKE values, which can result from the integration due to numerical
+!--          inaccuracies. In such cases the TKE 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) )
+                   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
@@ -556,6 +543,5 @@
                 ENDIF
 
-             ELSEIF ( intermediate_timestep_count < &
-                         intermediate_timestep_count_max )  THEN
+             ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
 
                 DO  k = nzb+1, nzt
@@ -570,6 +556,5 @@
                    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
@@ -581,8 +566,8 @@
 !
 !--       Boundary conditions for the prognostic variables.
-!--       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 and diss (e(nzb)=e(nzb+1)).
+!--       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 and diss (e(nzb)=e(nzb+1)).
           u1d_p(nzb) = 0.0_wp
           v1d_p(nzb) = 0.0_wp
@@ -611,7 +596,7 @@
 !
 !--                Stable stratification
-                   ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /          &
-                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *  &
-                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1)    &
+                   ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) /                              &
+                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *                      &
+                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1)                        &
                           )
                 ELSE
@@ -619,10 +604,8 @@
 !--                Unstable stratification
                    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 )
-
-                   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) )
+                   b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / zu(nzb+1) * z0h1d )
+
+                   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
 
@@ -634,6 +617,6 @@
 !
 !--          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.
+!--          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 Ri numbers of the previous time step are used.
 
@@ -646,6 +629,5 @@
                    flux = ts1d + 0.61_wp * pt_init(k) * qs1d
                 ENDIF
-                ri1d(nzb+1) = zu(nzb+1) * kappa * g * flux /                   &
-                               ( pt_0 * ( us1d**2 + 1E-30_wp ) )
+                ri1d(nzb+1) = zu(nzb+1) * kappa * g * flux / ( pt_0 * ( us1d**2 + 1E-30_wp ) )
              ENDIF
 
@@ -656,28 +638,27 @@
                 ELSE
                    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+1) * q_init(k+1)                   &
-                                - pt_init(k-1) * q_init(k-1) )                 &
+                   flux = ( ( pt_init(k+1) - pt_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
                 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     &
-                               + 1E-30_wp                                      &
+                   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                                                          &
                               )
                 ELSE
-                   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                                      &
+                   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 * ri1d(k) )**0.25_wp
                 ENDIF
              ENDDO
 !
-!--          Richardson numbers must remain restricted to a realistic value
-!--          range. It is exceeded excessively for very small velocities
-!--          (u,v --> 0).
+!--          Richardson numbers must remain restricted to a realistic value range. It is exceeded
+!--           excessively for very small velocities (u,v --> 0).
              WHERE ( ri1d < -5.0_wp )  ri1d = -5.0_wp
              WHERE ( ri1d > 1.0_wp )  ri1d = 1.0_wp
@@ -691,7 +672,7 @@
 !
 !--                Stable stratification
-                   us1d = kappa * uv_total / (                                 &
-                             LOG( zu(nzb+1) / z01d ) + 5.0_wp * ri1d(nzb+1) *  &
-                                              ( zu(nzb+1) - z01d ) / zu(nzb+1) &
+                   us1d = kappa * uv_total / ( LOG( zu(nzb+1) / z01d )                             &
+                                               + 5.0_wp * ri1d(nzb+1) * ( zu(nzb+1) - z01d )       &
+                                                 / zu(nzb+1)                                       &
                                              )
                 ELSE
@@ -699,10 +680,8 @@
 !--                Unstable stratification
                    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 / (                                 &
-                              LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) /      &
-                                   (1.0_wp+a) ) +                              &
-                              2.0_wp * ( ATAN( b ) - ATAN( a ) )               &
+                   b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / zu(nzb+1) * z01d ) )
+                   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
@@ -714,8 +693,8 @@
 
 !
-!--             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.
+!--             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_0 )**2
@@ -734,17 +713,15 @@
 !
 !--                   Stable stratification
-                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /         &
-                          ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *  &
-                                          ( zu(nzb+1) - z0h1d ) / zu(nzb+1)    &
-                          )
+                      qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) /                             &
+                             ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *                   &
+                                             ( zu(nzb+1) - z0h1d ) / zu(nzb+1)                     &
+                             )
                    ELSE
 !
 !--                   Unstable stratification
                       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) ) /         &
-                             LOG( (a-1.0_wp) / (a+1.0_wp) *                    &
-                                  (b+1.0_wp) / (b-1.0_wp) )
+                      b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / zu(nzb+1) * z0h1d )
+                      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
@@ -755,8 +732,7 @@
 
 !
-!--          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.
+!--          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 ( dissipation_1d /= 'prognostic' )  THEN
                 IF ( mixing_length_1d == 'blackadar' )  THEN
@@ -767,6 +743,5 @@
                       ELSE
                          l1d(k) = l1d_init(k)
-                         l1d_diss(k) = l1d_init(k) *                           &
-                                       SQRT( 1.0_wp - 16.0_wp * ri1d(k) )
+                         l1d_diss(k) = l1d_init(k) * SQRT( 1.0_wp - 16.0_wp * ri1d(k) )
                       ENDIF
                    ENDDO
@@ -775,6 +750,6 @@
                       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
+                         l_stable = 0.76_wp * SQRT( e1d(k) )                                       &
+                                    / SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
                       ELSE
                          l_stable = l1d_init(k)
@@ -786,19 +761,17 @@
              ELSE
                 DO  k = nzb+1, nzt
-                   l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) )                   &
-                          / ( diss1d(k) + 1.0E-30_wp )
+                   l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / ( diss1d(k) + 1.0E-30_wp )
                 ENDDO
              ENDIF
 
 !
-!--          Compute the diffusion coefficients for momentum via the
-!--          corresponding Prandtl-layer relationship and according to
-!--          Prandtl-Kolmogorov, respectively
+!--          Compute the diffusion coefficients for momentum via the corresponding Prandtl-layer
+!--          relationship and according to Prandtl-Kolmogorov, respectively
              IF ( constant_flux_layer )  THEN
                 IF ( ri1d(nzb+1) >= 0.0_wp )  THEN
-                   km1d(nzb+1) = us1d * kappa * zu(nzb+1) /                    &
+                   km1d(nzb+1) = us1d * kappa * zu(nzb+1) /                                        &
                                  ( 1.0_wp + 5.0_wp * ri1d(nzb+1) )
                 ELSE
-                   km1d(nzb+1) = us1d * kappa * zu(nzb+1) *                    &
+                   km1d(nzb+1) = us1d * kappa * zu(nzb+1) *                                        &
                                  ( 1.0_wp - 16.0_wp * ri1d(nzb+1) )**0.25_wp
                 ENDIF
@@ -823,6 +796,5 @@
 
 !
-!--          Compute the diffusion coefficient for heat via the relationship
-!--          kh = phim / phih * km
+!--          Compute the diffusion coefficient for heat via the relationship kh = phim / phih * km
              DO  k = nzb+1, nzt
                 IF ( ri1d(k) >= 0.0_wp )  THEN
@@ -865,6 +837,6 @@
     ENDDO   ! time loop
 !
-!-- Set intermediate_timestep_count back to zero. This is required e.g. for
-!-- initial calls of calc_mean_profile.
+!-- Set intermediate_timestep_count back to zero. This is required e.g. for initial calls of
+!-- calc_mean_profile.
     intermediate_timestep_count = 0
 
@@ -872,9 +844,9 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Compute and print out quantities for run control of the 1D model.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 
  SUBROUTINE run_control_1d
@@ -922,6 +894,6 @@
        alpha = alpha / ( 2.0_wp * pi ) * 360.0_wp
 
-       WRITE ( 15, 101 )  current_timestep_number_1d, simulated_time_chr, &
-                          dt_1d, umax, vmax, us1d, alpha, energy
+       WRITE ( 15, 101 )  current_timestep_number_1d, simulated_time_chr, dt_1d, umax, vmax, us1d, &
+                          alpha, energy
 !
 !--    Write buffer contents to disc immediately
@@ -932,8 +904,8 @@
 !
 !-- formats
-100 FORMAT (///'1D run control output:'/ &
-              &'------------------------------'// &
-           &'ITER.   HH:MM:SS    DT      UMAX   VMAX    U*   ALPHA   ENERG.'/ &
-           &'-------------------------------------------------------------')
+100 FORMAT (///'1D run control output:'/                                                           &
+               '------------------------------'//                                                  &
+            'ITER.   HH:MM:SS    DT      UMAX   VMAX    U*   ALPHA   ENERG.'/                      &
+            '-------------------------------------------------------------')
 101 FORMAT (I7,1X,A9,1X,F6.2,2X,F6.2,1X,F6.2,1X,F6.3,2X,F5.1,2X,F7.2)
 
@@ -943,9 +915,9 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Compute the time step w.r.t. the diffusion criterion
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 
  SUBROUTINE timestep_1d
@@ -965,6 +937,6 @@
 
 !
-!-- Compute the currently feasible time step according to the diffusion
-!-- criterion. At nzb+1 the half grid length is used.
+!-- Compute the currently feasible time step according to the diffusion criterion. At nzb+1 the half
+!-- grid length is used.
     fac = 0.125
     dt_diff = dt_max_1d
@@ -985,6 +957,6 @@
        stop_dt_1d = .TRUE.
 
-       WRITE( message_string, * ) 'timestep has exceeded the lower limit&',    &
-                                  'dt_1d = ',dt_1d,' s   simulation stopped!'
+       WRITE( message_string, * ) 'timestep has exceeded the lower limit&', 'dt_1d = ',dt_1d,      &
+                                  ' s   simulation stopped!'
        CALL message( 'timestep_1d', 'PA0192', 1, 2, 0, 6, 0 )
 
@@ -995,9 +967,9 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> List output of profiles from the 1D-model
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 
  SUBROUTINE print_1d_model
@@ -1029,6 +1001,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),      &
-                            ri1d(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 )
