Index: palm/trunk/SOURCE/ocean_mod.f90
===================================================================
--- palm/trunk/SOURCE/ocean_mod.f90	(revision 4796)
+++ palm/trunk/SOURCE/ocean_mod.f90	(revision 4797)
@@ -1,20 +1,19 @@
 !> @file ocean_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 2017-2020 Leibniz Universitaet Hannover
-!--------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 !
 ! Current revisions:
@@ -25,4 +24,7 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+! 4768 2020-11-02 19:11:23Z suehring
 ! Enable 3D data output also with 64-bit precision
 !
@@ -32,51 +34,49 @@
 ! 4671 2020-09-09 20:27:58Z pavelkrc
 ! Implementation of downward facing USM and LSM surfaces
-! 
+!
 ! 4535 2020-05-15 12:07:23Z raasch
 ! bugfix for restart data format query
-! 
+!
 ! 4517 2020-05-03 14:29:30Z raasch
 ! added restart with MPI-IO for reading local arrays
-! 
+!
 ! 4495 2020-04-13 20:11:20Z raasch
 ! restart data handling with MPI-IO added
-! 
+!
 ! 4481 2020-03-31 18:55:54Z maronga
 ! vector directives added to force vectorization on Intel19 compiler
-! 
+!
 ! 4346 2019-12-18 11:55:56Z motisi
-! Introduction of wall_flags_total_0, which currently sets bits based on static
-! topography information used in wall_flags_static_0
-! 
+! Introduction of wall_flags_total_0, which currently sets bits based on static topography
+! information used in wall_flags_static_0
+!
 ! 4329 2019-12-10 15:46:36Z motisi
 ! Renamed wall_flags_0 to wall_flags_static_0
-! 
+!
 ! 4272 2019-10-23 15:18:57Z schwenkel
-! Further modularization of boundary conditions: moved boundary conditions to
-! respective modules
+! Further modularization of boundary conditions: moved boundary conditions to respective modules
 !
 ! 4196 2019-08-29 11:02:06Z gronemeier
 ! Consider rotation of model domain for calculating the Stokes drift
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 4110 2019-07-22 17:05:21Z suehring
-! Pass integer flag array as well as boundary flags to WS scalar advection 
-! routine
-! 
+! Pass integer flag array as well as boundary flags to WS scalar advection routine
+!
 ! 4109 2019-07-22 17:00:34Z suehring
 ! implemented ocean_actions
-! 
+!
 ! 3767 2019-02-27 08:18:02Z raasch
-! unused variable for file index and tmp_2d removed from rrd-subroutine parameter
-! list
-! 
+! unused variable for file index and tmp_2d removed from rrd-subroutine parameter list
+!
 ! 3719 2019-02-06 13:10:18Z kanani
-! Changed log_point to log_point_s, otherwise this overlaps with 
-! 'all progn.equations' cpu measurement.
-! 
+! Changed log_point to log_point_s, otherwise this overlaps with 'all progn.equations'
+! cpu measurement.
+!
 ! 3684 2019-01-20 20:20:58Z knoop
 ! nopointer option removed
+!
 ! 3294 2018-10-01 02:37:10Z raasch
 ! initial revision
@@ -89,30 +89,30 @@
 ! Description:
 ! ------------
-!> This module contains relevant code for PALM's ocean mode, e.g. the
-!> prognostic equation for salinity, the equation of state for seawater,
-!> the Craik Leibovich force (Stokes force), and wave breaking effects
-!------------------------------------------------------------------------------!
+!> This module contains relevant code for PALM's ocean mode, e.g. the prognostic equation for
+!> salinity, the equation of state for seawater, the Craik Leibovich force (Stokes force), and wave
+!> breaking effects
+!--------------------------------------------------------------------------------------------------!
  MODULE ocean_mod
- 
-
-    USE arrays_3d,                                                             &
-        ONLY:  prho, prho_1, rho_ocean, rho_1, sa, sa_init, sa_1, sa_2, sa_3,  &
-               sa_p, tsa_m, flux_l_sa, flux_s_sa, diss_l_sa, diss_s_sa
-
-    USE control_parameters,                                                    &
-        ONLY:  atmos_ocean_sign, bottom_salinityflux,                          &
-               constant_top_salinityflux, restart_data_format_output, ocean_mode, top_salinityflux,&
-               wall_salinityflux, loop_optimization, ws_scheme_sca
+
+
+    USE arrays_3d,                                                                                 &
+        ONLY:  diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, prho, prho_1, rho_1, rho_ocean, sa,     &
+               sa_1, sa_2, sa_3, sa_init, sa_p, tsa_m
+
+    USE control_parameters,                                                                        &
+        ONLY:  atmos_ocean_sign, bottom_salinityflux, constant_top_salinityflux, loop_optimization,&
+               ocean_mode, restart_data_format_output, top_salinityflux, wall_salinityflux,        &
+               ws_scheme_sca
 
     USE kinds
 
-    USE pegrid,                                                                &
+    USE pegrid,                                                                                    &
         ONLY:  threads_per_task
 
-    USE statistics,                                                            &
+    USE statistics,                                                                                &
         ONLY:  sums_wssas_ws_l
 
-    USE indices,                                                               &
-        ONLY:  advc_flags_s, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0, nbgp
+    USE indices,                                                                                   &
+        ONLY:  advc_flags_s, nbgp, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
 
     USE restart_data_mpi_io_mod,                                                                   &
@@ -120,7 +120,6 @@
                wrd_mpi_io_global_array
 
-    USE surface_mod,                                                           &
-        ONLY:  bc_h, surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v,           &
-               surf_usm_h, surf_usm_v
+    USE surface_mod,                                                                               &
+        ONLY:  bc_h, surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
 
        USE exchange_horiz_mod,                                                                     &
@@ -132,13 +131,12 @@
     CHARACTER (LEN=20) ::  bc_sa_t = 'neumann'  !< namelist parameter
 
-    INTEGER(iwp) ::  ibc_sa_t   !< integer flag for bc_sa_t
-    INTEGER(iwp) ::  iran_ocean = -1234567  !< random number used for wave breaking effects
-
+    INTEGER(iwp) ::  ibc_sa_t                                    !< integer flag for bc_sa_t
+    INTEGER(iwp) ::  iran_ocean = -1234567                       !< random number used for wave breaking effects
     INTEGER(iwp) ::  sa_vertical_gradient_level_ind(10) = -9999  !< grid index values of sa_vertical_gradient_level(s)
 
     LOGICAL ::  salinity = .TRUE.             !< switch for using salinity
     LOGICAL ::  stokes_force = .FALSE.        !< switch to switch on the Stokes force
+    LOGICAL ::  surface_cooling_switched_off = .FALSE.  !< variable to check if surface heat flux has been switched off
     LOGICAL ::  wave_breaking = .FALSE.       !< switch to switch on wave breaking effects
-    LOGICAL ::  surface_cooling_switched_off = .FALSE.  !< variable to check if surface heat flux has been switched off
 
     REAL(wp) ::  alpha_wave_breaking = 3.0_wp !< coefficient for wave breaking generated turbulence from Noh et al. (2004), JPO
@@ -153,28 +151,28 @@
     REAL(wp) ::  u_star_wave_breaking        !< to store the absolute value of friction velocity at the ocean surface
 
-    REAL(wp), DIMENSION(12), PARAMETER ::  nom =                               &
-                          (/ 9.9984085444849347D2,   7.3471625860981584D0,     &
-                            -5.3211231792841769D-2,  3.6492439109814549D-4,    &
-                             2.5880571023991390D0,  -6.7168282786692354D-3,    &
-                             1.9203202055760151D-3,  1.1798263740430364D-2,    &
-                             9.8920219266399117D-8,  4.6996642771754730D-6,    &
-                            -2.5862187075154352D-8, -3.2921414007960662D-12 /)
-                          !< constants used in equation of state for seawater
-
-    REAL(wp), DIMENSION(13), PARAMETER ::  den =                               &
-                          (/ 1.0D0,                  7.2815210113327091D-3,    &
-                            -4.4787265461983921D-5,  3.3851002965802430D-7,    &
-                             1.3651202389758572D-10, 1.7632126669040377D-3,    &
-                            -8.8066583251206474D-6, -1.8832689434804897D-10,   &
-                             5.7463776745432097D-6,  1.4716275472242334D-9,    &
-                             6.7103246285651894D-6, -2.4461698007024582D-17,   &
-                            -9.1534417604289062D-18 /)
-                          !< constants used in equation of state for seawater
+    REAL(wp), DIMENSION(12), PARAMETER ::  nom =                                                   &
+                                              (/ 9.9984085444849347D2,   7.3471625860981584D0,     &
+                                                -5.3211231792841769D-2,  3.6492439109814549D-4,    &
+                                                 2.5880571023991390D0,  -6.7168282786692354D-3,    &
+                                                 1.9203202055760151D-3,  1.1798263740430364D-2,    &
+                                                 9.8920219266399117D-8,  4.6996642771754730D-6,    &
+                                                -2.5862187075154352D-8, -3.2921414007960662D-12 /)
+                                               !< constants used in equation of state for seawater
+
+    REAL(wp), DIMENSION(13), PARAMETER ::  den =                                                   &
+                                              (/ 1.0D0,                  7.2815210113327091D-3,    &
+                                                -4.4787265461983921D-5,  3.3851002965802430D-7,    &
+                                                 1.3651202389758572D-10, 1.7632126669040377D-3,    &
+                                                -8.8066583251206474D-6, -1.8832689434804897D-10,   &
+                                                 5.7463776745432097D-6,  1.4716275472242334D-9,    &
+                                                 6.7103246285651894D-6, -2.4461698007024582D-17,   &
+                                                -9.1534417604289062D-18 /)
+                                              !< constants used in equation of state for seawater
 
     SAVE
 
-    PUBLIC ::  bc_sa_t, ibc_sa_t, prho_reference, sa_surface,                  &
-               sa_vertical_gradient, sa_vertical_gradient_level,               &
-               sa_vertical_gradient_level_ind, stokes_force, wave_breaking
+    PUBLIC ::  bc_sa_t, ibc_sa_t, prho_reference, sa_surface, sa_vertical_gradient,                &
+               sa_vertical_gradient_level, sa_vertical_gradient_level_ind, stokes_force,           &
+               wave_breaking
 
 
@@ -285,31 +283,29 @@
 !
 !-- Add INTERFACES that must be available to other modules (alphabetical order)
-    PUBLIC eqn_state_seawater, ocean_actions, ocean_check_data_output,         &
-           ocean_check_data_output_pr, ocean_check_parameters,                 &
-           ocean_data_output_2d, ocean_data_output_3d, ocean_exchange_horiz,   &
-           ocean_define_netcdf_grid, ocean_header, ocean_init,                 &
-           ocean_init_arrays, ocean_parin, ocean_prognostic_equations,         &
-           ocean_swap_timelevel, ocean_rrd_global, ocean_rrd_local,            &
-           ocean_wrd_global, ocean_wrd_local, ocean_3d_data_averaging,         &
-           ocean_boundary_conditions, stokes_drift_terms, wave_breaking_term
+    PUBLIC eqn_state_seawater, ocean_actions, ocean_check_data_output, ocean_check_data_output_pr, &
+           ocean_check_parameters, ocean_data_output_2d, ocean_data_output_3d,                     &
+           ocean_exchange_horiz, ocean_define_netcdf_grid, ocean_header, ocean_init,               &
+           ocean_init_arrays, ocean_parin, ocean_prognostic_equations, ocean_swap_timelevel,       &
+           ocean_rrd_global, ocean_rrd_local, ocean_wrd_global, ocean_wrd_local,                   &
+           ocean_3d_data_averaging, ocean_boundary_conditions, stokes_drift_terms,                 &
+           wave_breaking_term
 
 
  CONTAINS
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Equation of state for seawater as a function of potential temperature,
-!> salinity, and pressure.
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Equation of state for seawater as a function of potential temperature, salinity, and pressure.
 !> For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech.
 !> eqn_state_seawater calculates the potential density referred at hyp(0).
 !> eqn_state_seawater_func calculates density.
 !> TODO: so far, routine is not adjusted to use topography
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE eqn_state_seawater
 
-    USE arrays_3d,                                                             &
+    USE arrays_3d,                                                                                 &
         ONLY:  hyp, prho, pt_p, rho_ocean, sa_p
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nxl, nxr, nyn, nys, nzb, nzt
 
@@ -355,11 +351,11 @@
              sa2  = sa1 * sa1
 
-             pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     +       &
-                    nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 +       &
+             pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     +                           &
+                    nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 +                           &
                     nom(7)*sa2
 
-             pden = den(1)           + den(2)*pt1     + den(3)*pt2     +       &
-                    den(4)*pt3       + den(5)*pt4     + den(6)*sa1     +       &
-                    den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    +       &
+             pden = den(1)           + den(2)*pt1     + den(3)*pt2     +                           &
+                    den(4)*pt3       + den(5)*pt4     + den(6)*sa1     +                           &
+                    den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    +                           &
                     den(10)*sa15*pt2
 !
@@ -367,8 +363,8 @@
              prho(k,j,i) = pnom / pden
 
-             pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  +       &
+             pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  +                           &
                     nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
 
-             pden = pden +             den(11)*p1     + den(12)*p2*pt3 +       &
+             pden = pden +             den(11)*p1     + den(12)*p2*pt3 +                           &
                     den(13)*p3*pt1
 
@@ -380,6 +376,6 @@
 
 !
-!--       Neumann conditions are assumed at top boundary and currently also at
-!--       bottom boundary (see comment lines below)
+!--       Neumann conditions are assumed at top boundary and currently also at bottom boundary
+!--       (see comment lines below)
           prho(nzt+1,j,i)      = prho(nzt,j,i)
           rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i)
@@ -411,15 +407,15 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Same as above, but for grid point i,j
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE eqn_state_seawater_ij( i, j )
 
-    USE arrays_3d,                                                             &
+    USE arrays_3d,                                                                                 &
         ONLY:  hyp, prho, pt_p, rho_ocean, sa_p
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nzb, nzt
 
@@ -446,4 +442,5 @@
     REAL(wp) ::  sa2    !< temporary scalar user for calculating density
 
+
     DO  k = nzb+1, nzt
 !
@@ -463,10 +460,10 @@
        sa2  = sa1 * sa1
 
-       pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     +             &
+       pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     +                                 &
               nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + nom(7)*sa2
 
-       pden = den(1)           + den(2)*pt1     + den(3)*pt2     +             &
-              den(4)*pt3       + den(5)*pt4     + den(6)*sa1     +             &
-              den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    +             &
+       pden = den(1)           + den(2)*pt1     + den(3)*pt2     +                                 &
+              den(4)*pt3       + den(5)*pt4     + den(6)*sa1     +                                 &
+              den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    +                                 &
               den(10)*sa15*pt2
 !
@@ -474,7 +471,8 @@
        prho(k,j,i) = pnom / pden
 
-       pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  +             &
+       pnom = pnom             + nom(8)*p1      + nom(9)*p1*pt2  +                                 &
               nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
-       pden = pden +             den(11)*p1     + den(12)*p2*pt3 +             &
+
+       pden = pden             + den(11)*p1     + den(12)*p2*pt3 +                                 &
               den(13)*p3*pt1
 
@@ -510,9 +508,9 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Equation of state for seawater as a function
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  REAL(wp) FUNCTION eqn_state_seawater_func( p, pt, sa )
 
@@ -532,4 +530,5 @@
     REAL(wp) ::  sa2    !< temporary scalar user for calculating density
 
+
 !
 !-- Pressure is needed in dbar
@@ -549,14 +548,14 @@
 
 
-    eqn_state_seawater_func =                                                  &
-         ( nom(1)        + nom(2)*pt1       + nom(3)*pt2    + nom(4)*pt3     + &
-           nom(5)*sa     + nom(6)*sa*pt1    + nom(7)*sa2    + nom(8)*p1      + &
-           nom(9)*p1*pt2 + nom(10)*p1*sa    + nom(11)*p2    + nom(12)*p2*pt2   &
-         ) /                                                                   &
-         ( den(1)        + den(2)*pt1       + den(3)*pt2    + den(4)*pt3     + &
-           den(5)*pt4    + den(6)*sa        + den(7)*sa*pt1 + den(8)*sa*pt3  + &
-           den(9)*sa15   + den(10)*sa15*pt2 + den(11)*p1    + den(12)*p2*pt3 + &
-           den(13)*p3*pt1                                                      &
-         )
+    eqn_state_seawater_func =                                                                      &
+                             ( nom(1)        + nom(2)*pt1       + nom(3)*pt2    + nom(4)*pt3     + &
+                               nom(5)*sa     + nom(6)*sa*pt1    + nom(7)*sa2    + nom(8)*p1      + &
+                               nom(9)*p1*pt2 + nom(10)*p1*sa    + nom(11)*p2    + nom(12)*p2*pt2   &
+                             ) /                                                                   &
+                             ( den(1)        + den(2)*pt1       + den(3)*pt2    + den(4)*pt3     + &
+                               den(5)*pt4    + den(6)*sa        + den(7)*sa*pt1 + den(8)*sa*pt3  + &
+                               den(9)*sa15   + den(10)*sa15*pt2 + den(11)*p1    + den(12)*p2*pt3 + &
+                               den(13)*p3*pt1                                                      &
+                             )
 
 
@@ -564,9 +563,9 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Reads the ocean parameters namelist
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_parin
 
@@ -576,8 +575,8 @@
 
 
-    NAMELIST /ocean_parameters/  bc_sa_t, bottom_salinityflux, salinity,       &
-             sa_surface, sa_vertical_gradient, sa_vertical_gradient_level,     &
-             stokes_waveheight, stokes_wavelength, surface_cooling_spinup_time,&
-             top_salinityflux, wall_salinityflux, wave_breaking
+    NAMELIST /ocean_parameters/  bc_sa_t, bottom_salinityflux, salinity, sa_surface,               &
+                                 sa_vertical_gradient, sa_vertical_gradient_level,                 &
+                                 stokes_waveheight, stokes_wavelength, surface_cooling_spinup_time,&
+                                 top_salinityflux, wall_salinityflux, wave_breaking
 
 !
@@ -607,16 +606,15 @@
  END SUBROUTINE ocean_parin
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Check parameters routine for the ocean mode
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_check_parameters
 
-    USE control_parameters,                                                    &
-        ONLY:  coupling_char, coupling_mode, initializing_actions,             &
-               message_string, use_top_fluxes
-
-    USE pmc_interface,                                                         &
+    USE control_parameters,                                                                        &
+        ONLY:  coupling_char, coupling_mode, initializing_actions, message_string, use_top_fluxes
+
+    USE pmc_interface,                                                                             &
         ONLY:  nested_run
 
@@ -638,13 +636,12 @@
 !
 !-- Check ocean setting
-    IF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.                           &
-         TRIM( coupling_char ) == '_O' .AND.                                   &
+    IF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.  TRIM( coupling_char ) == '_O' .AND.          &
          .NOT. ocean_mode )  THEN
 
 !
