Index: palm/trunk/SOURCE/ocean_mod.f90
===================================================================
--- palm/trunk/SOURCE/ocean_mod.f90	(revision 4768)
+++ 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 )
 