-!--    Check whether an (uncoupled) atmospheric run has been declared as an
-!--    ocean run (this setting is done via palmrun-option -y)
-       message_string = 'ocean mode does not allow coupling_char = "' //       &
-                        TRIM( coupling_char ) // '" set by palmrun-option "-y"'
+!--    Check whether an (uncoupled) atmospheric run has been declared as an ocean run (this setting
+!--    is done via palmrun-option -y)
+       message_string = 'ocean mode does not allow coupling_char = "' // TRIM( coupling_char ) //  &
+                        '" set by palmrun-option "-y"'
        CALL message( 'ocean_check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
 
@@ -665,6 +662,5 @@
        ibc_sa_t = 1
     ELSE
-       message_string = 'unknown boundary condition: bc_sa_t = "' //           &
-                        TRIM( bc_sa_t ) // '"'
+       message_string = 'unknown boundary condition: bc_sa_t = "' // TRIM( bc_sa_t ) // '"'
        CALL message( 'ocean_check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
     ENDIF
@@ -673,11 +669,8 @@
 
     IF ( .NOT. salinity )  THEN
-       IF ( ( bottom_salinityflux /= 0.0_wp  .AND.                             &
-              bottom_salinityflux /= 9999999.9_wp )  .OR.                      &
-            ( top_salinityflux /= 0.0_wp     .AND.                             &
-              top_salinityflux /= 9999999.9_wp ) )                             &
+       IF ( ( bottom_salinityflux /= 0.0_wp  .AND.  bottom_salinityflux /= 9999999.9_wp )  .OR.    &
+            ( top_salinityflux /= 0.0_wp     .AND.  top_salinityflux /= 9999999.9_wp    ) )        &
        THEN
-          message_string = 'salinityflux must not be set for ocean run ' //    &
-                           'without salinity'
+          message_string = 'salinityflux must not be set for ocean run ' // 'without salinity'
           CALL message( 'ocean_check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
        ENDIF
@@ -685,18 +678,15 @@
 
     IF ( ibc_sa_t == 1  .AND.  top_salinityflux == 9999999.9_wp )  THEN
-       message_string = 'boundary condition: bc_sa_t = "' //                   &
-                        TRIM( bc_sa_t ) // '" requires to set top_salinityflux'
+       message_string = 'boundary condition: bc_sa_t = "' // TRIM( bc_sa_t ) //                    &
+                        '" requires to set top_salinityflux'
        CALL message( 'ocean_check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
     ENDIF
 
 !
-!-- A fixed salinity at the top implies Dirichlet boundary condition for
-!-- salinity. In this case specification of a constant salinity flux is
-!-- forbidden.
-    IF ( ibc_sa_t == 0  .AND.  constant_top_salinityflux  .AND.                &
-         top_salinityflux /= 0.0_wp )  THEN
-       message_string = 'boundary condition: bc_sa_t = "' //                   &
-                        TRIM( bc_sa_t ) // '" is not allowed with ' //         &
-                        'top_salinityflux /= 0.0'
+!-- A fixed salinity at the top implies Dirichlet boundary condition for salinity. In this case
+!-- specification of a constant salinity flux is forbidden.
+    IF ( ibc_sa_t == 0  .AND.  constant_top_salinityflux  .AND.  top_salinityflux /= 0.0_wp )  THEN
+       message_string = 'boundary condition: bc_sa_t = "' // TRIM( bc_sa_t ) //                    &
+                        '" is not allowed with top_salinityflux /= 0.0'
        CALL message( 'ocean_check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
     ENDIF
@@ -707,12 +697,9 @@
        stokes_force = .TRUE.
     ELSE
-       IF ( ( stokes_waveheight <= 0.0_wp .AND. stokes_wavelength > 0.0_wp ) &
-            .OR.                                                               &
-            ( stokes_waveheight > 0.0_wp .AND. stokes_wavelength <= 0.0_wp ) &
-            .OR.                                                               &
-            ( stokes_waveheight < 0.0_wp .AND. stokes_wavelength < 0.0_wp  ) ) &
+       IF ( ( stokes_waveheight <= 0.0_wp .AND. stokes_wavelength >  0.0_wp )  .OR.                &
+            ( stokes_waveheight >  0.0_wp .AND. stokes_wavelength <= 0.0_wp )  .OR.                &
+            ( stokes_waveheight <  0.0_wp .AND. stokes_wavelength <  0.0_wp ) )                    &
        THEN
-          message_string = 'wrong settings for stokes_wavelength and/or ' //   &
-                           'stokes_waveheight'
+          message_string = 'wrong settings for stokes_wavelength and/or stokes_waveheight'
           CALL message( 'ocean_check_parameters', 'PA0460', 1, 2, 0, 6, 0 )
        ENDIF
@@ -722,11 +709,11 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Check data output.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_check_data_output( var, unit )
- 
+
     IMPLICIT NONE
 
@@ -751,15 +738,15 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Check data output of profiles
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_check_data_output_pr( variable, var_count, unit, dopr_unit )
 
-    USE arrays_3d,                                                             &
+    USE arrays_3d,                                                                                 &
         ONLY:  zu, zw
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  data_output_pr
 
@@ -772,9 +759,10 @@
     IMPLICIT NONE
 
+    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
     CHARACTER (LEN=*) ::  unit      !<
     CHARACTER (LEN=*) ::  variable  !<
-    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
 
     INTEGER(iwp) ::  var_count     !<
+
 
     SELECT CASE ( TRIM( variable ) )
@@ -837,12 +825,12 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Define appropriate grid for netcdf variables.
 !> It is called out from subroutine netcdf.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
-    
+
     IMPLICIT NONE
 
@@ -854,4 +842,5 @@
     LOGICAL, INTENT(OUT) ::  found   !< flag if output variable is found
 
+
     found  = .TRUE.
 
@@ -860,6 +849,6 @@
     SELECT CASE ( TRIM( var ) )
 
-       CASE ( 'rho_sea_water', 'rho_sea_water_xy', 'rho_sea_water_xz', &
-              'rho_sea_water_yz', 'sa', 'sa_xy', 'sa_xz', 'sa_yz' )
+       CASE ( 'rho_sea_water', 'rho_sea_water_xy', 'rho_sea_water_xz', 'rho_sea_water_yz', 'sa',   &
+              'sa_xy', 'sa_xz', 'sa_yz' )
           grid_x = 'x'
           grid_y = 'y'
@@ -877,22 +866,22 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Average 3D data.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_3d_data_averaging( mode, variable )
- 
-
-    USE arrays_3d,                                                             &
+
+
+    USE arrays_3d,                                                                                 &
         ONLY:  rho_ocean, sa
 
-    USE averaging,                                                             &
+    USE averaging,                                                                                 &
         ONLY:  rho_ocean_av, sa_av
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  average_count_3d
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
 
@@ -905,4 +894,5 @@
     INTEGER(iwp) ::  j   !< loop index
     INTEGER(iwp) ::  k   !< loop index
+
 
     IF ( mode == 'allocate' )  THEN
@@ -936,6 +926,5 @@
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) +           &
-                                               rho_ocean(k,j,i)
+                         rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) + rho_ocean(k,j,i)
                       ENDDO
                    ENDDO
@@ -964,9 +953,9 @@
 
           CASE ( 'rho_sea_water' )
-             IF ( ALLOCATED( rho_ocean_av ) ) THEN
+             IF ( ALLOCATED( rho_ocean_av ) )  THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) /           &
+                         rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) /                               &
                                                REAL( average_count_3d, KIND=wp )
                       ENDDO
@@ -976,10 +965,9 @@
 
           CASE ( 'sa' )
-             IF ( ALLOCATED( sa_av ) ) THEN
+             IF ( ALLOCATED( sa_av ) )  THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
-                         sa_av(k,j,i) = sa_av(k,j,i) /                         &
-                                        REAL( average_count_3d, KIND=wp )
+                         sa_av(k,j,i) = sa_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                       ENDDO
                    ENDDO
@@ -994,21 +982,19 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Define 2D output variables.
-!------------------------------------------------------------------------------!
- SUBROUTINE ocean_data_output_2d( av, variable, found, grid, mode, local_pf,   &
-                                  nzb_do, nzt_do )
- 
-    USE arrays_3d,                                                             &
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE ocean_data_output_2d( av, variable, found, grid, mode, local_pf, nzb_do, nzt_do )
+
+    USE arrays_3d,                                                                                 &
         ONLY:  rho_ocean, sa
 
-    USE averaging,                                                             &
+    USE averaging,                                                                                 &
         ONLY:  rho_ocean_av, sa_av
 
-    USE indices,                                                               &
-        ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,            &
-              wall_flags_total_0
+    USE indices,                                                                                   &
+        ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_total_0
 
     IMPLICIT NONE
@@ -1031,9 +1017,9 @@
     REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
 
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< local
-       !< array to which output data is resorted to
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< local array to which output data is resorted to
 
     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
-    
+
+
     found = .TRUE.
     resorted = .FALSE.
@@ -1048,5 +1034,5 @@
              to_be_resorted => rho_ocean
           ELSE
-             IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
+             IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
                 ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 rho_ocean_av = REAL( fill_value, KIND = wp )
@@ -1059,5 +1045,5 @@
              to_be_resorted => sa
           ELSE
-             IF ( .NOT. ALLOCATED( sa_av ) ) THEN
+             IF ( .NOT. ALLOCATED( sa_av ) )  THEN
                 ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 sa_av = REAL( fill_value, KIND = wp )
@@ -1077,7 +1063,6 @@
           DO  j = nys, nyn
              DO  k = nzb_do, nzt_do
-                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
-                                   REAL( fill_value, KIND = wp ),              &
-                                   BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
+                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, KIND = wp ),     &
+                                         BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
              ENDDO
           ENDDO
@@ -1085,25 +1070,24 @@
        resorted = .TRUE.
     ENDIF
- 
+
  END SUBROUTINE ocean_data_output_2d
 
- 
-!------------------------------------------------------------------------------!
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Define 3D output variables.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
- 
-
-    USE arrays_3d,                                                             &
+
+
+    USE arrays_3d,                                                                                 &
         ONLY:  rho_ocean, sa
 
-    USE averaging,                                                             &
+    USE averaging,                                                                                 &
         ONLY:  rho_ocean_av, sa_av
 
-    USE indices,                                                               &
-        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,           &
-               wall_flags_total_0
+    USE indices,                                                                                   &
+        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_total_0
 
     IMPLICIT NONE
@@ -1129,4 +1113,5 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
 
+
     found = .TRUE.
     resorted = .FALSE.
@@ -1141,5 +1126,5 @@
              to_be_resorted => rho_ocean
           ELSE
-             IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN
+             IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
                 ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 rho_ocean_av = REAL( fill_value, KIND = wp )
@@ -1152,5 +1137,5 @@
              to_be_resorted => sa
           ELSE
-             IF ( .NOT. ALLOCATED( sa_av ) ) THEN
+             IF ( .NOT. ALLOCATED( sa_av ) )  THEN
                 ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 sa_av = REAL( fill_value, KIND = wp )
@@ -1169,7 +1154,6 @@
           DO  j = nys, nyn
              DO  k = nzb_do, nzt_do
-                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
-                                   REAL( fill_value, KIND = wp ),              &
-                                   BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
+                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, KIND = wp ),     &
+                                         BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
              ENDDO
           ENDDO
@@ -1206,15 +1190,15 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Header output for ocean parameters
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_header( io )
 
-
     IMPLICIT NONE
 
     INTEGER(iwp), INTENT(IN) ::  io   !< Unit of the output file
+
 
 !
@@ -1230,12 +1214,9 @@
     ENDIF
 
-1   FORMAT (//' Ocean settings:'/                                              &
-              ' ------------------------------------------'/)
-2   FORMAT ('    --> Craik-Leibovich vortex force and Stokes drift switched',  &
-                     ' on'/                                                    &
+1   FORMAT (//' Ocean settings:'/' ------------------------------------------'/)
+2   FORMAT ('    --> Craik-Leibovich vortex force and Stokes drift switched',' on'/                &
             '        waveheight: ',F4.1,' m   wavelength: ',F6.1,' m')
-3   FORMAT ('    --> wave breaking generated turbulence switched on'/          &
-            '        alpha:    ',F4.1/                                         &
-            '        timescale:',F5.1,' s')
+3   FORMAT ('    --> wave breaking generated turbulence switched on'/                              &
+            '        alpha:    ',F4.1/'        timescale:',F5.1,' s')
 4   FORMAT ('    --> prognostic salinity equation is switched off' )
 5   FORMAT ('    --> surface heat flux is switched off after ',F8.1,' s')
@@ -1244,20 +1225,21 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Allocate arrays and assign pointers.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_init_arrays
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nxlg, nxrg, nyn, nyng, nys, nysg, nzb, nzt
 
     IMPLICIT NONE
 
-    ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
-              rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
-              sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+    ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                               &
+              rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                &
+              sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                                 &
               sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
 
     IF (  salinity )  THEN
@@ -1282,6 +1264,5 @@
 
     prho => prho_1
-    rho_ocean  => rho_1  ! routines calc_mean_profile and diffusion_e require
-                         ! density to be a pointer
+    rho_ocean  => rho_1  ! routines calc_mean_profile and diffusion_e require density to be a pointer
 
 !
@@ -1296,36 +1277,35 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Initialization of quantities needed for the ocean mode
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_init
 
 
-    USE arrays_3d,                                                             &
-        ONLY:  dzu, dzw, hyp, pt_init, ref_state, u_stokes_zu, u_stokes_zw,    &
-               v_stokes_zu, v_stokes_zw, zu, zw
-
-    USE basic_constants_and_equations_mod,                                     &
+    USE arrays_3d,                                                                                 &
+        ONLY:  dzu, dzw, hyp, pt_init, ref_state, u_stokes_zu, u_stokes_zw, v_stokes_zu,           &
+               v_stokes_zw, zu, zw
+
+    USE basic_constants_and_equations_mod,                                                         &
         ONLY:  g
 
-    USE basic_constants_and_equations_mod,                                     &
+    USE basic_constants_and_equations_mod,                                                         &
         ONLY:  pi
 
-    USE control_parameters,                                                    &
-        ONLY:  initializing_actions, molecular_viscosity, rho_surface,         &
-               rho_reference, surface_pressure, top_momentumflux_u,            &
-               top_momentumflux_v, use_single_reference_value
-
-    USE indices,                                                               &
+    USE control_parameters,                                                                        &
+        ONLY:  initializing_actions, molecular_viscosity, rho_surface, rho_reference,              &
+               surface_pressure, top_momentumflux_u, top_momentumflux_v, use_single_reference_value
+
+    USE indices,                                                                                   &
         ONLY:  nxl, nxlg, nxrg, nyng, nys, nysg, nzb, nzt
 
     USE kinds
 
-    USE pegrid,                                                                &
+    USE pegrid,                                                                                    &
         ONLY:  myid
 
-    USE statistics,                                                            &
+    USE statistics,                                                                                &
         ONLY:  hom, statistic_regions
 
@@ -1337,11 +1317,11 @@
     INTEGER(iwp) ::  n  !< loop index
 
-    REAL(wp) ::  alpha !< angle of surface stress
-    REAL(wp) ::  dum   !< dummy argument
-    REAL(wp) ::  pt_l  !< local scalar for pt used in equation of state function
-    REAL(wp) ::  sa_l  !< local scalar for sa used in equation of state function
+    REAL(wp) ::  alpha               !< angle of surface stress
+    REAL(wp) ::  dum                 !< dummy argument
+    REAL(wp) ::  pt_l                !< local scalar for pt used in equation of state function
+    REAL(wp) ::  sa_l                !< local scalar for sa used in equation of state function
     REAL(wp) ::  velocity_amplitude  !< local scalar for amplitude of Stokes drift velocity
-    REAL(wp) ::  x     !< temporary variable to store surface stress along x
-    REAL(wp) ::  y     !< temporary variable to store surface stress along y
+    REAL(wp) ::  x                   !< temporary variable to store surface stress along x
+    REAL(wp) ::  y                   !< temporary variable to store surface stress along y
 
     REAL(wp), DIMENSION(nzb:nzt+1) ::  rho_ocean_init  !< local array for initial density
@@ -1351,17 +1331,16 @@
 
 !
-!-- In case of no restart run, calculate the inital salinity profilevcusing the
-!-- given salinity gradients
+!-- In case of no restart run, calculate the inital salinity profile using the given salinity
+!-- gradients.
     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 
        sa_init = sa_surface
 !
-!--    Last arguments gives back the gradient at top level to be used as
-!--    possible Neumann boundary condition. This is not realized for the ocean
-!--    mode, therefore a dummy argument is used.
+!--    Last arguments gives back the gradient at top level to be used as possible Neumann boundary
+!--    condition. This is not realized for the ocean mode, therefore a dummy argument is used.
        IF ( salinity )  THEN
-          CALL init_vertical_profiles( sa_vertical_gradient_level_ind,          &
-                                       sa_vertical_gradient_level,              &
-                                       sa_vertical_gradient, sa_init,           &
+          CALL init_vertical_profiles( sa_vertical_gradient_level_ind,                             &
+                                       sa_vertical_gradient_level,                                 &
+                                       sa_vertical_gradient, sa_init,                              &
                                        sa_surface, dum )
        ENDIF
@@ -1370,5 +1349,5 @@
 !
 !-- Initialize required 3d-arrays
-    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
          TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
 !
@@ -1395,11 +1374,9 @@
 
 !
-!--    Initialize new time levels (only done in order to set boundary values
-!--    including ghost points)
+!--    Initialize new time levels (only done in order to set boundary values including ghost points)
        sa_p = sa
 !
-!--    Allthough tendency arrays are set in prognostic_equations, they have
-!--    have to be predefined here because they are used (but multiplied with 0)
-!--    there before they are set.
+!--    Allthough tendency arrays are set in prognostic_equations, they have to be predefined here
+!--    because they are used (but multiplied with 0) there before they are set.
        tsa_m = 0.0_wp
 
@@ -1416,12 +1393,12 @@
 
 !
-!-- Change sign of buoyancy/stability terms because density gradient is used
-!-- instead of the potential temperature gradient to calculate the buoyancy
+!-- Change sign of buoyancy/stability terms because density gradient is used instead of the
+!-- potential temperature gradient to calculate the buoyancy.
     atmos_ocean_sign = -1.0_wp
 
 !
-!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
-!-- and the reference density (used later in buoyancy term)
-!-- First step: Calculate pressure using reference density
+!-- Calculate initial vertical profile of hydrostatic pressure (in Pa) and the reference density
+!-- (used later in buoyancy term).
+!-- First step: Calculate pressure using reference density.
     hyp(nzt+1) = surface_pressure * 100.0_wp
     hyp(nzt)   = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
@@ -1435,6 +1412,6 @@
 
 !
-!-- Second step: Iteratively calculate in situ density (based on presssure)
-!-- and pressure (based on in situ density)
+!-- Second step: Iteratively calculate in situ density (based on presssure) and pressure
+!-- (based on in situ density)
     DO  n = 1, 5
 
@@ -1456,6 +1433,5 @@
        hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1)
        DO  k = nzt-1, 0, -1
-          hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k)                 &
-                                           + rho_ocean_init(k+1) ) * dzu(k+1)
+          hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k) + rho_ocean_init(k+1) ) * dzu(k+1)
        ENDDO
 
@@ -1470,6 +1446,5 @@
        pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) )
 
-       prho_reference = prho_reference + dzu(k+1) * &
-                        eqn_state_seawater_func( 0.0_wp, pt_l, sa_l )
+       prho_reference = prho_reference + dzu(k+1) * eqn_state_seawater_func( 0.0_wp, pt_l, sa_l )
 
     ENDDO
@@ -1478,6 +1453,6 @@
 
 !
-!-- Calculate the 3d array of initial in situ and potential density,
-!-- based on the initial temperature and salinity profile
+!-- Calculate the 3d array of initial in situ and potential density, based on the initial
+!-- temperature and salinity profile.
     CALL eqn_state_seawater
 
@@ -1520,13 +1495,13 @@
        ENDIF
 
-       velocity_amplitude = ( pi * stokes_waveheight / stokes_wavelength )**2 *&
+       velocity_amplitude = ( pi * stokes_waveheight / stokes_wavelength )**2 *                    &
                             SQRT( g * stokes_wavelength / ( 2.0_wp * pi ) )
 
        DO  k = nzb, nzt
-          u_stokes_zu(k) = velocity_amplitude * COS( alpha ) *                 &
+          u_stokes_zu(k) = velocity_amplitude * COS( alpha ) *                                     &
                            EXP( 4.0_wp * pi * zu(k) / stokes_wavelength )
-          u_stokes_zw(k) = velocity_amplitude * COS( alpha ) *                 &
+          u_stokes_zw(k) = velocity_amplitude * COS( alpha ) *                                     &
                            EXP( 4.0_wp * pi * zw(k) / stokes_wavelength )
-          v_stokes_zu(k) = velocity_amplitude * SIN( alpha ) *                 &
+          v_stokes_zu(k) = velocity_amplitude * SIN( alpha ) *                                     &
                            EXP( 4.0_wp * pi * zu(k) / stokes_wavelength )
           v_stokes_zw(k) = velocity_amplitude * SIN( alpha ) *                 &
@@ -1545,15 +1520,13 @@
 !
 !--    Calculate friction velocity at ocean surface
-       u_star_wave_breaking = SQRT( SQRT( top_momentumflux_u**2 +              &
-                                          top_momentumflux_v**2 ) )
-!
-!--    Set the time scale of random forcing. The vertical grid spacing at the
-!--    ocean surface is assumed as the length scale of turbulence.
+       u_star_wave_breaking = SQRT( SQRT( top_momentumflux_u**2 + top_momentumflux_v**2 ) )
+!
+!--    Set the time scale of random forcing. The vertical grid spacing at the ocean surface is
+!--    assumed as the length scale of turbulence.
 !--    Formula follows Noh et al. (2004), JPO
-       timescale_wave_breaking = 0.1_wp * dzw(nzt) / alpha_wave_breaking /     &
-                                 u_star_wave_breaking
-!
-!--    Set random number seeds differently on the processor cores in order to
-!--    create different random number sequences
+       timescale_wave_breaking = 0.1_wp * dzw(nzt) / alpha_wave_breaking / u_star_wave_breaking
+!
+!--    Set random number seeds differently on the processor cores in order to create different
+!--    random number sequences
        iran_ocean = iran_ocean + myid
     ENDIF
@@ -1562,9 +1535,9 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for all grid points
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_actions( location )
 
@@ -1588,16 +1561,17 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Call for grid points i,j
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_actions_ij( i, j, location )
 
-
-    INTEGER(iwp),      INTENT(IN) ::  i         !< grid index in x-direction
-    INTEGER(iwp),      INTENT(IN) ::  j         !< grid index in y-direction
     CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
-    INTEGER(iwp)  ::  dummy  !< call location string
+
+    INTEGER(iwp)             ::  dummy     !< call location string
+    INTEGER(iwp), INTENT(IN) ::  i         !< grid index in x-direction
+    INTEGER(iwp), INTENT(IN) ::  j         !< grid index in y-direction
+
 
     IF ( ocean_mode )   dummy = i + j
@@ -1619,44 +1593,49 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Prognostic equation for salinity.
 !> Vector-optimized version
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_prognostic_equations
 
-    USE advec_s_bc_mod,                                                        &
+    USE advec_s_bc_mod,                                                                            &
         ONLY:  advec_s_bc
 
-    USE advec_s_pw_mod,                                                        &
+    USE advec_s_pw_mod,                                                                            &
         ONLY:  advec_s_pw
 
-    USE advec_s_up_mod,                                                        &
+    USE advec_s_up_mod,                                                                            &
         ONLY:  advec_s_up
 
-    USE advec_ws,                                                              &
+    USE advec_ws,                                                                                  &
         ONLY:  advec_s_ws
 
-    USE arrays_3d,                                                             &
+    USE arrays_3d,                                                                                 &
         ONLY:  rdf_sc, tend, tsa_m
 
-    USE control_parameters,                                                    &
-        ONLY:  bc_dirichlet_l,                                                 &
-               bc_dirichlet_n,                                                 &
-               bc_dirichlet_r,                                                 &
-               bc_dirichlet_s,                                                 &
-               bc_radiation_l,                                                 &
-               bc_radiation_n,                                                 &
-               bc_radiation_r,                                                 &
-               bc_radiation_s,                                                 &
-               dt_3d, intermediate_timestep_count,                             &
-               intermediate_timestep_count_max, scalar_advec, simulated_time,  &
-               timestep_scheme, tsc, ws_scheme_sca
-
-    USE cpulog,                                                                &
+    USE control_parameters,                                                                        &
+        ONLY:  bc_dirichlet_l,                                                                     &
+               bc_dirichlet_n,                                                                     &
+               bc_dirichlet_r,                                                                     &
+               bc_dirichlet_s,                                                                     &
+               bc_radiation_l,                                                                     &
+               bc_radiation_n,                                                                     &
+               bc_radiation_r,                                                                     &
+               bc_radiation_s,                                                                     &
+               dt_3d,                                                                              &
+               intermediate_timestep_count,                                                        &
+               intermediate_timestep_count_max,                                                    &
+               scalar_advec,                                                                       &
+               simulated_time,                                                                     &
+               timestep_scheme,                                                                    &
+               tsc,                                                                                &
+               ws_scheme_sca
+
+    USE cpulog,                                                                                    &
         ONLY:  cpu_log, log_point_s
 
-    USE diffusion_s_mod,                                                       &
+    USE diffusion_s_mod,                                                                           &
         ONLY:  diffusion_s
 
@@ -1669,8 +1648,9 @@
     REAL(wp)     ::  sbt     !< weighting factor for sub-time step
 
+
 !
 !-- Switch of the surface heat flux, if requested
     IF ( surface_cooling_spinup_time /= 999999.9_wp )  THEN
-       IF ( .NOT. surface_cooling_switched_off  .AND.                          &
+       IF ( .NOT. surface_cooling_switched_off  .AND.                                              &
             simulated_time >= surface_cooling_spinup_time )  THEN
 
@@ -1682,6 +1662,6 @@
 
 !
-!-- Compute prognostic equations for the ocean mode
-!-- First, start with salinity
+!-- Compute prognostic equations for the ocean mode.
+!-- First, start with salinity.
     IF ( salinity )  THEN
 
@@ -1709,8 +1689,8 @@
           IF ( timestep_scheme(1:5) == 'runge' )  THEN
              IF ( ws_scheme_sca )  THEN
-                CALL advec_s_ws( advc_flags_s, sa, 'sa',                       &
-                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
-                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
-                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
+                CALL advec_s_ws( advc_flags_s, sa, 'sa',                                           &
+                                 bc_dirichlet_l  .OR.  bc_radiation_l,                             &
+                                 bc_dirichlet_n  .OR.  bc_radiation_n,                             &
+                                 bc_dirichlet_r  .OR.  bc_radiation_r,                             &
                                  bc_dirichlet_s  .OR.  bc_radiation_s )
              ELSE
@@ -1722,14 +1702,14 @@
        ENDIF
 
-       CALL diffusion_s( sa,                                                   &
-                         surf_def_h(0)%sasws, surf_def_h(1)%sasws,             &
-                         surf_def_h(2)%sasws,                                  &
-                         surf_lsm_h(0)%sasws, surf_lsm_h(1)%sasws,             &
-                         surf_usm_h(0)%sasws, surf_usm_h(1)%sasws,             &
-                         surf_def_v(0)%sasws, surf_def_v(1)%sasws,             &
-                         surf_def_v(2)%sasws, surf_def_v(3)%sasws,             &
-                         surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,             &
-                         surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,             &
-                         surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,             &
+       CALL diffusion_s( sa,                                                                       &
+                         surf_def_h(0)%sasws, surf_def_h(1)%sasws,                                 &
+                         surf_def_h(2)%sasws,                                                      &
+                         surf_lsm_h(0)%sasws, surf_lsm_h(1)%sasws,                                 &
+                         surf_usm_h(0)%sasws, surf_usm_h(1)%sasws,                                 &
+                         surf_def_v(0)%sasws, surf_def_v(1)%sasws,                                 &
+                         surf_def_v(2)%sasws, surf_def_v(3)%sasws,                                 &
+                         surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,                                 &
+                         surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,                                 &
+                         surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,                                 &
                          surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
 
@@ -1743,11 +1723,9 @@
              !DIR$ IVDEP
              DO  k = nzb+1, nzt
-                sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +            &
-                                                      tsc(3) * tsa_m(k,j,i) )        &
-                                                  - tsc(5) * rdf_sc(k) *             &
-                                                    ( sa(k,j,i) - sa_init(k) )       &
-                                          )                                          &
-                                            * MERGE( 1.0_wp, 0.0_wp,                 &
-                                               BTEST( wall_flags_total_0(k,j,i), 0 ) &
+                sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) )  &
+                                            - tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) )      &
+                                          )                                                        &
+                                            * MERGE( 1.0_wp, 0.0_wp,                               &
+                                                     BTEST( wall_flags_total_0(k,j,i), 0 )         &
                                                    )
                 IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
@@ -1767,11 +1745,9 @@
                 ENDDO
              ENDDO
-          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) &
-          THEN
+          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    DO  k = nzb+1, nzt
-                      tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
-                                        5.3125_wp * tsa_m(k,j,i)
+                      tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tsa_m(k,j,i)
                    ENDDO
                 ENDDO
@@ -1793,38 +1769,42 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Prognostic equations for ocean mode (so far, salinity only)
 !> Cache-optimized version
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_prognostic_equations_ij( i, j, i_omp_start, tn )
 
-    USE advec_s_pw_mod,                                                        &
+    USE advec_s_pw_mod,                                                                            &
         ONLY:  advec_s_pw
 
-    USE advec_s_up_mod,                                                        &
+    USE advec_s_up_mod,                                                                            &
         ONLY:  advec_s_up
 
-    USE advec_ws,                                                              &
+    USE advec_ws,                                                                                  &
         ONLY:  advec_s_ws
 
-    USE arrays_3d,                                                             &
+    USE arrays_3d,                                                                                 &
         ONLY:  diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, rdf_sc, tend, tsa_m
 
-    USE control_parameters,                                                    &
-        ONLY:  bc_dirichlet_l,                                                 &
-               bc_dirichlet_n,                                                 &
-               bc_dirichlet_r,                                                 &
-               bc_dirichlet_s,                                                 &
-               bc_radiation_l,                                                 &
-               bc_radiation_n,                                                 &
-               bc_radiation_r,                                                 &
-               bc_radiation_s,                                                 &
-               dt_3d, intermediate_timestep_count,                             &
-               intermediate_timestep_count_max, simulated_time,                &
-               timestep_scheme, tsc, ws_scheme_sca
-
-    USE diffusion_s_mod,                                                       &
+    USE control_parameters,                                                                        &
+        ONLY:  bc_dirichlet_l,                                                                     &
+               bc_dirichlet_n,                                                                     &
+               bc_dirichlet_r,                                                                     &
+               bc_dirichlet_s,                                                                     &
+               bc_radiation_l,                                                                     &
+               bc_radiation_n,                                                                     &
+               bc_radiation_r,                                                                     &
+               bc_radiation_s,                                                                     &
+               dt_3d,                                                                              &
+               intermediate_timestep_count,                                                        &
+               intermediate_timestep_count_max,                                                    &
+               simulated_time,                                                                     &
+               timestep_scheme,                                                                    &
+               tsc,                                                                                &
+               ws_scheme_sca
+
+    USE diffusion_s_mod,                                                                           &
         ONLY:  diffusion_s
 
@@ -1832,6 +1812,5 @@
 
     INTEGER(iwp) ::  i             !< loop index x direction
-    INTEGER(iwp) ::  i_omp_start   !< first loop index of i-loop in calling    &
-                                   !< routine prognostic_equations
+    INTEGER(iwp) ::  i_omp_start   !< first loop index of i-loop in calling routine prognostic_equations
     INTEGER(iwp) ::  j             !< loop index y direction
     INTEGER(iwp) ::  k             !< loop index z direction
@@ -1842,5 +1821,5 @@
 !-- Switch of the surface heat flux, if requested
     IF ( surface_cooling_spinup_time /= 999999.9_wp )  THEN
-       IF ( .NOT. surface_cooling_switched_off  .AND.                          &
+       IF ( .NOT. surface_cooling_switched_off  .AND.                                              &
             simulated_time >= surface_cooling_spinup_time )  THEN
 
@@ -1852,18 +1831,17 @@
 
 !
-!-- Compute prognostic equations for the ocean mode
-!-- First, start with tendency-terms for salinity
+!-- Compute prognostic equations for the ocean mode.
+!-- First, start with tendency-terms for salinity.
     IF ( salinity )  THEN
 
        tend(:,j,i) = 0.0_wp
-       IF ( timestep_scheme(1:5) == 'runge' ) &
-       THEN
+       IF ( timestep_scheme(1:5) == 'runge' )  THEN
           IF ( ws_scheme_sca )  THEN
-             CALL advec_s_ws( advc_flags_s,                                    &
-                              i, j, sa, 'sa', flux_s_sa,  diss_s_sa, flux_l_sa,&
-                              diss_l_sa, i_omp_start, tn,                      &
-                              bc_dirichlet_l  .OR.  bc_radiation_l,            &
-                              bc_dirichlet_n  .OR.  bc_radiation_n,            &
-                              bc_dirichlet_r  .OR.  bc_radiation_r,            &
+             CALL advec_s_ws( advc_flags_s,                                                        &
+                              i, j, sa, 'sa', flux_s_sa,  diss_s_sa, flux_l_sa, diss_l_sa,         &
+                              i_omp_start, tn,                                                     &
+                              bc_dirichlet_l  .OR.  bc_radiation_l,                                &
+                              bc_dirichlet_n  .OR.  bc_radiation_n,                                &
+                              bc_dirichlet_r  .OR.  bc_radiation_r,                                &
                               bc_dirichlet_s  .OR.  bc_radiation_s )
           ELSE
@@ -1873,15 +1851,14 @@
           CALL advec_s_up( i, j, sa )
        ENDIF
-       CALL diffusion_s( i, j, sa,                                             &
-                         surf_def_h(0)%sasws, surf_def_h(1)%sasws,             &
-                         surf_def_h(2)%sasws,                                  &
-                         surf_lsm_h(0)%sasws, surf_lsm_h(1)%sasws,             &
-                         surf_usm_h(0)%sasws, surf_usm_h(1)%sasws,             &
-                         surf_def_v(0)%sasws, surf_def_v(1)%sasws,             &
-                         surf_def_v(2)%sasws, surf_def_v(3)%sasws,             &
-                         surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,             &
-                         surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,             &
-                         surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,             &
-                         surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
+       CALL diffusion_s( i, j, sa,                                                                 &
+                         surf_def_h(0)%sasws, surf_def_h(1)%sasws, surf_def_h(2)%sasws,            &
+                         surf_lsm_h(0)%sasws, surf_lsm_h(1)%sasws,                                 &
+                         surf_usm_h(0)%sasws, surf_usm_h(1)%sasws,                                 &
+                         surf_def_v(0)%sasws, surf_def_v(1)%sasws, surf_def_v(2)%sasws,            &
+                         surf_def_v(3)%sasws,             &
+                         surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, surf_lsm_v(2)%sasws,            &
+                         surf_lsm_v(3)%sasws,             &
+                         surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, surf_usm_v(2)%sasws,            &
+                         surf_usm_v(3)%sasws )
 
 !       CALL user_actions( i, j, 'sa-tendency' ) ToDo: find general solution for dependency between modules
@@ -1891,11 +1868,9 @@
        DO  k = nzb+1, nzt
 
-          sa_p(k,j,i) = sa(k,j,i) + ( dt_3d *                                  &
-                                              ( tsc(2) * tend(k,j,i) +         &
-                                                tsc(3) * tsa_m(k,j,i) )        &
-                                    - tsc(5) * rdf_sc(k)                       &
-                                             * ( sa(k,j,i) - sa_init(k) )      &
-                                    ) * MERGE( 1.0_wp, 0.0_wp,                 &
-                                         BTEST( wall_flags_total_0(k,j,i), 0 ) )
+          sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) )     &
+                                    - tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) )              &
+                                    ) * MERGE( 1.0_wp, 0.0_wp,                                     &
+                                               BTEST( wall_flags_total_0(k,j,i), 0 )               &
+                                             )
 
           IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
@@ -1910,9 +1885,7 @@
                 tsa_m(k,j,i) = tend(k,j,i)
              ENDDO
-          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) &
-          THEN
+          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
              DO  k = nzb+1, nzt
-                tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +                    &
-                                  5.3125_wp * tsa_m(k,j,i)
+                tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) + 5.3125_wp * tsa_m(k,j,i)
              ENDDO
           ENDIF
@@ -1927,9 +1900,9 @@
  END SUBROUTINE ocean_prognostic_equations_ij
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Boundary conditions for ocean model
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_boundary_conditions
 
@@ -1943,7 +1916,6 @@
 
 !
-!--    Boundary conditions for salinity
-!--    Bottom boundary: Neumann condition because salinity flux is always
-!--    given.
+!--    Boundary conditions for salinity.
+!--    Bottom boundary: Neumann condition because salinity flux is always given.
        DO  l = 0, 1
           !$OMP PARALLEL DO PRIVATE( i, j, k )
@@ -1965,9 +1937,9 @@
  END SUBROUTINE ocean_boundary_conditions
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Swapping of timelevels.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_swap_timelevel( mod_count )
 
@@ -1994,13 +1966,13 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Read module-specific global restart data (Fortran binary format).
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_rrd_global_ftn( found )
 
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY: length, restart_string
 
@@ -2063,9 +2035,9 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Read module-specific global restart data (MPI-IO).
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_rrd_global_mpi
 
@@ -2088,20 +2060,19 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Read module-specific local restart data arrays (Fortran binary format).
-!------------------------------------------------------------------------------!
- SUBROUTINE ocean_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,           &
-                                 nxr_on_file, nynf, nync, nyn_on_file, nysf,       &
-                                 nysc, nys_on_file, tmp_3d, found )
-
-    USE averaging,                                                             &
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE ocean_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,  &
+                                 nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found )
+
+    USE averaging,                                                                                 &
         ONLY:  rho_ocean_av, sa_av
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  length, restart_string
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzt
 
@@ -2127,5 +2098,7 @@
     LOGICAL, INTENT(OUT)  ::  found
 
-    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
+    REAL(wp),                                                                                      &
+       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp)    &
+       :: tmp_3d   !<
 
 
@@ -2139,11 +2112,11 @@
           ENDIF
           IF ( k == 1 )  READ ( 13 )  tmp_3d
-          rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =            &
-                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+          rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
        CASE ( 'sa' )
           IF ( k == 1 )  READ ( 13 )  tmp_3d
-          sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                      &
-                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+          sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                          &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
        CASE ( 'sa_av' )
@@ -2152,6 +2125,6 @@
           ENDIF
           IF ( k == 1 )  READ ( 13 )  tmp_3d
-          sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                   &
-                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+          sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
+                                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
        CASE DEFAULT
@@ -2163,15 +2136,15 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Read module-specific local restart data arrays (MPI-IO).
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_rrd_local_mpi
 
-    USE averaging,                                                             &
+    USE averaging,                                                                                 &
         ONLY:  rho_ocean_av, sa_av
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
 
@@ -2183,5 +2156,6 @@
     CALL rd_mpi_io_check_array( 'rho_ocean_av' , found = array_found )
     IF ( array_found )  THEN
-       IF ( .NOT. ALLOCATED( rho_ocean_av ) )  ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       IF ( .NOT. ALLOCATED( rho_ocean_av ) )                                                      &
+          ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
        CALL rrd_mpi_io( 'rho_ocean_av', rho_ocean_av )
     ENDIF
@@ -2198,9 +2172,9 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> This routine writes the respective restart data for the ocean module.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_wrd_global
 
@@ -2258,5 +2232,6 @@
        CALL wrd_mpi_io_global_array( 'sa_vertical_gradient', sa_vertical_gradient )
        CALL wrd_mpi_io_global_array( 'sa_vertical_gradient_level', sa_vertical_gradient_level )
-       CALL wrd_mpi_io_global_array( 'sa_vertical_gradient_level_ind', sa_vertical_gradient_level_ind )
+       CALL wrd_mpi_io_global_array( 'sa_vertical_gradient_level_ind',                             &
+                                     sa_vertical_gradient_level_ind )
        CALL wrd_mpi_io( 'stokes_waveheight', stokes_waveheight )
        CALL wrd_mpi_io( 'stokes_wavelength', stokes_wavelength )
@@ -2271,12 +2246,12 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> This routine writes the respective restart data for the ocean module.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE ocean_wrd_local
 
-    USE averaging,                                                             &
+    USE averaging,                                                                                 &
         ONLY:  rho_ocean_av, sa_av
 
@@ -2310,27 +2285,26 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> This routine calculates the Craik Leibovich vortex force and the additional
-!> effect of the Stokes drift on the Coriolis force
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine calculates the Craik Leibovich vortex force and the additional effect of the Stokes
+!> drift on the Coriolis force.
 !> Call for all gridpoints.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE stokes_drift_terms( component )
 
-    USE arrays_3d,                                                             &
-        ONLY:  ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu,              &
-               v_stokes_zw, w, tend
-
-    USE basic_constants_and_equations_mod,                                     &
+    USE arrays_3d,                                                                                 &
+        ONLY:  ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu, v_stokes_zw, w, tend
+
+    USE basic_constants_and_equations_mod,                                                         &
         ONLY:  pi
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  f, fs, message_string, rotation_angle
 
-    USE grid_variables,                                                        &
+    USE grid_variables,                                                                            &
         ONLY:  ddx, ddy
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nxl, nxr, nys, nysv, nyn, nzb, nzt
 
@@ -2344,4 +2318,5 @@
     REAL(wp)     ::  cos_rot_angle  !< cosine of model rotation angle
     REAL(wp)     ::  sin_rot_angle  !< sine of model rotation angle
+
 
 !
@@ -2355,9 +2330,9 @@
              DO  j = nysv, nyn
                 DO  k = nzb+1, nzt
-                   tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * (              &
-                                   0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)           &
-                                         + v(k,j,i)   - v(k,j,i-1)   ) * ddx   &
-                                 - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) )   * ddy   &
-                                                                )              &
+                   tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * (                                  &
+                                                  0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)                &
+                                                        + v(k,j,i)   - v(k,j,i-1)   ) * ddx        &
+                                                - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) )   * ddy        &
+                                                                )                                  &
                                  + f * v_stokes_zu(k)
                 ENDDO
@@ -2371,8 +2346,8 @@
              DO  j = nysv, nyn
                 DO  k = nzb+1, nzt
-                   tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * (              &
-                                   0.5 * ( v(k,j,i+1) - v(k,j,i-1) )   * ddx   &
-                                 - 0.5 * ( u(k,j,i) - u(k,j-1,i)               &
-                                         + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
+                   tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * (                                  &
+                                                  0.5 * ( v(k,j,i+1) - v(k,j,i-1) )   * ddx        &
+                                                - 0.5 * ( u(k,j,i) - u(k,j-1,i)                    &
+                                                        + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy        &
                                                                 )              &
                                  - f * u_stokes_zu(k)
@@ -2393,20 +2368,19 @@
              DO  j = nys, nyn
                 DO  k = nzb+1, nzt
-                   tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * (              &
-                                             0.5 * ( u(k+1,j,i) - u(k,j,i)     &
-                                                   + u(k+1,j,i+1) - u(k,j,i+1) &
-                                                   ) * ddzu(k+1)               &
-                                           - 0.5 * ( w(k,j,i+1) - w(k,j,i-1)   &
-                                                   ) * ddx      )              &
-                                             - v_stokes_zw(k) * (              &
-                                             0.5 * ( w(k,j+1,i) - w(k,j-1,i)   &
-                                                   ) * ddy                     &
-                                           - 0.5 * ( v(k+1,j,i) - v(k,j,i)     &
-                                                   + v(k+1,j+1,i) - v(k,j+1,i) &
-                                                   ) * ddzu(k)  )              &
-                                           + fs * (                            &
-                                               sin_rot_angle * v_stokes_zw(k)  &
-                                             + cos_rot_angle * u_stokes_zw(k)  &
-                                                  )
+                   tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) *   (                                &
+                                                   0.5 * ( u(k+1,j,i) - u(k,j,i)                   &
+                                                         + u(k+1,j,i+1) - u(k,j,i+1)               &
+                                                         ) * ddzu(k+1)                             &
+                                                 - 0.5 * ( w(k,j,i+1) - w(k,j,i-1)                 &
+                                                         ) * ddx  )                                &
+                                             - v_stokes_zw(k) *      (                             &
+                                                   0.5 * ( w(k,j+1,i) - w(k,j-1,i)                 &
+                                                         ) * ddy                                   &
+                                                 - 0.5 * ( v(k+1,j,i) - v(k,j,i)                   &
+                                                         + v(k+1,j+1,i) - v(k,j+1,i)               &
+                                                         ) * ddzu(k) )                             &
+                                              + fs * (   sin_rot_angle * v_stokes_zw(k)            &
+                                                       + cos_rot_angle * u_stokes_zw(k)            &
+                                                     )
                 ENDDO
              ENDDO
@@ -2414,6 +2388,5 @@
 
        CASE DEFAULT
-          WRITE( message_string, * ) 'wrong component of Stokes force: ',      &
-                                     component
+          WRITE( message_string, * ) 'wrong component of Stokes force: ', component
           CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 )
 
@@ -2423,28 +2396,27 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> This routine calculates the Craik Leibovich vortex force and the additional
-!> effect of the Stokes drift on the Coriolis force
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine calculates the Craik Leibovich vortex force and the additional effect of the Stokes
+!> drift on the Coriolis force.
 !> Call for gridpoints i,j.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 
  SUBROUTINE stokes_drift_terms_ij( i, j, component )
 
-    USE arrays_3d,                                                             &
-        ONLY:  ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu,              &
-               v_stokes_zw, w, tend
-
-    USE basic_constants_and_equations_mod,                                     &
+    USE arrays_3d,                                                                                 &
+        ONLY:  ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu, v_stokes_zw, w, tend
+
+    USE basic_constants_and_equations_mod,                                                         &
         ONLY:  pi
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  f, fs, message_string, rotation_angle
 
-    USE grid_variables,                                                        &
+    USE grid_variables,                                                                            &
         ONLY:  ddx, ddy
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nzb, nzt
 
@@ -2459,4 +2431,5 @@
     REAL(wp)     ::  sin_rot_angle  !< sine of model rotation angle
 
+
 !
 !-- Compute Stokes terms for the respective velocity components
@@ -2467,9 +2440,9 @@
        CASE ( 1 )
           DO  k = nzb+1, nzt
-             tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * (                    &
-                                     0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)         &
-                                           + v(k,j,i)   - v(k,j,i-1)   ) * ddx &
-                                   - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) )   * ddy &
-                                                          )                    &
+             tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * (                                        &
+                                              0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)                    &
+                                                    + v(k,j,i)   - v(k,j,i-1)   ) * ddx            &
+                                            - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) )   * ddy            &
+                                                          )                                        &
                                        + f * v_stokes_zu(k)
           ENDDO
@@ -2478,9 +2451,9 @@
        CASE ( 2 )
           DO  k = nzb+1, nzt
-             tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * (                    &
-                                     0.5 * ( v(k,j,i+1) - v(k,j,i-1) )   * ddx &
-                                   - 0.5 * ( u(k,j,i) - u(k,j-1,i)             &
-                                           + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
-                                                          )                    &
+             tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * (                                        &
+                                              0.5 * ( v(k,j,i+1) - v(k,j,i-1) )   * ddx            &
+                                            - 0.5 * ( u(k,j,i) - u(k,j-1,i)                        &
+                                                    + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy            &
+                                                          )                                        &
                                        - f * u_stokes_zu(k)
           ENDDO
@@ -2496,18 +2469,16 @@
 
           DO  k = nzb+1, nzt
-             tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * (              &
-                                     0.5 * ( u(k+1,j,i) - u(k,j,i)     &
-                                                   + u(k+1,j,i+1) - u(k,j,i+1) &
-                                                   ) * ddzu(k+1)               &
-                                           - 0.5 * ( w(k,j,i+1) - w(k,j,i-1)   &
-                                                   ) * ddx )                   &
-                                       - v_stokes_zw(k) * (                    &
-                                             0.5 * ( w(k,j+1,i) - w(k,j-1,i)   &
-                                                   ) * ddy                     &
-                                           - 0.5 * ( v(k+1,j,i) - v(k,j,i)     &
-                                                   + v(k+1,j+1,i) - v(k,j+1,i) &
-                                                   ) * ddzu(k)  )              &
-                                       + fs * ( sin_rot_angle * v_stokes_zw(k) &
-                                              + cos_rot_angle * u_stokes_zw(k) &
+             tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * ( 0.5 * ( u(k+1,j,i) - u(k,j,i)          &
+                                                                  + u(k+1,j,i+1) - u(k,j,i+1)      &
+                                                                  ) * ddzu(k+1)                    &
+                                                          - 0.5 * ( w(k,j,i+1) - w(k,j,i-1)        &
+                                                                  ) * ddx      )                   &
+                                       - v_stokes_zw(k) * ( 0.5 * ( w(k,j+1,i) - w(k,j-1,i)        &
+                                                                  ) * ddy                          &
+                                                          - 0.5 * ( v(k+1,j,i) - v(k,j,i)          &
+                                                                  + v(k+1,j+1,i) - v(k,j+1,i)      &
+                                                                  ) * ddzu(k)  )                   &
+                                       + fs * ( sin_rot_angle * v_stokes_zw(k)                     &
+                                              + cos_rot_angle * u_stokes_zw(k)                     &
                                               )
           ENDDO
@@ -2522,22 +2493,22 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> This routine calculates turbulence generated by wave breaking near the ocean
-!> surface, following a parameterization given in Noh et al. (2004), JPO
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine calculates turbulence generated by wave breaking near the ocean surface, following
+!> a parameterization given in Noh et al. (2004), JPO
 !> Call for all gridpoints.
-!> TODO: so far, this routine only works if the model time step has about the
-!>       same value as the time scale of wave breaking!
-!------------------------------------------------------------------------------!
+!> TODO: so far, this routine only works if the model time step has about the same value as the time
+!>       scale of wave breaking!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE wave_breaking_term( component )
 
-    USE arrays_3d,                                                             &
+    USE arrays_3d,                                                                                 &
         ONLY:  u_p, v_p
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  dt_3d, message_string
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nxl, nxlu, nxr, nys, nysv, nyn, nzt
 
@@ -2548,6 +2519,5 @@
     INTEGER(iwp) ::  j          !< loop index along y
 
-    REAL(wp) ::  random_gauss  !< function that creates a random number with a
-                               !< Gaussian distribution
+    REAL(wp) ::  random_gauss  !< function that creates a random number with a Gaussian distribution
 
 
@@ -2562,7 +2532,7 @@
           DO  i = nxlu, nxr
              DO  j = nys, nyn
-                u_p(nzt,j,i) = u_p(nzt,j,i) +                                  &
-                               ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
-                               * alpha_wave_breaking * u_star_wave_breaking    &
+                u_p(nzt,j,i) = u_p(nzt,j,i) +                                                      &
+                               ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )                     &
+                               * alpha_wave_breaking * u_star_wave_breaking                        &
                                / timescale_wave_breaking * dt_3d
              ENDDO
@@ -2573,7 +2543,7 @@
           DO  i = nxl, nxr
              DO  j = nysv, nyn
-                v_p(nzt,j,i) = v_p(nzt,j,i) +                                  &
-                               ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &
-                               * alpha_wave_breaking * u_star_wave_breaking    &
+                v_p(nzt,j,i) = v_p(nzt,j,i) +                                                      &
+                               ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )                     &
+                               * alpha_wave_breaking * u_star_wave_breaking                        &
                                / timescale_wave_breaking * dt_3d
              ENDDO
@@ -2581,6 +2551,5 @@
 
        CASE DEFAULT
-          WRITE( message_string, * ) 'wrong component of wave breaking: ',     &
-                                     component
+          WRITE( message_string, * ) 'wrong component of wave breaking: ', component
           CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 )
 
@@ -2590,22 +2559,22 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> This routine calculates turbulence generated by wave breaking near the ocean
-!> surface, following a parameterization given in Noh et al. (2004), JPO
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine calculates turbulence generated by wave breaking near the ocean surface, following a
+!> parameterization given in Noh et al. (2004), JPO
 !> Call for gridpoint i,j.
-!> TODO: so far, this routine only works if the model time step has about the
-!>       same value as the time scale of wave breaking!
-!------------------------------------------------------------------------------!
+!> TODO: so far, this routine only works if the model time step has about the same value as the time
+!> scale of wave breaking!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE wave_breaking_term_ij( i, j, component )
 
-    USE arrays_3d,                                                             &
+    USE arrays_3d,                                                                                 &
         ONLY:  u_p, v_p
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  dt_3d, message_string
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nzt
 
@@ -2616,6 +2585,5 @@
     INTEGER(iwp) ::  j          !< loop index along y
 
-    REAL(wp) ::  random_gauss  !< function that creates a random number with a
-                               !< Gaussian distribution
+    REAL(wp) ::  random_gauss  !< function that creates a random number with a Gaussian distribution
 
 !
@@ -2626,18 +2594,17 @@
 !--    u-/v-component
        CASE ( 1 )
-          u_p(nzt,j,i) = u_p(nzt,j,i) +                                        &
-                         ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )       &
-                         * alpha_wave_breaking * u_star_wave_breaking          &
+          u_p(nzt,j,i) = u_p(nzt,j,i) +                                                            &
+                         ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )                           &
+                         * alpha_wave_breaking * u_star_wave_breaking                              &
                          / timescale_wave_breaking * dt_3d
 
        CASE ( 2 )
-          v_p(nzt,j,i) = v_p(nzt,j,i) +                                        &
-                         ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )       &
-                         * alpha_wave_breaking * u_star_wave_breaking          &
+          v_p(nzt,j,i) = v_p(nzt,j,i) +                                                            &
+                         ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp )                           &
+                         * alpha_wave_breaking * u_star_wave_breaking                              &
                          / timescale_wave_breaking * dt_3d
 
        CASE DEFAULT
-          WRITE( message_string, * ) 'wrong component of wave breaking: ',     &
-                                     component
+          WRITE( message_string, * ) 'wrong component of wave breaking: ', component
           CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 )
 
Index: palm/trunk/SOURCE/outflow_turbulence.f90
===================================================================
--- palm/trunk/SOURCE/outflow_turbulence.f90	(revision 4796)
+++ palm/trunk/SOURCE/outflow_turbulence.f90	(revision 4797)
@@ -1,30 +1,32 @@
 !> @file outflow_turbulence.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 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.
+! 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/>.
+! 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:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id: outflow_turbulence.f90 3241 2018-09-12 15:02:00Z raasch $
+! file re-formatted to follow the PALM coding standard
+!
+! 3241 2018-09-12 15:02:00Z raasch
 ! Corrected "Former revisions" section
-! 
+!
 ! 3241 2018-09-12 15:02:00Z raasch
 ! unused variables removed
@@ -32,32 +34,31 @@
 ! 2050 2016-11-08 15:00:55Z gronemeier
 ! Initial version
-! 
+!
 !
 ! Description:
 ! ------------
-!> Routine based on inflow_turbulence.f90. Copies values of 3d data from a 2d
-!> vertical source plane (defined by outflow_source_plane) to the outflow
-!> boundary.
-!------------------------------------------------------------------------------!
+!> Routine based on inflow_turbulence.f90. Copies values of 3d data from a 2d vertical source plane
+!> (defined by outflow_source_plane) to the outflow boundary.
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE outflow_turbulence
 
-    USE arrays_3d,                                                             &
+    USE arrays_3d,                                                                                 &
         ONLY:  e, pt, q, s, u, v, w
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                        &
         ONLY:  humidity, passive_scalar, outflow_source_plane
 
-    USE cpulog,                                                                &
+    USE cpulog,                                                                                    &
         ONLY:  cpu_log, log_point
 
-    USE grid_variables,                                                        &
+    USE grid_variables,                                                                            &
         ONLY:  ddx
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nbgp, nx, nxr, nyn, nys, nyng, nysg, nzb, nzt
 
     USE kinds
 
-    USE pegrid!,                                                                &
+    USE pegrid!,                                                                                    &
         !ONLY:  comm1dx, id_outflow, id_outflow_source, ierr, myidx, status
 
@@ -70,6 +71,5 @@
     INTEGER(iwp) ::  ngp_ofv  !< number of grid points stored in outflow_val
 
-    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,7,nbgp) ::                         &
-       outflow_val            !< values to be copied to the outflow boundary
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,7,nbgp) ::  outflow_val  !< values to be copied to the outflow boundary
 
 
@@ -94,7 +94,7 @@
                 outflow_val(k,j,4,l) = pt(k,j,i)
                 outflow_val(k,j,5,l) = e(k,j,i)
-                IF ( humidity  )                                               &
+                IF ( humidity  )                                                                   &
                    outflow_val(k,j,6,l) = q(k,j,i)
-                IF ( passive_scalar )                                          &
+                IF ( passive_scalar )                                                              &
                    outflow_val(k,j,7,l) = s(k,j,i)
 
@@ -115,7 +115,7 @@
              outflow_val(k,j,4,l) = pt(k,j,i)
              outflow_val(k,j,5,l) = e(k,j,i)
-             IF ( humidity  )                                                  &
+             IF ( humidity  )                                                                      &
                 outflow_val(k,j,6,l) = q(k,j,i)
-             IF ( passive_scalar )                                             &
+             IF ( passive_scalar )                                                                 &
                 outflow_val(k,j,7,l) = s(k,j,i)
 
@@ -131,12 +131,11 @@
     IF ( myidx == id_outflow_source  .AND.  myidx /= id_outflow )  THEN
 
-       CALL MPI_SEND( outflow_val(nzb,nysg,1,1), ngp_ofv, MPI_REAL,            &
-                      id_outflow, 1, comm1dx, ierr )
+       CALL MPI_SEND( outflow_val(nzb,nysg,1,1), ngp_ofv, MPI_REAL, id_outflow, 1, comm1dx, ierr )
 
     ELSEIF ( myidx /= id_outflow_source  .AND.  myidx == id_outflow )  THEN
 
        outflow_val = 0.0_wp
-       CALL MPI_RECV( outflow_val(nzb,nysg,1,1), ngp_ofv, MPI_REAL,            &
-                      id_outflow_source, 1, comm1dx, status, ierr )
+       CALL MPI_RECV( outflow_val(nzb,nysg,1,1), ngp_ofv, MPI_REAL, id_outflow_source, 1, comm1dx, &
+                      status, ierr )
 
     ENDIF
@@ -157,7 +156,7 @@
              e(k,j,nx+1:nx+nbgp)  = MAX( e(k,j,nx+1:nx+nbgp), 0.0_wp )
 
-             IF ( humidity )                                                   &
+             IF ( humidity )                                                                       &
                 q(k,j,nx+1:nx+nbgp)  = outflow_val(k,j,6,1:nbgp)
-             IF ( passive_scalar )                                             &
+             IF ( passive_scalar )                                                                 &
                 s(k,j,nx+1:nx+nbgp)  = outflow_val(k,j,7,1:nbgp)
 
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 4796)
+++ palm/trunk/SOURCE/palm.f90	(revision 4797)
@@ -1,76 +1,78 @@
 ! !> @file palm.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:
 ! -----------------
-! 
+!
 !
 ! Former revisions:
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+! 4539 2020-05-18 14:05:17Z raasch
 ! log point name changed
-! 
+!
 ! 4535 2020-05-15 12:07:23Z raasch
 ! bugfix for restart data format query
-! 
+!
 ! 4496 2020-04-15 08:37:26Z raasch
 ! bugfix: coupling character added to restart output filename
-! 
+!
 ! 4495 2020-04-13 20:11:20Z raasch
 ! restart data handling with MPI-IO added
-! 
+!
 ! 4457 2020-03-11 14:20:43Z raasch
 ! use statement for exchange horiz added
-! 
+!
 ! 4444 2020-03-05 15:59:50Z raasch
 ! bugfix: cpp-directives for serial mode added
-! 
+!
 ! 4414 2020-02-19 20:16:04Z suehring
 ! Call to module_interface_init_numerics
-! 
+!
 ! 4400 2020-02-10 20:32:41Z suehring
 ! Add interface to initialize data output with dom
-! 
+!
 ! 4360 2020-01-07 11:25:50Z suehring
 ! implement new palm_date_time_mod
-! 
+!
 ! 4094 2019-07-12 09:24:21Z gronemeier
 ! Corrected "Former revisions" section
-! 
+!
 ! 4039 2019-06-18 10:32:41Z suehring
 ! Rename subroutines in module for diagnostic quantities
-! 
+!
 ! 4017 2019-06-06 12:16:46Z schwenkel
 ! new module for calculation and output of diagnostic quantities added
-! 
+!
 ! 3885 2019-04-11 11:29:34Z kanani
-! Changes related to global restructuring of location messages and introduction 
+! Changes related to global restructuring of location messages and introduction
 ! of additional debug messages
-! 
+!
 ! 3761 2019-02-25 15:31:42Z raasch
 ! unused variable removed
-! 
+!
 ! 3719 2019-02-06 13:10:18Z kanani
 ! Included cpu measurement for wall/soil spinup
-! 
+!
 ! 3703 2019-01-29 16:43:53Z knoop
 ! Some interface calls moved to module_interface + cleanup
-! 
+!
 ! 3648 2019-01-02 16:35:46Z suehring
 ! Rename subroutines for surface-data output
@@ -82,46 +84,45 @@
 ! Description:
 ! ------------
-!> Large-Eddy Simulation (LES) model for atmospheric and oceanic boundary-layer
-!> flows
+!> Large-Eddy Simulation (LES) model for atmospheric and oceanic boundary-layer flows,
 !> see the PALM homepage https://palm-model.org for further information
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  PROGRAM palm
- 
+
 
     USE arrays_3d
 
 #if defined( __parallel )
-    USE bulk_cloud_model_mod,                                                  &
+    USE bulk_cloud_model_mod,                                                                      &
         ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
 #endif
 
-    USE control_parameters,                                                    &
-        ONLY:  coupling_char, do2d_at_begin, do3d_at_begin, io_blocks,         &
-               io_group, message_string, restart_data_format_output, runnr, simulated_time_chr, spinup,   &
-               time_since_reference_point, user_interface_current_revision,    &
+    USE control_parameters,                                                                        &
+        ONLY:  coupling_char, do2d_at_begin, do3d_at_begin, io_blocks, io_group, message_string,   &
+               restart_data_format_output, runnr, simulated_time_chr, spinup,                      &
+               time_since_reference_point, user_interface_current_revision,                        &
                user_interface_required_revision, version, write_binary
 
 #if defined( __parallel )
-    USE control_parameters,                                                    &
-        ONLY:  child_domain, constant_diffusion, humidity,                     &
-               initializing_actions, neutral, passive_scalar
-#endif
-
-    USE cpulog,                                                                &
-        ONLY:  cpu_log, log_point, cpu_statistics
-
-#if defined( __parallel )
-    USE cpulog,                                                                &
+    USE control_parameters,                                                                        &
+        ONLY:  child_domain, constant_diffusion, humidity, initializing_actions, neutral,          &
+               passive_scalar
+#endif
+
+    USE cpulog,                                                                                    &
+        ONLY:  cpu_log, cpu_statistics, log_point
+
+#if defined( __parallel )
+    USE cpulog,                                                                                    &
         ONLY:  log_point_s
 #endif
 
-    USE diagnostic_output_quantities_mod,                                      &
+    USE diagnostic_output_quantities_mod,                                                          &
         ONLY:  doq_calculate
 
 #if defined( __parallel )
-    USE exchange_horiz_mod,                                                    &
+    USE exchange_horiz_mod,                                                                        &
         ONLY:  exchange_horiz
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nbgp
 #endif
@@ -129,15 +130,15 @@
     USE kinds
 
-    USE module_interface,                                                      &
-        ONLY:  module_interface_init_numerics,                                 &
-               module_interface_init_output,                                   &
+    USE module_interface,                                                                          &
+        ONLY:  module_interface_init_numerics,                                                     &
+               module_interface_init_output,                                                       &
                module_interface_last_actions
 
 
-    USE multi_agent_system_mod,                                                &
+    USE multi_agent_system_mod,                                                                    &
         ONLY:  agents_active, mas_last_actions
 
-    USE netcdf_data_input_mod,                                                 &
-        ONLY:  netcdf_data_input_inquire_file, netcdf_data_input_init,         &
+    USE netcdf_data_input_mod,                                                                     &
+        ONLY:  netcdf_data_input_inquire_file, netcdf_data_input_init,                             &
                netcdf_data_input_surface_data, netcdf_data_input_topo
 
@@ -145,19 +146,19 @@
 
 #if defined( __parallel )
-    USE pmc_particle_interface,                                                &
+    USE pmc_particle_interface,                                                                    &
         ONLY: pmcp_g_alloc_win
 
-    USE pmc_interface,                                                         &
-        ONLY:  nested_run, pmci_child_initialize, pmci_init,                   &
-               pmci_modelconfiguration, pmci_parent_initialize
-#endif
-
-    USE restart_data_mpi_io_mod,                                               &
+    USE pmc_interface,                                                                             &
+        ONLY:  nested_run, pmci_child_initialize, pmci_init, pmci_modelconfiguration,              &
+               pmci_parent_initialize
+#endif
+
+    USE restart_data_mpi_io_mod,                                                                   &
         ONLY:  rd_mpi_io_close, rd_mpi_io_open
 
-    USE surface_data_output_mod,                                               &
+    USE surface_data_output_mod,                                                                   &
         ONLY:  surface_data_output_last_action
 
-    USE write_restart_data_mod,                                                &
+    USE write_restart_data_mod,                                                                    &
         ONLY:  wrd_global, wrd_local
 
@@ -172,12 +173,13 @@
 !-- Local variables
     CHARACTER(LEN=9) ::  time_to_string  !<
-    INTEGER(iwp)     ::  i               !< loop counter for blocked I/O
+
+    INTEGER(iwp)     ::  i                   !< loop counter for blocked I/O
 #if defined( __parallel) && defined( _OPENACC )
-    INTEGER(iwp)     :: local_comm       !< local communicator (shared memory)
-    INTEGER(iwp)     :: local_num_procs  !< local number of processes
-    INTEGER(iwp)     :: local_id         !< local id
-    INTEGER(acc_device_kind) :: device_type !< device type for OpenACC
-    INTEGER(iwp)     ::  num_devices     !< number of devices visible to OpenACC
-    INTEGER(iwp)     ::  my_device       !< device used by this process
+    INTEGER(acc_device_kind) :: device_type      !< device type for OpenACC
+    INTEGER(iwp)             :: local_comm       !< local communicator (shared memory)
+    INTEGER(iwp)             :: local_num_procs  !< local number of processes
+    INTEGER(iwp)             :: local_id         !< local id
+    INTEGER(iwp)             ::  num_devices     !< number of devices visible to OpenACC
+    INTEGER(iwp)             ::  my_device       !< device used by this process
 #endif
 
@@ -187,13 +189,12 @@
 #if defined( __parallel )
 !
-!-- MPI initialisation. comm2d is preliminary set, because
-!-- it will be defined in init_pegrid but is used before in cpu_log.
+!-- MPI initialisation. comm2d is preliminary set, because it will be defined in init_pegrid but is
+!-- used before in cpu_log.
     CALL MPI_INIT( ierr )
 
 !
-!-- Initialize the coupling for nested-domain runs
-!-- comm_palm is the communicator which includes all PEs (MPI processes)
-!-- available for this (nested) model. If it is not a nested run, comm_palm
-!-- is returned as MPI_COMM_WORLD
+!-- Initialize the coupling for nested-domain runs comm_palm is the communicator which includes all
+!-- PEs (MPI processes) available for this (nested) model. If it is not a nested run, comm_palm is
+!-- returned as MPI_COMM_WORLD.
     CALL cpu_log( log_point_s(70), 'pmci_init', 'start' )
     CALL pmci_init( comm_palm )
@@ -201,7 +202,6 @@
     comm2d = comm_palm
 !
-!-- Get the (preliminary) number of MPI processes and the local PE-id (in case
-!-- of a further communicator splitting in init_coupling, these numbers will
-!-- be changed in init_pegrid).
+!-- Get the (preliminary) number of MPI processes and the local PE-id (in case of a further
+!-- communicator splitting in init_coupling, these numbers will be changed in init_pegrid).
     IF ( nested_run )  THEN
 
@@ -214,6 +214,6 @@
        CALL MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
 !
-!--    Initialize PE topology in case of coupled atmosphere-ocean runs (comm_palm
-!--    will be splitted in init_coupling)
+!--    Initialize PE topology in case of coupled atmosphere-ocean runs (comm_palm will be splitted
+!--    in init_coupling)
        CALL init_coupling
     ENDIF
@@ -221,12 +221,12 @@
 #ifdef _OPENACC
 !
-!-- Select OpenACC device to use in this process. For this find out how many
-!-- neighbors there are running on the same node and which id this process is.
+!-- Select OpenACC device to use in this process. For this find out how many neighbors there are
+!-- running on the same node and which id this process is.
     IF ( nested_run )  THEN
-       CALL MPI_COMM_SPLIT_TYPE( comm_palm, MPI_COMM_TYPE_SHARED, 0,           &
-                                 MPI_INFO_NULL, local_comm, ierr )
+       CALL MPI_COMM_SPLIT_TYPE( comm_palm, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, local_comm,    &
+                                 ierr )
     ELSE
-       CALL MPI_COMM_SPLIT_TYPE( MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, 0,      &
-                                 MPI_INFO_NULL, local_comm, ierr )
+       CALL MPI_COMM_SPLIT_TYPE( MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL,           &
+                                 local_comm, ierr )
     ENDIF
     CALL MPI_COMM_SIZE( local_comm, local_num_procs, ierr )
@@ -234,7 +234,7 @@
 
 !
-!-- This loop including the barrier is a workaround for PGI compiler versions
-!-- up to and including 18.4. Later releases are able to select their GPUs in
-!-- parallel, without running into spurious errors.
+!-- This loop including the barrier is a workaround for PGI compiler versions up to and including
+!-- 18.4. Later releases are able to select their GPUs in parallel, without running into spurious
+!-- errors.
     DO i = 0, local_num_procs-1
        CALL MPI_BARRIER( local_comm, ierr )
@@ -267,8 +267,7 @@
 
 !
-!-- Initialize dvrp logging. Also, one PE maybe split from the global
-!-- communicator for doing the dvrp output. In that case, the number of
-!-- PEs available for PALM is reduced by one and communicator comm_palm
-!-- is changed respectively.
+!-- Initialize dvrp logging. Also, one PE maybe split from the global communicator for doing the
+!-- dvrp output. In that case, the number of PEs available for PALM is reduced by one and
+!-- communicator comm_palm is changed respectively.
 #if defined( __parallel )
     CALL MPI_COMM_RANK( comm_palm, myid, ierr )
@@ -281,9 +280,8 @@
 !
 !-- Check for the user's interface version
-    IF ( user_interface_current_revision /= user_interface_required_revision )  &
-    THEN
-       message_string = 'current user-interface revision "' //                  &
-                        TRIM( user_interface_current_revision ) // '" does ' // &
-                        'not match the required revision ' //                   &
+    IF ( user_interface_current_revision /= user_interface_required_revision )  THEN
+       message_string = 'current user-interface revision "' //                                     &
+                        TRIM( user_interface_current_revision ) // '" does ' //                    &
+                        'not match the required revision ' //                                      &
                         TRIM( user_interface_required_revision )
         CALL message( 'palm', 'PA0169', 1, 2, 0, 6, 0 )
@@ -297,22 +295,21 @@
     CALL netcdf_data_input_inquire_file
 !
-!-- Read topography input data if required. This is required before the 
-!-- numerical grid is finally created in init_grid
-    CALL netcdf_data_input_topo  
-!
-!-- Generate grid parameters, initialize generic topography and further process
-!-- topography information if required
+!-- Read topography input data if required. This is required before the numerical grid is finally
+!-- created in init_grid.
+    CALL netcdf_data_input_topo
+!
+!-- Generate grid parameters, initialize generic topography and further process topography
+!-- information if required.
     CALL init_grid
 !
-!-- Initialize boundary conditions and numerics such as the multigrid solver or 
-!-- the advection routine
+!-- Initialize boundary conditions and numerics such as the multigrid solver or the advection
+!-- routine
     CALL module_interface_init_numerics
 !
-!-- Read global attributes if available.  
-    CALL netcdf_data_input_init 
-!
-!-- Read surface classification data, e.g. vegetation and soil types, water 
-!-- surfaces, etc., if available. Some of these data is required before 
-!-- check parameters is invoked.     
+!-- Read global attributes if available.
+    CALL netcdf_data_input_init
+!
+!-- Read surface classification data, e.g. vegetation and soil types, water surfaces, etc., if
+!-- available. Some of these data is required before check parameters is invoked.
     CALL netcdf_data_input_surface_data
 !
@@ -331,6 +328,5 @@
 !
 !--    Receive and interpolate initial data on children.
-!--    Child initialization must be made first if the model is both child and
-!--    parent if necessary
+!--    Child initialization must be made first if the model is both child and parent if necessary.
        IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
           CALL pmci_child_initialize
@@ -355,5 +351,5 @@
                 ENDIF
                 IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-                   CALL exchange_horiz( qr, nbgp ) 
+                   CALL exchange_horiz( qr, nbgp )
                    CALL exchange_horiz( nr, nbgp )
                 ENDIF
@@ -374,6 +370,5 @@
 
 !
-!-- Integration of the non-atmospheric equations (land surface model, urban
-!-- surface model)
+!-- Integration of the non-atmospheric equations (land surface model, urban surface model)
     IF ( spinup )  THEN
        CALL cpu_log( log_point(41), 'wall/soil spinup', 'start' )
@@ -457,5 +452,5 @@
 
        CALL cpu_log( log_point(22), 'write-restart-data', 'stop' )
-       
+
     ENDIF
 !
@@ -471,5 +466,5 @@
     CALL cpu_log( log_point(4), 'last actions', 'start' )
 
-    IF ( myid == 0 .AND. agents_active ) CALL mas_last_actions ! ToDo: move to module_interface
+    IF ( myid == 0 .AND. agents_active )  CALL mas_last_actions ! ToDo: move to module_interface
 
     CALL module_interface_last_actions
@@ -482,6 +477,5 @@
 
 !
-!-- Write run number to file (used by palmrun to create unified cycle numbers
-!-- for output files
+!-- Write run number to file (used by palmrun to create unified cycle numbers for output files).
     IF ( myid == 0  .AND.  runnr > 0 )  THEN
        OPEN( 90, FILE='RUN_NUMBER', FORM='FORMATTED' )
Index: palm/trunk/SOURCE/palm_date_time_mod.f90
===================================================================
--- palm/trunk/SOURCE/palm_date_time_mod.f90	(revision 4796)
+++ palm/trunk/SOURCE/palm_date_time_mod.f90	(revision 4797)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+! 4680 2020-09-16 10:20:34Z gronemeier
 ! Add option to fix date and time; renamed set_reference_date_time to init_date_time
 !
@@ -31,6 +34,5 @@
 ! 4227 2019-09-10 18:04:34Z gronemeier
 ! Complete rework of module date_and_time_mod:
-!  - renamed module to prevent confusion with
-!    FORTRAN Standard routine date_and_time
+!  - renamed module to prevent confusion with FORTRAN Standard routine date_and_time
 !  - introduce date_time_type
 !  - add set_reference_date_time
@@ -46,6 +48,5 @@
 ! Description:
 ! ------------
-!> This routine calculates all needed information on date and time used by
-!> other modules
+!> This routine calculates all needed information on date and time used by other modules
 !>
 !> @todo Consider leap seconds
@@ -74,34 +75,34 @@
     INTEGER(iwp), PARAMETER ::  southward_equinox  = 264_iwp                                !< Sep 21 (leap year: Sep 20)
 
+    CHARACTER(LEN=3), DIMENSION(days_per_week), PARAMETER ::                                       &
+       weekdays = (/"Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"/)                       !< names of weekdays
+
+    INTEGER, DIMENSION(months_per_year), PARAMETER ::                                              &
+       days_per_month_noleapyear = (/31,28,31,30,31,30,31,31,30,31,30,31/)                  !< days for each month (no leap year)
+
+    INTEGER, DIMENSION(months_per_year), PARAMETER ::                                              &
+       days_per_month_leapyear = (/31,29,31,30,31,30,31,31,30,31,30,31/)                    !< days for each month (leap year)
+
+    INTEGER, DIMENSION(121), PARAMETER ::  leap_years = &                                   !< list of leap years
+        (/ 1804_iwp, 1808_iwp, 1812_iwp, 1816_iwp, 1820_iwp, 1824_iwp, 1828_iwp, 1832_iwp,         &
+           1836_iwp, 1840_iwp, 1844_iwp, 1848_iwp, 1852_iwp, 1856_iwp, 1860_iwp, 1864_iwp,         &
+           1868_iwp, 1872_iwp, 1876_iwp, 1880_iwp, 1884_iwp, 1888_iwp, 1892_iwp, 1896_iwp,         &
+           1904_iwp, 1908_iwp, 1912_iwp, 1916_iwp, 1920_iwp, 1924_iwp, 1928_iwp, 1932_iwp,         &
+           1936_iwp, 1940_iwp, 1944_iwp, 1948_iwp, 1952_iwp, 1956_iwp, 1960_iwp, 1964_iwp,         &
+           1968_iwp, 1972_iwp, 1976_iwp, 1980_iwp, 1984_iwp, 1988_iwp, 1992_iwp, 1996_iwp,         &
+           2000_iwp, 2004_iwp, 2008_iwp, 2012_iwp, 2016_iwp, 2020_iwp, 2024_iwp, 2028_iwp,         &
+           2032_iwp, 2036_iwp, 2040_iwp, 2044_iwp, 2048_iwp, 2052_iwp, 2056_iwp, 2060_iwp,         &
+           2064_iwp, 2068_iwp, 2072_iwp, 2076_iwp, 2080_iwp, 2084_iwp, 2088_iwp, 2092_iwp,         &
+           2096_iwp, 2104_iwp, 2108_iwp, 2112_iwp, 2116_iwp, 2120_iwp, 2124_iwp, 2128_iwp,         &
+           2132_iwp, 2136_iwp, 2140_iwp, 2144_iwp, 2148_iwp, 2152_iwp, 2156_iwp, 2160_iwp,         &
+           2164_iwp, 2168_iwp, 2172_iwp, 2176_iwp, 2180_iwp, 2184_iwp, 2188_iwp, 2192_iwp,         &
+           2196_iwp, 2204_iwp, 2208_iwp, 2212_iwp, 2216_iwp, 2220_iwp, 2224_iwp, 2228_iwp,         &
+           2232_iwp, 2236_iwp, 2240_iwp, 2244_iwp, 2248_iwp, 2252_iwp, 2256_iwp, 2260_iwp,         &
+           2264_iwp, 2268_iwp, 2272_iwp, 2276_iwp, 2280_iwp, 2284_iwp, 2288_iwp, 2292_iwp,         &
+           2296_iwp /)
+
     REAL(wp),     PARAMETER ::  seconds_per_minute = 60.0_wp                                !< seconds in a minute
     REAL(wp),     PARAMETER ::  seconds_per_hour   = seconds_per_minute * minutes_per_hour  !< seconds in an hour
     REAL(wp),     PARAMETER ::  seconds_per_day    = seconds_per_hour * hours_per_day       !< seconds in a day
-
-    CHARACTER(LEN=3), DIMENSION(days_per_week), PARAMETER ::  &
-       weekdays = (/"Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"/)                       !< names of weekdays
-
-    INTEGER, DIMENSION(months_per_year), PARAMETER ::  &
-       days_per_month_noleapyear = (/31,28,31,30,31,30,31,31,30,31,30,31/)                  !< days for each month (no leap year)
-
-    INTEGER, DIMENSION(months_per_year), PARAMETER ::  &
-       days_per_month_leapyear = (/31,29,31,30,31,30,31,31,30,31,30,31/)                    !< days for each month (leap year)
-
-    INTEGER, DIMENSION(121), PARAMETER ::  leap_years = &                                   !< list of leap years
-        (/1804_iwp, 1808_iwp, 1812_iwp, 1816_iwp, 1820_iwp, 1824_iwp, 1828_iwp, 1832_iwp, &
-          1836_iwp, 1840_iwp, 1844_iwp, 1848_iwp, 1852_iwp, 1856_iwp, 1860_iwp, 1864_iwp, &
-          1868_iwp, 1872_iwp, 1876_iwp, 1880_iwp, 1884_iwp, 1888_iwp, 1892_iwp, 1896_iwp, &
-          1904_iwp, 1908_iwp, 1912_iwp, 1916_iwp, 1920_iwp, 1924_iwp, 1928_iwp, 1932_iwp, &
-          1936_iwp, 1940_iwp, 1944_iwp, 1948_iwp, 1952_iwp, 1956_iwp, 1960_iwp, 1964_iwp, &
-          1968_iwp, 1972_iwp, 1976_iwp, 1980_iwp, 1984_iwp, 1988_iwp, 1992_iwp, 1996_iwp, &
-          2000_iwp, 2004_iwp, 2008_iwp, 2012_iwp, 2016_iwp, 2020_iwp, 2024_iwp, 2028_iwp, &
-          2032_iwp, 2036_iwp, 2040_iwp, 2044_iwp, 2048_iwp, 2052_iwp, 2056_iwp, 2060_iwp, &
-          2064_iwp, 2068_iwp, 2072_iwp, 2076_iwp, 2080_iwp, 2084_iwp, 2088_iwp, 2092_iwp, &
-          2096_iwp, 2104_iwp, 2108_iwp, 2112_iwp, 2116_iwp, 2120_iwp, 2124_iwp, 2128_iwp, &
-          2132_iwp, 2136_iwp, 2140_iwp, 2144_iwp, 2148_iwp, 2152_iwp, 2156_iwp, 2160_iwp, &
-          2164_iwp, 2168_iwp, 2172_iwp, 2176_iwp, 2180_iwp, 2184_iwp, 2188_iwp, 2192_iwp, &
-          2196_iwp, 2204_iwp, 2208_iwp, 2212_iwp, 2216_iwp, 2220_iwp, 2224_iwp, 2228_iwp, &
-          2232_iwp, 2236_iwp, 2240_iwp, 2244_iwp, 2248_iwp, 2252_iwp, 2256_iwp, 2260_iwp, &
-          2264_iwp, 2268_iwp, 2272_iwp, 2276_iwp, 2280_iwp, 2284_iwp, 2288_iwp, 2292_iwp, &
-          2296_iwp /)
 
 !
@@ -146,20 +147,20 @@
 !
 !-- Public Interfaces
-    PUBLIC &
-       get_date_time, &
+    PUBLIC                                                                                         &
+       get_date_time,                                                                              &
        init_date_time
 !
 !-- Public variables
-    PUBLIC                 &
-       date_time_str_len,  &
-       days_per_week,      &
-       hours_per_day,      &
-       minutes_per_hour,   &
-       months_per_year,    &
-       northward_equinox,  &
-       seconds_per_minute, &
-       seconds_per_hour,   &
-       seconds_per_day,    &
-       southward_equinox,  &
+    PUBLIC                                                                                         &
+       date_time_str_len,                                                                          &
+       days_per_week,                                                                              &
+       hours_per_day,                                                                              &
+       minutes_per_hour,                                                                           &
+       months_per_year,                                                                            &
+       northward_equinox,                                                                          &
+       seconds_per_minute,                                                                         &
+       seconds_per_hour,                                                                           &
+       seconds_per_day,                                                                            &
+       southward_equinox,                                                                          &
        weekdays
 
@@ -180,4 +181,5 @@
     LOGICAL, INTENT(IN), OPTIONAL ::  use_fixed_date  !< flag to fix date
     LOGICAL, INTENT(IN), OPTIONAL ::  use_fixed_time  !< flag to fix time
+
 !
 !-- Check if date and time are already set
@@ -185,5 +187,5 @@
        !> @note This error should never be observed by a user.
        !>       It can only appear if the code was modified.
-       WRITE( message_string, * ) 'Multiple calls to init_date_time detected.&' //  &
+       WRITE( message_string, * ) 'Multiple calls to init_date_time detected.&' //                 &
                                   'This routine must not be called more than once.'
        CALL message( 'init_date_time', 'PA0680', 2, 2, 0, 6, 0 )
@@ -212,9 +214,7 @@
 !> via 'reference_date_time_str' or previously set by calling routine 'init_date_time'.
 !--------------------------------------------------------------------------------------------------!
- SUBROUTINE get_date_time( time_since_reference, reference_date_time_str,    &
-                           year, month, day, hour, minute, second, zone,     &
-                           second_of_day, second_of_year,                    &
-                           day_of_year, day_of_week, weekday, date_time_str, &
-                           days_per_month, days_per_year                     )
+ SUBROUTINE get_date_time( time_since_reference, reference_date_time_str, year, month, day, hour,  &
+                           minute, second, zone, second_of_day, second_of_year, day_of_year,       &
+                           day_of_week, weekday, date_time_str, days_per_month, days_per_year )
 
     CHARACTER(LEN=date_time_str_len), INTENT(OUT), OPTIONAL ::  date_time_str            !< date-time as string
@@ -223,13 +223,14 @@
     CHARACTER(LEN=3),                 INTENT(OUT), OPTIONAL ::  weekday                  !< weekday
 
-    INTEGER(iwp),                             INTENT(OUT), OPTIONAL ::  day              !< day of month
-    INTEGER(iwp),                             INTENT(OUT), OPTIONAL ::  day_of_week      !< number of weekday
-    INTEGER(iwp),                             INTENT(OUT), OPTIONAL ::  day_of_year      !< day of the year
-    INTEGER(iwp),                             INTENT(OUT), OPTIONAL ::  hour             !< hour of day
-    INTEGER(iwp),                             INTENT(OUT), OPTIONAL ::  minute           !< minute of hour
-    INTEGER(iwp),                             INTENT(OUT), OPTIONAL ::  month            !< month of year
-    INTEGER(iwp),                             INTENT(OUT), OPTIONAL ::  year             !< year
-    INTEGER(iwp),                             INTENT(OUT), OPTIONAL ::  zone             !< time zone
-    INTEGER(iwp),                             INTENT(OUT), OPTIONAL ::  days_per_year    !< days per year
+    INTEGER(iwp), INTENT(OUT), OPTIONAL ::  day              !< day of month
+    INTEGER(iwp), INTENT(OUT), OPTIONAL ::  day_of_week      !< number of weekday
+    INTEGER(iwp), INTENT(OUT), OPTIONAL ::  day_of_year      !< day of the year
+    INTEGER(iwp), INTENT(OUT), OPTIONAL ::  days_per_year    !< days per year
+    INTEGER(iwp), INTENT(OUT), OPTIONAL ::  hour             !< hour of day
+    INTEGER(iwp), INTENT(OUT), OPTIONAL ::  minute           !< minute of hour
+    INTEGER(iwp), INTENT(OUT), OPTIONAL ::  month            !< month of year
+    INTEGER(iwp), INTENT(OUT), OPTIONAL ::  year             !< year
+    INTEGER(iwp), INTENT(OUT), OPTIONAL ::  zone             !< time zone
+
     INTEGER(iwp), DIMENSION(months_per_year), INTENT(OUT), OPTIONAL ::  days_per_month   !< days per year
 
@@ -247,7 +248,7 @@
        !> @note This error should never be observed by a user.
        !>       It can only appear if the code was modified.
-       WRITE( message_string, * ) 'No reference date-time defined. '//                   &
-                                  'Returning date-time information is not possible. ' // &
-                                  'Either specify reference_date_time_str ' //           &
+       WRITE( message_string, * ) 'No reference date-time defined. '//                             &
+                                  'Returning date-time information is not possible. ' //           &
+                                  'Either specify reference_date_time_str ' //                     &
                                   'or set a reference via init_date_time.'
        CALL message( 'get_date_time', 'PA0677', 2, 2, 0, 6, 0 )
@@ -273,14 +274,12 @@
 !--    If date shall be fixed, revert it to the reference date if changed
        IF ( date_is_fixed )  THEN
-          IF ( date_time%year /= internal_reference_date_time%year  .OR.                          &
-               get_day_of_year( date_time ) /= get_day_of_year( internal_reference_date_time ) )  &
+          IF ( date_time%year /= internal_reference_date_time%year  .OR.                           &
+               get_day_of_year( date_time ) /= get_day_of_year( internal_reference_date_time ) )   &
           THEN
 
-             date_time%year           = internal_reference_date_time%year
-             date_time%month          = internal_reference_date_time%month
-             date_time%day            = internal_reference_date_time%day
-
-             date_time = update_leapyear_setting( date_time )
-
+             date_time%year  = internal_reference_date_time%year
+             date_time%month = internal_reference_date_time%month
+             date_time%day   = internal_reference_date_time%day
+             date_time       = update_leapyear_setting( date_time )
              date_time%second_of_year = get_second_of_year( date_time )
 
@@ -311,8 +310,8 @@
           plus_minus = '+'
        ENDIF
-       WRITE( UNIT = date_time_str,                                                 &
-              FMT = '(I4,"-",I2.2,"-",I2.2,1X,I2.2,":",I2.2,":",I2.2,1X,A1,I2.2)' ) &
-          date_time%year, date_time%month, date_time%day,                           &
-          date_time%hour, date_time%minute, INT( date_time%second ),                &
+       WRITE( UNIT = date_time_str,                                                                &
+              FMT = '(I4,"-",I2.2,"-",I2.2,1X,I2.2,":",I2.2,":",I2.2,1X,A1,I2.2)' )                &
+          date_time%year, date_time%month, date_time%day,                                          &
+          date_time%hour, date_time%minute, INT( date_time%second ),                               &
           plus_minus, ABS( date_time%zone )
     ENDIF
@@ -326,5 +325,5 @@
 !> Convert a date-time string into a date_time object.
 !--------------------------------------------------------------------------------------------------!
- FUNCTION convert_string_to_date_time( date_time_str ) RESULT( date_time )
+ FUNCTION convert_string_to_date_time( date_time_str )  RESULT( date_time )
 
     CHARACTER(LEN=date_time_str_len), INTENT(IN) ::  date_time_str  !< date-time as string
@@ -337,21 +336,27 @@
 !-- Decompose string to date-time information
     READ( UNIT = date_time_str( 1: 4), IOSTAT = read_status, FMT = '(I4)'   )  date_time%year
-    IF ( read_status == 0 )  &
+    IF ( read_status == 0 )  THEN
        READ( UNIT = date_time_str( 6: 7), IOSTAT = read_status, FMT = '(I2)'   )  date_time%month
-    IF ( read_status == 0 )  &
+    ENDIF
+    IF ( read_status == 0 )  THEN
        READ( UNIT = date_time_str( 9:10), IOSTAT = read_status, FMT = '(I2)'   )  date_time%day
-    IF ( read_status == 0 )  &
+    ENDIF
+    IF ( read_status == 0 )  THEN
        READ( UNIT = date_time_str(12:13), IOSTAT = read_status, FMT = '(I2)'   )  date_time%hour
-    IF ( read_status == 0 )  &
+    ENDIF
+    IF ( read_status == 0 )  THEN
        READ( UNIT = date_time_str(15:16), IOSTAT = read_status, FMT = '(I2)'   )  date_time%minute
-    IF ( read_status == 0 )  &
+    ENDIF
+    IF ( read_status == 0 )  THEN
        READ( UNIT = date_time_str(18:19), IOSTAT = read_status, FMT = '(F2.0)' )  date_time%second
-    IF ( read_status == 0 )  &
+    ENDIF
+    IF ( read_status == 0 )  THEN
        READ( UNIT = date_time_str(21:23), IOSTAT = read_status, FMT = '(I3)'   )  date_time%zone
+    ENDIF
 
     IF ( read_status /= 0 )  THEN
-       WRITE( message_string, * ) 'Error while converting date-time string. ' //  &
-                                  'Please check format of date-time string: "' // &
-                                  TRIM( date_time_str ) // '". ' //               &
+       WRITE( message_string, * ) 'Error while converting date-time string. ' //                   &
+                                  'Please check format of date-time string: "' //                  &
+                                  TRIM( date_time_str ) // '". ' //                                &
                                   'Format must be "YYYY-MM-DD hh:mm:ss ZZZ".'
        CALL message( 'convert_string_to_date_time', 'PA0678', 2, 2, 0, 6, 0 )
@@ -367,6 +372,5 @@
 !
 !--    Shift time to UTC and set zone to UTC
-       date_time = add_date_time( REAL( -1 * date_time%zone, KIND = wp ) &
-                                  * seconds_per_hour,                    &
+       date_time = add_date_time( REAL( -1 * date_time%zone, KIND = wp ) * seconds_per_hour,       &
                                   date_time )
        date_time%zone = 0_iwp
@@ -401,10 +405,10 @@
     date_time%second_of_year = date_time%second_of_year + inc_seconds
 !
-!-- Check if year changes
+!-- Check if year changes.
 !-- First, if year is reduced
     DO WHILE ( date_time%second_of_year < 0.0_wp )
-       date_time%year = date_time%year - 1_iwp
-       date_time = update_leapyear_setting( date_time )
-       seconds_per_year = REAL( date_time%days_per_year * seconds_per_day, KIND = wp )
+       date_time%year           = date_time%year - 1_iwp
+       date_time                = update_leapyear_setting( date_time )
+       seconds_per_year         = REAL( date_time%days_per_year * seconds_per_day, KIND = wp )
        date_time%second_of_year = date_time%second_of_year + seconds_per_year
     ENDDO
@@ -412,21 +416,21 @@
 !-- Now, if year is increased
     DO WHILE ( date_time%second_of_year > seconds_per_year )
-       date_time%year = date_time%year + 1_iwp
-       date_time = update_leapyear_setting( date_time )
+       date_time%year           = date_time%year + 1_iwp
+       date_time                = update_leapyear_setting( date_time )
        date_time%second_of_year = date_time%second_of_year - seconds_per_year
-       seconds_per_year = REAL( date_time%days_per_year * seconds_per_day, KIND = wp )
+       seconds_per_year         = REAL( date_time%days_per_year * seconds_per_day, KIND = wp )
     ENDDO
 !
 !-- Based on updated year and second_of_year, update month, day, hour, minute, second
     DO  i = 1, months_per_year
-       IF ( date_time%second_of_year < SUM( date_time%days_per_month(1:i) ) * seconds_per_day ) &
+       IF ( date_time%second_of_year < SUM( date_time%days_per_month(1:i) ) * seconds_per_day )    &
        THEN
           date_time%month  = i
-          seconds_left     = date_time%second_of_year                                &
-                           - REAL( SUM( date_time%days_per_month(1:i-1) ), KIND=wp ) &
-                           * seconds_per_day
+          seconds_left     = date_time%second_of_year                                              &
+                              - REAL( SUM( date_time%days_per_month(1:i-1) ), KIND=wp )            &
+                              * seconds_per_day
           date_time%day    = INT( seconds_left / seconds_per_day, KIND = iwp ) + 1_iwp
-          seconds_left     = seconds_left &
-                           - REAL( date_time%day - 1_iwp, KIND = wp ) * seconds_per_day
+          seconds_left     = seconds_left                                                          &
+                              - REAL( date_time%day - 1_iwp, KIND = wp ) * seconds_per_day
           date_time%hour   = INT( seconds_left / seconds_per_hour, KIND = iwp )
           seconds_left     = seconds_left - REAL( date_time%hour, KIND = wp ) * seconds_per_hour
@@ -445,5 +449,5 @@
 !> Return day of year of given date.
 !--------------------------------------------------------------------------------------------------!
- FUNCTION get_day_of_year( date_time ) RESULT( day_of_year )
+ FUNCTION get_day_of_year( date_time )  RESULT( day_of_year )
 
     INTEGER(iwp)                     ::  day_of_year         !< day of the year
@@ -466,5 +470,5 @@
 !> Return second of day of given time.
 !--------------------------------------------------------------------------------------------------!
- FUNCTION get_second_of_day( date_time ) RESULT( second_of_day )
+ FUNCTION get_second_of_day( date_time )  RESULT( second_of_day )
 
     REAL(wp)                         ::  second_of_day  !< second of the day
@@ -473,7 +477,7 @@
 
 
-    second_of_day = date_time%second                                                            &
-                  + REAL( ( date_time%hour * minutes_per_hour ) + date_time%minute, KIND = wp ) &
-                  * seconds_per_minute
+    second_of_day = date_time%second                                                               &
+                     + REAL( ( date_time%hour * minutes_per_hour ) + date_time%minute, KIND = wp ) &
+                     * seconds_per_minute
 
  END FUNCTION get_second_of_day
@@ -485,5 +489,5 @@
 !> Return second of year of given date-time.
 !--------------------------------------------------------------------------------------------------!
- FUNCTION get_second_of_year( date_time ) RESULT( second_of_year )
+ FUNCTION get_second_of_year( date_time )  RESULT( second_of_year )
 
     REAL(wp)                         ::  second_of_year  !< second of the year
@@ -492,6 +496,6 @@
 
 
-    second_of_year = get_second_of_day( date_time ) &
-                   + REAL( get_day_of_year( date_time ) - 1_iwp, KIND = wp ) * seconds_per_day
+    second_of_year = get_second_of_day( date_time )                                                &
+                     + REAL( get_day_of_year( date_time ) - 1_iwp, KIND = wp ) * seconds_per_day
 
  END FUNCTION get_second_of_year
@@ -503,9 +507,9 @@
 !> Return index of the day of the week of the given date-time.
 !--------------------------------------------------------------------------------------------------!
- FUNCTION get_day_of_week( date_time_in ) RESULT( day_of_week )
-
-    INTEGER(iwp)                     ::  date_time_internal_reference_day_of_week  !< day of week of reference date
+ FUNCTION get_day_of_week( date_time_in )  RESULT( day_of_week )
+
     INTEGER(iwp)                     ::  day_difference                            !< day between given date and reference
     INTEGER(iwp)                     ::  day_of_week                               !< day of the week
+    INTEGER(iwp)                     ::  date_time_internal_reference_day_of_week  !< day of week of reference date
 
     TYPE(date_time_type), INTENT(IN) ::  date_time_in                              !< date of which to get the weekday
@@ -548,6 +552,5 @@
     ENDIF
 !
-!-- Remove full weeks of day_difference and shift day_of_week of reference by
-!-- remaining days
+!-- Remove full weeks of day_difference and shift day_of_week of reference by remaining days.
     day_of_week = date_time_internal_reference_day_of_week + MOD( day_difference, days_per_week )
 
@@ -570,5 +573,5 @@
 !> Check if given year is a leap year and update days per month accordingly.
 !--------------------------------------------------------------------------------------------------!
- FUNCTION update_leapyear_setting( date_time_in ) RESULT( date_time_out )
+ FUNCTION update_leapyear_setting( date_time_in )  RESULT( date_time_out )
 
     TYPE(date_time_type), INTENT(IN) ::  date_time_in   !< input date-time
@@ -593,8 +596,8 @@
 ! ------------
 !> Check if given date and time are valid. Returns 0 if all checks are passed.
-!> @todo Revise error message. ATM, gives only last errorneous value even if
-!>       multiple values violate the bounds.
-!--------------------------------------------------------------------------------------------------!
- FUNCTION check_date( date_time, date_time_str ) RESULT( error_code )
+!> @todo Revise error message. ATM, gives only last errorneous value even if multiple values violate
+!>       the bounds.
+!--------------------------------------------------------------------------------------------------!
+ FUNCTION check_date( date_time, date_time_str )  RESULT( error_code )
 
 
@@ -612,9 +615,8 @@
 !
 !-- Check date
-    IF ( date_time%month < 1_iwp  .OR.     &
-         date_time%month > months_per_year )  THEN
+    IF ( date_time%month < 1_iwp  .OR.  date_time%month > months_per_year )  THEN
        error_code = 2
     ELSE
-       IF ( date_time%day < 1_iwp  .OR.                               &
+       IF ( date_time%day < 1_iwp  .OR.                                                            &
             date_time%day > date_time%days_per_month(date_time%month) )  THEN
           error_code = 3
@@ -623,14 +625,11 @@
 !
 !-- Check time
-    IF ( date_time%hour < 0_iwp  .OR.   &
-         date_time%hour > hours_per_day )  THEN
+    IF ( date_time%hour < 0_iwp  .OR.  date_time%hour > hours_per_day )  THEN
        error_code = 4
     ELSE
-        IF ( date_time%minute < 0_iwp  .OR.      &
-             date_time%minute > minutes_per_hour )  THEN
+        IF ( date_time%minute < 0_iwp  .OR.  date_time%minute > minutes_per_hour )  THEN
            error_code = 5
         ELSE
-           IF ( date_time%second < 0.0_wp  .OR.        &
-                date_time%second >= seconds_per_minute )  THEN
+           IF ( date_time%second < 0.0_wp  .OR.  date_time%second >= seconds_per_minute )  THEN
               error_code = 6
            ENDIF
@@ -638,8 +637,7 @@
     ENDIF
 !
-!-- Check time zone
+!-- Check time zone.
 !-- Bounds defined by maximum and minimum time zone present on earth
-    IF ( date_time%zone < -12_iwp  .OR.  &
-         date_time%zone > 14_iwp )  THEN
+    IF ( date_time%zone < -12_iwp  .OR.  date_time%zone > 14_iwp )  THEN
        error_code = 7
     ENDIF
@@ -647,8 +645,8 @@
 !-- Raise error if any check is marked invalid
     IF ( error_code /= 0 )  THEN
-       WRITE( message_string, * ) 'Date-time values out of bounds: "' //                    &
-                                  TRIM( error_str_list(error_code) ) //                     &
-                                  '" is out of bounds. Please check date-time string: "' // &
-                                  TRIM( date_time_str ) // '". ' //                         &
+       WRITE( message_string, * ) 'Date-time values out of bounds: "' //                           &
+                                  TRIM( error_str_list(error_code) ) //                            &
+                                  '" is out of bounds. Please check date-time string: "' //        &
+                                  TRIM( date_time_str ) // '". ' //                                &
                                   'Format must be "YYYY-MM-DD hh:mm:ss ZZZ".'
        CALL message( 'check_date', 'PA0679', 2, 2, 0, 6, 0 )
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 4796)
+++ palm/trunk/SOURCE/parin.f90	(revision 4797)
@@ -1,20 +1,19 @@
 !> @file parin.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
+!
+! 4783 2020-11-13 13:58:45Z raasch
 ! bugfix for reading restart data with MPI-I/O (does not work with blockwise I/O)
 !
@@ -76,6 +78,5 @@
 !
 ! 4022 2019-06-12 11:52:39Z suehring
-! Change default top boundary condition for pressure to Neumann in offline
-! nesting case
+! Change default top boundary condition for pressure to Neumann in offline nesting case
 !
 ! 4017 2019-06-06 12:16:46Z schwenkel
@@ -83,6 +84,6 @@
 !
 ! 3885 2019-04-11 11:29:34Z kanani
-! Changes related to global restructuring of location messages and introduction
-! of additional debug messages
+! Changes related to global restructuring of location messages and introduction of additional debug
+! messages
 !
 ! 3806 2019-03-21 12:45:50Z raasch
@@ -107,11 +108,10 @@
 !>
 !> @todo: Revise max_pr_cs (profiles for chemistry)
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE parin
 
 
-    USE arrays_3d,                                                             &
-        ONLY:  pt_init, q_init, ref_state, s_init, sa_init,                    &
-               ug, u_init, v_init, vg
+    USE arrays_3d,                                                                                 &
+        ONLY:  pt_init, q_init, ref_state, s_init, sa_init, ug, u_init, v_init, vg
 
     USE chem_modules
@@ -119,42 +119,42 @@
     USE control_parameters
 
-    USE cpulog,                                                                &
+    USE cpulog,                                                                                    &
         ONLY:  cpu_log_barrierwait
 
-    USE grid_variables,                                                        &
+    USE grid_variables,                                                                            &
         ONLY:  dx, dy
 
-    USE indices,                                                               &
+    USE indices,                                                                                   &
         ONLY:  nx, ny, nz
 
     USE kinds
 
-    USE model_1d_mod,                                                          &
+    USE model_1d_mod,                                                                              &
         ONLY:  damp_level_1d, dt_pr_1d, dt_run_control_1d, end_time_1d
 
-    USE module_interface,                                                      &
+    USE module_interface,                                                                          &
         ONLY:  module_interface_parin
 
-    USE netcdf_interface,                                                      &
+    USE netcdf_interface,                                                                          &
         ONLY:  netcdf_data_format, netcdf_deflate, netcdf_precision
 
     USE pegrid
 
-    USE pmc_interface,                                                         &
+    USE pmc_interface,                                                                             &
         ONLY:  nested_run, nesting_mode
 
-    USE profil_parameter,                                                      &
+    USE profil_parameter,                                                                          &
         ONLY:  cross_profiles, profile_columns, profile_rows
 
-    USE progress_bar,                                                          &
+    USE progress_bar,                                                                              &
         ONLY :  progress_bar_disabled
 
-    USE read_restart_data_mod,                                                 &
+    USE read_restart_data_mod,                                                                     &
         ONLY:  rrd_global
 
-    USE statistics,                                                            &
+    USE statistics,                                                                                &
         ONLY:  hom, hom_sum, pr_palm, statistic_regions
 
-    USE turbulence_closure_mod,                                                &
+    USE turbulence_closure_mod,                                                                    &
         ONLY:  rans_const_c, rans_const_sigma
 
@@ -169,228 +169,553 @@
     INTEGER(iwp) ::  ioerr          !< error flag for open/read/write
 
-    NAMELIST /inipar/  alpha_surface, approximation, bc_e_b,     &
-                       bc_lr, bc_ns, bc_p_b, bc_p_t, bc_pt_b, bc_pt_t, bc_q_b, &
-             bc_q_t,bc_s_b, bc_s_t, bc_uv_b, bc_uv_t,                 &
-             building_height, building_length_x,          &
-             building_length_y, building_wall_left, building_wall_south,       &
-             calc_soil_moisture_during_spinup,                                 &
-             call_psolver_at_all_substeps,  &
-             canyon_height,                                                    &
-             canyon_width_x, canyon_width_y, canyon_wall_left,                 &
-             canyon_wall_south, cfl_factor, check_realistic_q, cloud_droplets, &
-             collective_wait, complex_terrain,           &
-             conserve_volume_flow,                                             &
-             conserve_volume_flow_mode, constant_flux_layer,                   &
-             coupling_start_time,             &
-             cycle_mg, damp_level_1d,                                          &
-             data_output_during_spinup,                                        &
-             origin_date_time,                                                 &
-             dissipation_1d,                                                   &
-             dp_external, dp_level_b, dp_smooth, dpdxy,    &
-             dt, dt_pr_1d, dt_run_control_1d, dt_spinup, dx, dy, dz, dz_max,   &
-             dz_stretch_factor, dz_stretch_level, dz_stretch_level_start,      &
-             dz_stretch_level_end, end_time_1d, ensemble_member_nr, e_init,    &
-             e_min, fft_method, flux_input_mode, flux_output_mode,             &
-             galilei_transformation, humidity,                                 &
-             inflow_damping_height, inflow_damping_width,                      &
-             inflow_disturbance_begin, inflow_disturbance_end,                 &
-             initializing_actions, km_constant,                                &
-             large_scale_forcing, large_scale_subsidence, latitude,            &
-             longitude,                                 &
-             loop_optimization, lsf_exception, masking_method, mg_cycles,      &
-             mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             monotonic_limiter_z,                                              &
-             netcdf_precision, neutral, ngsrb,                                 &
-             nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
-             omega_sor, outflow_source_plane, passive_scalar,                  &
-             prandtl_number, psolver, pt_damping_factor,  pt_damping_width,    &
-             pt_reference, pt_surface, pt_surface_heating_rate,                &
-             pt_surface_initial_change, pt_vertical_gradient,                  &
-             pt_vertical_gradient_level, q_surface, q_surface_initial_change,  &
-             q_vertical_gradient, q_vertical_gradient_level,                   &
-             random_generator, random_heatflux, rans_const_c, rans_const_sigma,&
-             rayleigh_damping_factor, rayleigh_damping_height,                 &
-             recycling_method_for_thermodynamic_quantities, recycling_width,   &
-             reference_state, residual_limit,                                  &
-             rotation_angle,                                                   &
-             roughness_length,                                                 &
-             scalar_advec,   &
-             scalar_rayleigh_damping,                              &
-             spinup_time, spinup_pt_amplitude, spinup_pt_mean,                 &
-             statistic_regions, subs_vertical_gradient,                        &
-             subs_vertical_gradient_level, surface_heatflux, surface_pressure, &
-             surface_scalarflux, surface_waterflux,                            &
-             s_surface, s_surface_initial_change, s_vertical_gradient,         &
-             s_vertical_gradient_level, timestep_scheme,                       &
-             topography, topography_grid_convention, top_heatflux,             &
-             top_momentumflux_u, top_momentumflux_v,                           &
-             top_scalarflux, transpose_compute_overlap,                        &
-             tunnel_height, tunnel_length, tunnel_width_x, tunnel_width_y,     &
-             tunnel_wall_depth, turbulence_closure,                            &
-             turbulent_inflow, turbulent_outflow,                              &
-             use_subsidence_tendencies, ug_surface, ug_vertical_gradient,      &
-             use_fixed_date, use_fixed_time,                                   &
-             use_free_convection_scaling,                                      &
-             ug_vertical_gradient_level, use_surface_fluxes, use_cmax,         &
-             use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke,      &
-             uv_heights, u_bulk, u_profile, vdi_checks, vg_surface,            &
-             vg_vertical_gradient,  &
-             vg_vertical_gradient_level, v_bulk, v_profile,&
-             wall_adjustment, wall_heatflux, wall_humidityflux,                &
-             wall_scalarflux, y_shift, zeta_max, zeta_min,  &
-             z0h_factor
-
-    NAMELIST /initialization_parameters/  alpha_surface,                       &
-             approximation, bc_e_b,                                            &
-             bc_lr, bc_ns, bc_p_b, bc_p_t, bc_pt_b, bc_pt_t, bc_q_b,           &
-             bc_q_t,bc_s_b, bc_s_t, bc_uv_b, bc_uv_t,                          &
-             building_height, building_length_x,                               &
-             building_length_y, building_wall_left, building_wall_south,       &
-             calc_soil_moisture_during_spinup,                                 &
-             call_psolver_at_all_substeps,                                     &
-             canyon_height,                                                    &
-             canyon_width_x, canyon_width_y, canyon_wall_left,                 &
-             canyon_wall_south, cfl_factor, check_realistic_q, cloud_droplets, &
-             collective_wait, complex_terrain,                                 &
-             conserve_volume_flow,                                             &
-             conserve_volume_flow_mode, constant_flux_layer,                   &
-             coupling_start_time,                                              &
-             cycle_mg, damp_level_1d,                                          &
-             data_output_during_spinup,                                        &
-             origin_date_time,                                                 &
-             dissipation_1d,                                                   &
-             dp_external, dp_level_b, dp_smooth, dpdxy,                        &
-             dt, dt_pr_1d, dt_run_control_1d, dt_spinup, dx, dy, dz, dz_max,   &
-             dz_stretch_factor, dz_stretch_level, dz_stretch_level_start,      &
-             dz_stretch_level_end, end_time_1d, ensemble_member_nr, e_init,    &
-             e_min, fft_method, flux_input_mode, flux_output_mode,             &
-             galilei_transformation, humidity,                                 &
-             inflow_damping_height, inflow_damping_width,                      &
-             inflow_disturbance_begin, inflow_disturbance_end,                 &
-             initializing_actions, km_constant,                                &
-             large_scale_forcing, large_scale_subsidence, latitude,            &
-             longitude,                                                        &
-             loop_optimization, lsf_exception, masking_method, mg_cycles,      &
-             mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             monotonic_limiter_z,                                              &
-             netcdf_precision, neutral, ngsrb,                                 &
-             nsor, nsor_ini, nudging, nx, ny, nz, ocean_mode, omega,           &
-             omega_sor, outflow_source_plane, passive_scalar,                  &
-             prandtl_number, psolver, pt_damping_factor, pt_damping_width,     &
-             pt_surface_heating_rate, pt_reference, pt_surface,                &
-             pt_surface_initial_change, pt_vertical_gradient,                  &
-             pt_vertical_gradient_level, q_surface, q_surface_initial_change,  &
-             q_vertical_gradient, q_vertical_gradient_level,                   &
-             random_generator, random_heatflux, rans_const_c, rans_const_sigma,&
-             rayleigh_damping_factor, rayleigh_damping_height,                 &
-             recycling_method_for_thermodynamic_quantities, recycling_width,   &
-             reference_state, residual_limit,                                  &
-             restart_data_format, restart_data_format_input, restart_data_format_output,           &
-             rotation_angle,                                                   &
-             roughness_length, scalar_advec,                                   &
-             scalar_rayleigh_damping,                                          &
-             spinup_time, spinup_pt_amplitude, spinup_pt_mean,                 &
-             statistic_regions, subs_vertical_gradient,                        &
-             subs_vertical_gradient_level, surface_heatflux, surface_pressure, &
-             surface_scalarflux, surface_waterflux,                            &
-             s_surface, s_surface_initial_change, s_vertical_gradient,         &
-             s_vertical_gradient_level, timestep_scheme,                       &
-             topography, topography_grid_convention, top_heatflux,             &
-             top_momentumflux_u, top_momentumflux_v,                           &
-             top_scalarflux, transpose_compute_overlap,                        &
-             tunnel_height, tunnel_length, tunnel_width_x, tunnel_width_y,     &
-             tunnel_wall_depth, turbulence_closure,                            &
-             turbulent_inflow, turbulent_outflow,                              &
-             use_subsidence_tendencies, ug_surface, ug_vertical_gradient,      &
-             ug_vertical_gradient_level, use_surface_fluxes, use_cmax,         &
-             use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke,      &
-             use_fixed_date, use_fixed_time,                                   &
-             use_free_convection_scaling,                                      &
-             uv_heights, u_bulk, u_profile, vdi_checks,                        &
-             vg_surface, vg_vertical_gradient,  &
-             vg_vertical_gradient_level, v_bulk, v_profile,                    &
-             wall_adjustment, wall_heatflux, wall_humidityflux,                &
-             wall_scalarflux, y_shift, zeta_max, zeta_min, z0h_factor
-
-    NAMELIST /d3par/  averaging_interval, averaging_interval_pr,               &
-             cpu_log_barrierwait, create_disturbances,                         &
-             cross_profiles, data_output, data_output_masks,                   &
-             data_output_pr, data_output_2d_on_each_pe,                        &
-             debug_output,                                                     &
-             debug_output_timestep,                                            &
-             disturbance_amplitude,                                            &
-             disturbance_energy_limit, disturbance_level_b,                    &
-             disturbance_level_t, do2d_at_begin, do3d_at_begin,                &
-             dt, dt_averaging_input, dt_averaging_input_pr,                    &
-             dt_coupling, dt_data_output, dt_data_output_av, dt_disturb,       &
-             dt_domask, dt_dopr, dt_dopr_listing, dt_dots, dt_do2d_xy,         &
-             dt_do2d_xz, dt_do2d_yz, dt_do3d, dt_max, dt_restart,              &
-             dt_run_control,end_time, force_print_header, mask_k_over_surface, &
-             mask_scale_x,                                                     &
-             mask_scale_y, mask_scale_z, mask_x, mask_y, mask_z, mask_x_loop,  &
-             mask_y_loop, mask_z_loop, netcdf_data_format, netcdf_deflate,     &
-             normalizing_region, npex, npey, nz_do3d,                          &
-             profile_columns, profile_rows,     &
-             restart_time, section_xy, section_xz, section_yz,                 &
-             skip_time_data_output, skip_time_data_output_av, skip_time_dopr,  &
-             skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,          &
-             skip_time_do3d, skip_time_domask, synchronous_exchange,           &
-             termination_time_needed
-
-    NAMELIST /runtime_parameters/  averaging_interval, averaging_interval_pr,  &
-             cpu_log_barrierwait, create_disturbances,                         &
-             cross_profiles, data_output, data_output_masks,                   &
-             data_output_pr, data_output_2d_on_each_pe,                        &
-             debug_output,                                                     &
-             debug_output_timestep,                                            &
-             disturbance_amplitude,                                            &
-             disturbance_energy_limit, disturbance_level_b,                    &
-             disturbance_level_t, do2d_at_begin, do3d_at_begin,                &
-             dt, dt_averaging_input, dt_averaging_input_pr,                    &
-             dt_coupling, dt_data_output, dt_data_output_av, dt_disturb,       &
-             dt_domask, dt_dopr, dt_dopr_listing, dt_dots, dt_do2d_xy,         &
-             dt_do2d_xz, dt_do2d_yz, dt_do3d, dt_max, dt_restart,              &
-             dt_run_control,end_time, force_print_header, mask_k_over_surface, &
-             mask_scale_x,                                                     &
-             mask_scale_y, mask_scale_z, mask_x, mask_y, mask_z, mask_x_loop,  &
-             mask_y_loop, mask_z_loop, netcdf_data_format, netcdf_deflate,     &
-             normalizing_region, npex, npey, nz_do3d,                          &
-             profile_columns, profile_rows,     &
-             restart_time, section_xy, section_xz, section_yz,                 &
-             restart_data_format, restart_data_format_input, restart_data_format_output,           &
-             skip_time_data_output, skip_time_data_output_av, skip_time_dopr,  &
-             skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,          &
-             skip_time_do3d, skip_time_domask, synchronous_exchange,           &
-             termination_time_needed
-
-    NAMELIST /envpar/  progress_bar_disabled, host,                            &
-                       maximum_cpu_time_allowed, maximum_parallel_io_streams,  &
-                       read_svf, revision, run_identifier, tasks_per_node,     &
-                       write_binary, write_svf
-
-!
-!-- First read values of environment variables (this NAMELIST file is
-!-- generated by palmrun)
+    NAMELIST /inipar/  alpha_surface,                                                              &
+                       approximation,                                                              &
+                       bc_e_b,                                                                     &
+                       bc_lr,                                                                      &
+                       bc_ns,                                                                      &
+                       bc_p_b,                                                                     &
+                       bc_p_t,                                                                     &
+                       bc_pt_b,                                                                    &
+                       bc_pt_t,                                                                    &
+                       bc_q_b,                                                                     &
+                       bc_q_t,                                                                     &
+                       bc_s_b,                                                                     &
+                       bc_s_t,                                                                     &
+                       bc_uv_b,                                                                    &
+                       bc_uv_t,                                                                    &
+                       building_height,                                                            &
+                       building_length_x,                                                          &
+                       building_length_y,                                                          &
+                       building_wall_left,                                                         &
+                       building_wall_south,                                                        &
+                       calc_soil_moisture_during_spinup,                                           &
+                       call_psolver_at_all_substeps,                                               &
+                       canyon_height,                                                              &
+                       canyon_wall_left,                                                           &
+                       canyon_wall_south,                                                          &
+                       canyon_width_x,                                                             &
+                       canyon_width_y,                                                             &
+                       cfl_factor,                                                                 &
+                       check_realistic_q,                                                          &
+                       cloud_droplets,                                                             &
+                       collective_wait,                                                            &
+                       complex_terrain,                                                            &
+                       conserve_volume_flow,                                                       &
+                       conserve_volume_flow_mode,                                                  &
+                       constant_flux_layer,                                                        &
+                       coupling_start_time,                                                        &
+                       cycle_mg,                                                                   &
+                       damp_level_1d,                                                              &
+                       data_output_during_spinup,                                                  &
+                       dissipation_1d,                                                             &
+                       dp_external,                                                                &
+                       dp_level_b,                                                                 &
+                       dp_smooth,                                                                  &
+                       dpdxy,                                                                      &
+                       dt,                                                                         &
+                       dt_pr_1d,                                                                   &
+                       dt_run_control_1d,                                                          &
+                       dt_spinup,                                                                  &
+                       dx,                                                                         &
+                       dy,                                                                         &
+                       dz,                                                                         &
+                       dz_max,                                                                     &
+                       dz_stretch_factor,                                                          &
+                       dz_stretch_level,                                                           &
+                       dz_stretch_level_end,                                                       &
+                       dz_stretch_level_start,                                                     &
+                       e_init,                                                                     &
+                       e_min,                                                                      &
+                       end_time_1d,                                                                &
+                       ensemble_member_nr,                                                         &
+                       fft_method,                                                                 &
+                       flux_input_mode,                                                            &
+                       flux_output_mode,                                                           &
+                       galilei_transformation,                                                     &
+                       humidity,                                                                   &
+                       inflow_damping_height,                                                      &
+                       inflow_damping_width,                                                       &
+                       inflow_disturbance_begin,                                                   &
+                       inflow_disturbance_end,                                                     &
+                       initializing_actions,                                                       &
+                       km_constant,                                                                &
+                       large_scale_forcing,                                                        &
+                       large_scale_subsidence,                                                     &
+                       latitude,                                                                   &
+                       longitude,                                                                  &
+                       loop_optimization,                                                          &
+                       lsf_exception,                                                              &
+                       masking_method,                                                             &
+                       mg_cycles,                                                                  &
+                       mg_switch_to_pe0_level,                                                     &
+                       mixing_length_1d,                                                           &
+                       momentum_advec,                                                             &
+                       monotonic_limiter_z,                                                        &
+                       netcdf_precision,                                                           &
+                       neutral,                                                                    &
+                       ngsrb,                                                                      &
+                       nsor,                                                                       &
+                       nsor_ini,                                                                   &
+                       nudging,                                                                    &
+                       nx,                                                                         &
+                       ny,                                                                         &
+                       nz,                                                                         &
+                       ocean_mode,                                                                 &
+                       omega,                                                                      &
+                       omega_sor,                                                                  &
+                       origin_date_time,                                                           &
+                       outflow_source_plane,                                                       &
+                       passive_scalar,                                                             &
+                       prandtl_number,                                                             &
+                       psolver,                                                                    &
+                       pt_damping_factor,                                                          &
+                       pt_damping_width,                                                           &
+                       pt_reference,                                                               &
+                       pt_surface,                                                                 &
+                       pt_surface_heating_rate,                                                    &
+                       pt_surface_initial_change,                                                  &
+                       pt_vertical_gradient,                                                       &
+                       pt_vertical_gradient_level,                                                 &
+                       q_surface,                                                                  &
+                       q_surface_initial_change,                                                   &
+                       q_vertical_gradient,                                                        &
+                       q_vertical_gradient_level,                                                  &
+                       random_generator,                                                           &
+                       random_heatflux,                                                            &
+                       rans_const_c,                                                               &
+                       rans_const_sigma,                                                           &
+                       rayleigh_damping_factor,                                                    &
+                       rayleigh_damping_height,                                                    &
+                       recycling_method_for_thermodynamic_quantities,                              &
+                       recycling_width,                                                            &
+                       reference_state,                                                            &
+                       residual_limit,                                                             &
+                       rotation_angle,                                                             &
+                       roughness_length,                                                           &
+                       scalar_advec,                                                               &
+                       scalar_rayleigh_damping,                                                    &
+                       spinup_time,                                                                &
+                       spinup_pt_amplitude,                                                        &
+                       spinup_pt_mean,                                                             &
+                       statistic_regions,                                                          &
+                       subs_vertical_gradient,                                                     &
+                       subs_vertical_gradient_level,                                               &
+                       surface_heatflux,                                                           &
+                       surface_pressure,                                                           &
+                       surface_scalarflux,                                                         &
+                       surface_waterflux,                                                          &
+                       s_surface,                                                                  &
+                       s_surface_initial_change,                                                   &
+                       s_vertical_gradient,                                                        &
+                       s_vertical_gradient_level,                                                  &
+                       timestep_scheme,                                                            &
+                       topography,                                                                 &
+                       topography_grid_convention,                                                 &
+                       top_heatflux,                                                               &
+                       top_momentumflux_u,                                                         &
+                       top_momentumflux_v,                                                         &
+                       top_scalarflux,                                                             &
+                       transpose_compute_overlap,                                                  &
+                       tunnel_height,                                                              &
+                       tunnel_length,                                                              &
+                       tunnel_wall_depth,                                                          &
+                       tunnel_width_x,                                                             &
+                       tunnel_width_y,                                                             &
+                       turbulence_closure,                                                         &
+                       turbulent_inflow,                                                           &
+                       turbulent_outflow,                                                          &
+                       u_bulk,                                                                     &
+                       u_profile,                                                                  &
+                       ug_surface,                                                                 &
+                       ug_vertical_gradient,                                                       &
+                       ug_vertical_gradient_level,                                                 &
+                       use_cmax,                                                                   &
+                       use_fixed_date,                                                             &
+                       use_fixed_time,                                                             &
+                       use_free_convection_scaling,                                                &
+                       use_ug_for_galilei_tr,                                                      &
+                       use_subsidence_tendencies,                                                  &
+                       use_surface_fluxes,                                                         &
+                       use_top_fluxes,                                                             &
+                       use_upstream_for_tke,                                                       &
+                       uv_heights,                                                                 &
+                       v_bulk,                                                                     &
+                       v_profile,                                                                  &
+                       vdi_checks,                                                                 &
+                       vg_surface,                                                                 &
+                       vg_vertical_gradient,                                                       &
+                       vg_vertical_gradient_level,                                                 &
+                       wall_adjustment,                                                            &
+                       wall_heatflux,                                                              &
+                       wall_humidityflux,                                                          &
+                       wall_scalarflux,                                                            &
+                       y_shift,                                                                    &
+                       zeta_max,                                                                   &
+                       zeta_min,                                                                   &
+                       z0h_factor
+
+    NAMELIST /initialization_parameters/  alpha_surface,                                           &
+                                          approximation,                                           &
+                                          bc_e_b,                                                  &
+                                          bc_lr,                                                   &
+                                          bc_ns,                                                   &
+                                          bc_p_b,                                                  &
+                                          bc_p_t,                                                  &
+                                          bc_pt_b,                                                 &
+                                          bc_pt_t,                                                 &
+                                          bc_q_b,                                                  &
+                                          bc_q_t,                                                  &
+                                          bc_s_b,                                                  &
+                                          bc_s_t,                                                  &
+                                          bc_uv_b,                                                 &
+                                          bc_uv_t,                                                 &
+                                          building_height,                                         &
+                                          building_length_x,                                       &
+                                          building_length_y,                                       &
+                                          building_wall_left,                                      &
+                                          building_wall_south,                                     &
+                                          calc_soil_moisture_during_spinup,                        &
+                                          call_psolver_at_all_substeps,                            &
+                                          canyon_height,                                           &
+                                          canyon_wall_left,                                        &
+                                          canyon_wall_south,                                       &
+                                          canyon_width_x,                                          &
+                                          canyon_width_y,                                          &
+                                          cfl_factor,                                              &
+                                          check_realistic_q,                                       &
+                                          cloud_droplets,                                          &
+                                          collective_wait,                                         &
+                                          complex_terrain,                                         &
+                                          conserve_volume_flow,                                    &
+                                          conserve_volume_flow_mode,                               &
+                                          constant_flux_layer,                                     &
+                                          coupling_start_time,                                     &
+                                          cycle_mg,                                                &
+                                          damp_level_1d,                                           &
+                                          data_output_during_spinup,                               &
+                                          dissipation_1d,                                          &
+                                          dp_external,                                             &
+                                          dp_level_b,                                              &
+                                          dp_smooth, dpdxy,                                        &
+                                          dt,                                                      &
+                                          dt_pr_1d,                                                &
+                                          dt_run_control_1d,                                       &
+                                          dt_spinup,                                               &
+                                          dx,                                                      &
+                                          dy,                                                      &
+                                          dz,                                                      &
+                                          dz_max,                                                  &
+                                          dz_stretch_factor,                                       &
+                                          dz_stretch_level,                                        &
+                                          dz_stretch_level_start,                                  &
+                                          dz_stretch_level_end,                                    &
+                                          e_init,                                                  &
+                                          e_min,                                                   &
+                                          end_time_1d,                                             &
+                                          ensemble_member_nr,                                      &
+                                          fft_method,                                              &
+                                          flux_input_mode,                                         &
+                                          flux_output_mode,                                        &
+                                          galilei_transformation,                                  &
+                                          humidity,                                                &
+                                          inflow_damping_height,                                   &
+                                          inflow_damping_width,                                    &
+                                          inflow_disturbance_begin,                                &
+                                          inflow_disturbance_end,                                  &
+                                          initializing_actions,                                    &
+                                          km_constant,                                             &
+                                          large_scale_forcing,                                     &
+                                          large_scale_subsidence,                                  &
+                                          latitude,                                                &
+                                          longitude,                                               &
+                                          loop_optimization,                                       &
+                                          lsf_exception,                                           &
+                                          masking_method,                                          &
+                                          mg_cycles,                                               &
+                                          mg_switch_to_pe0_level,                                  &
+                                          mixing_length_1d,                                        &
+                                          momentum_advec,                                          &
+                                          monotonic_limiter_z,                                     &
+                                          netcdf_precision,                                        &
+                                          neutral,                                                 &
+                                          ngsrb,                                                   &
+                                          nsor,                                                    &
+                                          nsor_ini,                                                &
+                                          nudging,                                                 &
+                                          nx,                                                      &
+                                          ny,                                                      &
+                                          nz,                                                      &
+                                          ocean_mode,                                              &
+                                          omega,                                                   &
+                                          omega_sor,                                               &
+                                          origin_date_time,                                        &
+                                          outflow_source_plane,                                    &
+                                          passive_scalar,                                          &
+                                          prandtl_number,                                          &
+                                          psolver,                                                 &
+                                          pt_damping_factor,                                       &
+                                          pt_damping_width,                                        &
+                                          pt_reference, pt_surface,                                &
+                                          pt_surface_heating_rate,                                 &
+                                          pt_surface_initial_change,                               &
+                                          pt_vertical_gradient,                                    &
+                                          pt_vertical_gradient_level,                              &
+                                          q_surface,                                               &
+                                          q_surface_initial_change,                                &
+                                          q_vertical_gradient,                                     &
+                                          q_vertical_gradient_level,                               &
+                                          random_generator,                                        &
+                                          random_heatflux,                                         &
+                                          rans_const_c,                                            &
+                                          rans_const_sigma,                                        &
+                                          rayleigh_damping_factor,                                 &
+                                          rayleigh_damping_height,                                 &
+                                          recycling_method_for_thermodynamic_quantities,           &
+                                          recycling_width,                                         &
+                                          reference_state,                                         &
+                                          residual_limit,                                          &
+                                          restart_data_format,                                     &
+                                          restart_data_format_input,                               &
+                                          restart_data_format_output,                              &
+                                          rotation_angle,                                          &
+                                          roughness_length,                                        &
+                                          s_surface,                                               &
+                                          s_surface_initial_change,                                &
+                                          s_vertical_gradient,                                     &
+                                          s_vertical_gradient_level,                               &
+                                          scalar_advec,                                            &
+                                          scalar_rayleigh_damping,                                 &
+                                          spinup_time,                                             &
+                                          spinup_pt_amplitude,                                     &
+                                          spinup_pt_mean,                                          &
+                                          statistic_regions,                                       &
+                                          subs_vertical_gradient,                                  &
+                                          subs_vertical_gradient_level,                            &
+                                          surface_heatflux,                                        &
+                                          surface_pressure,                                        &
+                                          surface_scalarflux,                                      &
+                                          surface_waterflux,                                       &
+                                          timestep_scheme,                                         &
+                                          topography,                                              &
+                                          topography_grid_convention,                              &
+                                          top_heatflux,                                            &
+                                          top_momentumflux_u,                                      &
+                                          top_momentumflux_v,                                      &
+                                          top_scalarflux,                                          &
+                                          transpose_compute_overlap,                               &
+                                          tunnel_height,                                           &
+                                          tunnel_length,                                           &
+                                          tunnel_wall_depth,                                       &
+                                          tunnel_width_x,                                          &
+                                          tunnel_width_y,                                          &
+                                          turbulence_closure,                                      &
+                                          turbulent_inflow,                                        &
+                                          turbulent_outflow,                                       &
+                                          u_bulk,                                                  &
+                                          u_profile,                                               &
+                                          ug_surface,                                              &
+                                          ug_vertical_gradient,                                    &
+                                          ug_vertical_gradient_level,                              &
+                                          use_fixed_date,                                          &
+                                          use_fixed_time,                                          &
+                                          use_free_convection_scaling,                             &
+                                          use_ug_for_galilei_tr,                                   &
+                                          use_subsidence_tendencies,                               &
+                                          use_surface_fluxes, use_cmax,                            &
+                                          use_top_fluxes,                                          &
+                                          use_upstream_for_tke,                                    &
+                                          uv_heights,                                              &
+                                          v_bulk,                                                  &
+                                          v_profile,                                               &
+                                          vdi_checks,                                              &
+                                          vg_surface,                                              &
+                                          vg_vertical_gradient,                                    &
+                                          vg_vertical_gradient_level,                              &
+                                          wall_adjustment,                                         &
+                                          wall_heatflux,                                           &
+                                          wall_humidityflux,                                       &
+                                          wall_scalarflux,                                         &
+                                          y_shift,                                                 &
+                                          zeta_max,                                                &
+                                          zeta_min,                                                &
+                                          z0h_factor
+
+    NAMELIST /d3par/  averaging_interval,                                                          &
+                      averaging_interval_pr,                                                       &
+                      cpu_log_barrierwait,                                                         &
+                      create_disturbances,                                                         &
+                      cross_profiles,                                                              &
+                      data_output,                                                                 &
+                      data_output_2d_on_each_pe,                                                   &
+                      data_output_masks,                                                           &
+                      data_output_pr,                                                              &
+                      debug_output,                                                                &
+                      debug_output_timestep,                                                       &
+                      disturbance_amplitude,                                                       &
+                      disturbance_energy_limit,                                                    &
+                      disturbance_level_b,                                                         &
+                      disturbance_level_t,                                                         &
+                      do2d_at_begin,                                                               &
+                      do3d_at_begin,                                                               &
+                      dt,                                                                          &
+                      dt_averaging_input,                                                          &
+                      dt_averaging_input_pr,                                                       &
+                      dt_coupling,                                                                 &
+                      dt_data_output,                                                              &
+                      dt_data_output_av,                                                           &
+                      dt_disturb,                                                                  &
+                      dt_domask,                                                                   &
+                      dt_dopr,                                                                     &
+                      dt_dopr_listing,                                                             &
+                      dt_dots,                                                                     &
+                      dt_do2d_xy,                                                                  &
+                      dt_do2d_xz,                                                                  &
+                      dt_do2d_yz,                                                                  &
+                      dt_do3d,                                                                     &
+                      dt_max,                                                                      &
+                      dt_restart,                                                                  &
+                      dt_run_control,                                                              &
+                      end_time,                                                                    &
+                      force_print_header,                                                          &
+                      mask_k_over_surface,                                                         &
+                      mask_scale_x,                                                                &
+                      mask_scale_y,                                                                &
+                      mask_scale_z,                                                                &
+                      mask_x,                                                                      &
+                      mask_y,                                                                      &
+                      mask_z,                                                                      &
+                      mask_x_loop,                                                                 &
+                      mask_y_loop,                                                                 &
+                      mask_z_loop,                                                                 &
+                      netcdf_data_format,                                                          &
+                      netcdf_deflate,                                                              &
+                      normalizing_region,                                                          &
+                      npex,                                                                        &
+                      npey,                                                                        &
+                      nz_do3d,                                                                     &
+                      profile_columns,                                                             &
+                      profile_rows,                                                                &
+                      restart_time,                                                                &
+                      section_xy,                                                                  &
+                      section_xz,                                                                  &
+                      section_yz,                                                                  &
+                      skip_time_data_output,                                                       &
+                      skip_time_data_output_av,                                                    &
+                      skip_time_dopr,                                                              &
+                      skip_time_do2d_xy,                                                           &
+                      skip_time_do2d_xz,                                                           &
+                      skip_time_do2d_yz,                                                           &
+                      skip_time_do3d,                                                              &
+                      skip_time_domask,                                                            &
+                      synchronous_exchange,                                                        &
+                      termination_time_needed
+
+    NAMELIST /runtime_parameters/  averaging_interval,                                             &
+                                   averaging_interval_pr,                                          &
+                                   cpu_log_barrierwait,                                            &
+                                   create_disturbances,                                            &
+                                   cross_profiles,                                                 &
+                                   data_output,                                                    &
+                                   data_output_2d_on_each_pe,                                      &
+                                   data_output_masks,                                              &
+                                   data_output_pr,                                                 &
+                                   debug_output,                                                   &
+                                   debug_output_timestep,                                          &
+                                   disturbance_amplitude,                                          &
+                                   disturbance_energy_limit,                                       &
+                                   disturbance_level_b,                                            &
+                                   disturbance_level_t,                                            &
+                                   do2d_at_begin,                                                  &
+                                   do3d_at_begin,                                                  &
+                                   dt,                                                             &
+                                   dt_averaging_input,                                             &
+                                   dt_averaging_input_pr,                                          &
+                                   dt_coupling,                                                    &
+                                   dt_data_output,                                                 &
+                                   dt_data_output_av,                                              &
+                                   dt_disturb,                                                     &
+                                   dt_domask,                                                      &
+                                   dt_dopr,                                                        &
+                                   dt_dopr_listing,                                                &
+                                   dt_dots,                                                        &
+                                   dt_do2d_xy,                                                     &
+                                   dt_do2d_xz,                                                     &
+                                   dt_do2d_yz,                                                     &
+                                   dt_do3d,                                                        &
+                                   dt_max,                                                         &
+                                   dt_restart,                                                     &
+                                   dt_run_control,                                                 &
+                                   end_time,                                                       &
+                                   force_print_header,                                             &
+                                   mask_k_over_surface,                                            &
+                                   mask_scale_x,                                                   &
+                                   mask_scale_y,                                                   &
+                                   mask_scale_z,                                                   &
+                                   mask_x,                                                         &
+                                   mask_y,                                                         &
+                                   mask_z,                                                         &
+                                   mask_x_loop,                                                    &
+                                   mask_y_loop,                                                    &
+                                   mask_z_loop,                                                    &
+                                   netcdf_data_format,                                             &
+                                   netcdf_deflate,                                                 &
+                                   normalizing_region,                                             &
+                                   npex,                                                           &
+                                   npey,                                                           &
+                                   nz_do3d,                                                        &
+                                   profile_columns,                                                &
+                                   profile_rows,                                                   &
+                                   restart_time,                                                   &
+                                   section_xy,                                                     &
+                                   section_xz,                                                     &
+                                   section_yz,                                                     &
+                                   restart_data_format,                                            &
+                                   restart_data_format_input,                                      &
+                                   restart_data_format_output,                                     &
+                                   skip_time_data_output,                                          &
+                                   skip_time_data_output_av,                                       &
+                                   skip_time_dopr,                                                 &
+                                   skip_time_do2d_xy,                                              &
+                                   skip_time_do2d_xz,                                              &
+                                   skip_time_do2d_yz,                                              &
+                                   skip_time_do3d,                                                 &
+                                   skip_time_domask,                                               &
+                                   synchronous_exchange,                                           &
+                                   termination_time_needed
+
+    NAMELIST /envpar/  host,                                                                       &
+                       maximum_cpu_time_allowed,                                                   &
+                       maximum_parallel_io_streams,                                                &
+                       progress_bar_disabled,                                                      &
+                       read_svf,                                                                   &
+                       revision,                                                                   &
+                       run_identifier,                                                             &
+                       tasks_per_node,                                                             &
+                       write_binary,                                                               &
+                       write_svf
+
+!
+!-- First read values of environment variables (this NAMELIST file is generated by palmrun)
     CALL location_message( 'reading environment parameters from ENVPAR', 'start' )
 
-    OPEN ( 90, FILE='ENVPAR', STATUS='OLD', FORM='FORMATTED', IOSTAT=ioerr )
+    OPEN( 90, FILE='ENVPAR', STATUS='OLD', FORM='FORMATTED', IOSTAT=ioerr )
 
     IF ( ioerr /= 0 )  THEN
-       message_string = 'local file ENVPAR not found' //                       &
+       message_string = 'local file ENVPAR not found' //                                           &
                         '&some variables for steering may not be properly set'
        CALL message( 'parin', 'PA0276', 0, 1, 0, 6, 0 )
     ELSE
-       READ ( 90, envpar, IOSTAT=ioerr )
+       READ( 90, envpar, IOSTAT=ioerr )
        IF ( ioerr < 0 )  THEN
-          message_string = 'no envpar-NAMELIST found in local file '  //       &
+          message_string = 'no envpar-NAMELIST found in local file '  //                           &
                            'ENVPAR& or some variables for steering may '  //   &
                            'not be properly set'
           CALL message( 'parin', 'PA0278', 0, 1, 0, 6, 0 )
        ELSEIF ( ioerr > 0 )  THEN
-          message_string = 'errors in local file ENVPAR' //                    &
+          message_string = 'errors in local file ENVPAR' //                                        &
                            '&some variables for steering may not be properly set'
           CALL message( 'parin', 'PA0277', 0, 1, 0, 6, 0 )
        ENDIF
-       CLOSE ( 90 )
+       CLOSE( 90 )
     ENDIF
 
@@ -414,8 +739,8 @@
 
  11 REWIND ( 11 )
-    READ ( 11, inipar, ERR=12, END=13 )
-
-    message_string = 'namelist inipar is deprecated and will be ' //    &
-                     'removed in near future. & Please use namelist ' // &
+    READ( 11, inipar, ERR=12, END=13 )
+
+    message_string = 'namelist inipar is deprecated and will be ' //                               &
+                     'removed in near future. & Please use namelist ' //                           &
                      'initialization_parameters instead'
     CALL message( 'parin', 'PA0017', 0, 1, 0, 6, 0 )
@@ -431,20 +756,19 @@
 
 !
-!-- Try to read runtime parameters given by the user for this run
-!-- (namelist "runtime_parameters"). The namelist "runtime_parmeters"
-!-- can be omitted. In that case default values are used for the
+!-- Try to read runtime parameters given by the user for this run (namelist "runtime_parameters").
+!-- The namelist "runtime_parmeters" can be omitted. In that case default values are used for the
 !-- parameters.
  14 line = ' '
 
-    REWIND ( 11 )
+    REWIND( 11 )
     line = ' '
     DO WHILE ( INDEX( line, '&runtime_parameters' ) == 0 )
-       READ ( 11, '(A)', END=16 )  line
+       READ( 11, '(A)', END=16 )  line
     ENDDO
-    BACKSPACE ( 11 )
+    BACKSPACE( 11 )
 
 !
 !-- Read namelist
-    READ ( 11, runtime_parameters, ERR = 15 )
+    READ( 11, runtime_parameters, ERR = 15 )
     GOTO 18
 
@@ -453,8 +777,8 @@
     CALL parin_fail_message( 'runtime_parameters', line )
 
- 16 REWIND ( 11 )
+ 16 REWIND( 11 )
     line = ' '
     DO WHILE ( INDEX( line, '&d3par' ) == 0 )
-       READ ( 11, '(A)', END=18 )  line
+       READ( 11, '(A)', END=18 )  line
     ENDDO
     BACKSPACE ( 11 )
@@ -462,8 +786,8 @@
 !
 !-- Read namelist
-    READ ( 11, d3par, ERR = 17, END = 18 )
-
-    message_string = 'namelist d3par is deprecated and will be ' //      &
-                     'removed in near future. &Please use namelist ' //   &
+    READ( 11, d3par, ERR = 17, END = 18 )
+
+    message_string = 'namelist d3par is deprecated and will be ' //                                &
+                     'removed in near future. &Please use namelist ' //                            &
                      'runtime_parameters instead'
     CALL message( 'parin', 'PA0487', 0, 1, 0, 6, 0 )
@@ -483,13 +807,11 @@
 !
 !-- Calculate the number of groups into which parallel I/O is split.
-!-- The default for files which are opened by all PEs (or where each
-!-- PE opens his own independent file) is, that all PEs are doing input/output
-!-- in parallel at the same time. This might cause performance or even more
-!-- severe problems depending on the configuration of the underlying file
+!-- The default for files which are opened by all PEs (or where each PE opens its own independent
+!-- file) is, that all PEs are doing input/output in parallel at the same time. This might cause
+!-- performance or even more severe problems depending on the configuration of the underlying file
 !-- system.
-!-- Calculation of the number of blocks and the I/O group must be based on all
-!-- PEs involved in this run. Since myid and numprocs are related to the
-!-- comm2d communicator, which gives only a subset of all PEs in case of
-!-- nested runs, that information must be inquired again from the global
+!-- Calculation of the number of blocks and the I/O group must be based on all PEs involved in this
+!-- run. Since myid and numprocs are related to the comm2d communicator, which gives only a subset
+!-- of all PEs in case of nested runs, that information must be inquired again from the global
 !-- communicator.
 !-- First, set the default:
@@ -501,26 +823,23 @@
     global_procs = 1
 #endif
-    IF ( maximum_parallel_io_streams == -1  .OR.                               &
-         maximum_parallel_io_streams > global_procs )  THEN
+    IF ( maximum_parallel_io_streams == -1  .OR.  maximum_parallel_io_streams > global_procs )  THEN
        maximum_parallel_io_streams = global_procs
     ENDIF
 !
-!-- Now calculate the number of io_blocks and the io_group to which the
-!-- respective PE belongs. I/O of the groups is done in serial, but in parallel
-!-- for all PEs belonging to the same group.
+!-- Now calculate the number of io_blocks and the io_group to which the respective PE belongs. I/O
+!-- of the groups is done in serial, but in parallel for all PEs belonging to the same group.
     io_blocks = global_procs / maximum_parallel_io_streams
     io_group  = MOD( global_id+1, io_blocks )
 
 !
-!-- If required, read control parameters from restart file (produced by
-!-- a prior run). All PEs are reading from file created by PE0 (see
-!-- check_open)
+!-- If required, read control parameters from restart file (produced by a prior run). All PEs are
+!-- reading from file created by PE0 (see check_open)
     IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
 
 !
-!--    If not set by the user in the runtime parameters, the data format for restart
-!--    input needs to be set now! This is normally done later in check parameters.
-       IF ( TRIM( restart_data_format ) /= 'fortran_binary'  .AND.                           &
-            TRIM( restart_data_format ) /= 'mpi'             .AND.                           &
+!--    If not set by the user in the runtime parameters, the data format for restart input needs to
+!--    be set now! This is normally done later in check parameters.
+       IF ( TRIM( restart_data_format ) /= 'fortran_binary'  .AND.                                 &
+            TRIM( restart_data_format ) /= 'mpi'             .AND.                                 &
             TRIM( restart_data_format ) /= 'mpi_shared_memory' )  THEN
           message_string = 'illegal restart data format "' // TRIM( restart_data_format ) // '"'
@@ -553,13 +872,12 @@
        runnr = runnr + 1
 !
-!--    In case of a restart run, the number of user-defined profiles on
-!--    the restart file (already stored in max_pr_user) has to match the
-!--    one given for the current run. max_pr_user_tmp is calculated in
-!--    user_parin and max_pr_user is read in via rrd_global.
+!--    In case of a restart run, the number of user-defined profiles on the restart file (already
+!--    stored in max_pr_user) has to match the one given for the current run. max_pr_user_tmp is
+!--    calculated in user_parin and max_pr_user is read in via rrd_global.
        IF ( max_pr_user /= max_pr_user_tmp )  THEN
-          WRITE( message_string, * ) 'the number of user-defined ',      &
-                 'profiles given in data_output_pr (', max_pr_user_tmp,   &
-                 ') does not match the one ',                             &
-                 'found in the restart file (', max_pr_user, ')'
+          WRITE( message_string, * ) 'the number of user-defined ',                                &
+                                     'profiles given in data_output_pr (', max_pr_user_tmp,        &
+                                     ') does not match the one ',                                  &
+                                     'found in the restart file (', max_pr_user, ')'
           CALL message( 'user_parin', 'UI0009', 1, 2, 0, 6, 0 )
        ENDIF
@@ -569,5 +887,5 @@
 !
 !-- Activate spinup
-    IF ( land_surface .OR. urban_surface )  THEN
+    IF ( land_surface  .OR.  urban_surface )  THEN
        IF ( spinup_time > 0.0_wp )  THEN
           coupling_start_time = spinup_time
@@ -579,6 +897,6 @@
           IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
              spinup = .TRUE.
-          ELSEIF (        TRIM( initializing_actions ) == 'read_restart_data'      &
-                   .AND.  time_since_reference_point > 0.0_wp )  THEN
+          ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .AND.                      &
+                   time_since_reference_point > 0.0_wp )  THEN
              data_output_during_spinup = .FALSE.  !< required for correct ntdim calculation
                                                   !< in check_parameters for restart run
@@ -593,7 +911,7 @@
     IF ( nested_run )  THEN
        IF ( nesting_mode == 'vertical' )  THEN
-          IF (bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' )  THEN
-             WRITE ( message_string, *) 'bc_lr and bc_ns were set to ,', &
-                     'cyclic for vertical nesting'
+          IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
+             WRITE ( message_string, *) 'bc_lr and bc_ns were set to ,',                           &
+                                        'cyclic for vertical nesting'
              CALL message( 'parin', 'PA0428', 0, 0, 0, 6, 0 )
              bc_lr   = 'cyclic'
@@ -609,6 +927,5 @@
           ENDIF
 !
-!--    For other nesting modes only set boundary conditions for
-!--    nested domains.
+!--    For other nesting modes only set boundary conditions for nested domains.
        ELSE
           IF ( child_domain )  THEN
@@ -625,6 +942,5 @@
     ENDIF
 !
-!-- Set boundary conditions also in case the model is offline-nested in
-!-- larger-scale models.
+!-- Set boundary conditions also in case the model is offline-nested in larger-scale models.
     IF ( nesting_offline )  THEN
        bc_lr    = 'nesting_offline'
@@ -639,19 +955,15 @@
 
 !
-!-- In case of nested runs, make sure that initializing_actions =
-!-- 'set_constant_profiles' even though the constant-profiles
-!-- initializations for the prognostic variables will be overwritten
-!-- by pmci_child_initialize and pmci_parent_initialize. This is,
-!-- however, important e.g. to make sure that diagnostic variables
-!-- are set properly. An exception is made in case of restart runs and
-!-- if user decides to do everything by its own.
-    IF ( child_domain  .AND.  .NOT. (                                    &
-         TRIM( initializing_actions ) == 'read_restart_data'      .OR.   &
-         TRIM( initializing_actions ) == 'set_constant_profiles'  .OR.   &
-         TRIM( initializing_actions ) == 'by_user' ) )  THEN
-       message_string = 'initializing_actions = ' //                     &
-                        TRIM( initializing_actions ) // ' has been ' //  &
-                        'changed to set_constant_profiles in child ' //  &
-                        'domain.'
+!-- In case of nested runs, make sure that initializing_actions = 'set_constant_profiles' even
+!-- though the constant-profiles initializations for the prognostic variables will be overwritten by
+!-- pmci_child_initialize and pmci_parent_initialize. This is, however, important e.g. to make sure
+!-- that diagnostic variables are set properly. An exception is made in case of restart runs and if
+!-- user decides to do everything by its own.
+    IF ( child_domain  .AND.  .NOT. ( TRIM( initializing_actions ) == 'read_restart_data'     .OR. &
+                                      TRIM( initializing_actions ) == 'set_constant_profiles' .OR. &
+                                      TRIM( initializing_actions ) == 'by_user' )                  &
+                                    )  THEN
+       message_string = 'initializing_actions = ' // TRIM( initializing_actions ) //               &
+                        ' has been ' // 'changed to set_constant_profiles in child ' // 'domain.'
        CALL message( 'parin', 'PA0492', 0, 0, 0, 6, 0 )
 
@@ -659,19 +971,16 @@
     ENDIF
 !
-!-- Check validity of lateral boundary conditions. This has to be done
-!-- here because they are already used in init_pegrid and init_grid and
-!-- therefore cannot be check in check_parameters
-    IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/radiation'  .AND. &
-         bc_lr /= 'radiation/dirichlet'  .AND.  bc_lr /= 'nested'  .AND. &
+!-- Check validity of lateral boundary conditions. This has to be done here because they are already
+!-- used in init_pegrid and init_grid and therefore cannot be check in check_parameters
+    IF ( bc_lr /= 'cyclic'               .AND.  bc_lr /= 'dirichlet/radiation'  .AND.              &
+         bc_lr /= 'radiation/dirichlet'  .AND.  bc_lr /= 'nested'               .AND.              &
          bc_lr /= 'nesting_offline' )  THEN
-       message_string = 'unknown boundary condition: bc_lr = "' // &
-                        TRIM( bc_lr ) // '"'
+       message_string = 'unknown boundary condition: bc_lr = "' // TRIM( bc_lr ) // '"'
        CALL message( 'parin', 'PA0049', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/radiation'  .AND. &
-         bc_ns /= 'radiation/dirichlet'  .AND.  bc_ns /= 'nested'  .AND. &
+    IF ( bc_ns /= 'cyclic'               .AND.  bc_ns /= 'dirichlet/radiation'  .AND.              &
+         bc_ns /= 'radiation/dirichlet'  .AND.  bc_ns /= 'nested'               .AND.              &
          bc_ns /= 'nesting_offline' )  THEN
-       message_string = 'unknown boundary condition: bc_ns = "' // &
-                        TRIM( bc_ns ) // '"'
+       message_string = 'unknown boundary condition: bc_ns = "' // TRIM( bc_ns ) // '"'
        CALL message( 'parin', 'PA0050', 1, 2, 0, 6, 0 )
     ENDIF
@@ -686,40 +995,35 @@
 !
 !-- Radiation-Dirichlet conditions are allowed along one of the horizontal directions only.
-!-- In general, such conditions along x and y may work, but require a) some code changes
-!-- (e.g. concerning allocation of c_u, c_v ... arrays), and b) a careful model setup by the
-!-- user, in order to guarantee that there is a clearly defined outflow at two sides.
+!-- In general, such conditions along x and y may work, but require a) some code changes (e.g.
+!-- concerning allocation of c_u, c_v ... arrays), and b) a careful model setup by the user, in
+!-- order to guarantee that there is a clearly defined outflow at two sides.
 !-- Otherwise, the radiation condition may produce wrong results.
     IF ( ( bc_lr_dirrad .OR. bc_lr_raddir )  .AND.  ( bc_ns_dirrad .OR. bc_ns_raddir ) )  THEN
-       message_string = 'bc_lr = "' // TRIM( bc_lr ) // '" and bc_ns = "' //                 &
-                        TRIM( bc_ns ) // '" are not allowed to be set at the same time'
+       message_string = 'bc_lr = "' // TRIM( bc_lr ) // '" and bc_ns = "' // TRIM( bc_ns ) //      &
+                        '" are not allowed to be set at the same time'
        CALL message( 'parin', 'PA0589', 1, 2, 0, 6, 0 )
     ENDIF
 !
-!-- Check in case of initial run, if the grid point numbers are well
-!-- defined and allocate some arrays which are already needed in
-!-- init_pegrid or check_parameters. During restart jobs, these arrays
-!-- will be allocated in rrd_global. All other arrays are allocated
-!-- in init_3d_model.
+!-- Check in case of initial run, if the grid point numbers are well defined and allocate some
+!-- arrays which are already needed in init_pegrid or check_parameters. During restart jobs, these
+!-- arrays will be allocated in rrd_global. All other arrays are allocated in init_3d_model.
     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 
        IF ( nx <= 0 )  THEN
-          WRITE( message_string, * ) 'no value or wrong value given',    &
-                                     ' for nx: nx=', nx
+          WRITE( message_string, * ) 'no value or wrong value given', ' for nx: nx=', nx
           CALL message( 'parin', 'PA0273', 1, 2, 0, 6, 0 )
        ENDIF
        IF ( ny <= 0 )  THEN
-          WRITE( message_string, * ) 'no value or wrong value given',    &
-                                     ' for ny: ny=', ny
+          WRITE( message_string, * ) 'no value or wrong value given', ' for ny: ny=', ny
           CALL message( 'parin', 'PA0274', 1, 2, 0, 6, 0 )
        ENDIF
        IF ( nz <= 0 )  THEN
-          WRITE( message_string, * ) 'no value or wrong value given',    &
-                                     ' for nz: nz=', nz
+          WRITE( message_string, * ) 'no value or wrong value given', ' for nz: nz=', nz
           CALL message( 'parin', 'PA0275', 1, 2, 0, 6, 0 )
        ENDIF
 
 !
-!--    As a side condition, routine flow_statistics require at least 14
-!--    vertical grid levels (they are used to store time-series data)
+!--    As a side condition, routine flow_statistics require at least 14 vertical grid levels (they
+!--    are used to store time-series data)
 !>     @todo   Remove this restriction
        IF ( nz < 14 )  THEN
@@ -729,10 +1033,10 @@
 
 !
-!--    ATTENTION: in case of changes to the following statement please
-!--               also check the allocate statement in routine rrd_global
-       ALLOCATE( pt_init(0:nz+1), q_init(0:nz+1), s_init(0:nz+1),        &
-                 ref_state(0:nz+1), sa_init(0:nz+1), ug(0:nz+1),         &
-                 u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1),             &
-                 hom(0:nz+1,2,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:statistic_regions),  &
+!--    ATTENTION: in case of changes to the following statement please also check the allocate
+!--               statement in routine rrd_global
+       ALLOCATE( pt_init(0:nz+1), q_init(0:nz+1), s_init(0:nz+1),                                  &
+                 ref_state(0:nz+1), sa_init(0:nz+1), ug(0:nz+1),                                   &
+                 u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1),                                       &
+                 hom(0:nz+1,2,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:statistic_regions),     &
                  hom_sum(0:nz+1,pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0:statistic_regions) )
 
