Changeset 4559

Timestamp:

Jun 11, 2020 8:51:48 AM (11 months ago)

Author:

raasch

Message:

files re-formatted to follow the PALM coding standard

Location:

palm/trunk/SOURCE

Files:

• check_parameters.f90 (modified) (166 diffs)
• chem_emissions_mod.f90 (modified) (106 diffs)
• chem_modules.f90 (modified) (8 diffs)
• chem_photolysis_mod.f90 (modified) (10 diffs)
• chemistry_model_mod.f90 (modified) (267 diffs)
• close_file.f90 (modified) (6 diffs)
• compute_vpt.f90 (modified) (4 diffs)
• coriolis.f90 (modified) (10 diffs)
• cpulog_mod.f90 (modified) (10 diffs)
• cuda_fft_interfaces.f90 (modified) (9 diffs)
• data_log.f90 (modified) (5 diffs)
• data_output_2d.f90 (modified) (118 diffs)
• data_output_3d.f90 (modified) (39 diffs)
• data_output_binary_module.f90 (modified) (27 diffs)
• data_output_flight.f90 (modified) (6 diffs)
• data_output_mask.f90 (modified) (35 diffs)
• surface_mod.f90 (modified) (22 diffs)
• synthetic_turbulence_generator_mod.f90 (modified) (89 diffs)
• temperton_fft_mod.f90 (modified) (163 diffs)

palm/trunk/SOURCE/check_parameters.f90

@@ -1 +1 @@
 !> @file check_parameters.f90
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! This file is part of the PALM model system.
 !
-! PALM is free software: you can redistribute it and/or modify it under the
-! terms of the GNU General Public License as published by the Free Software
-! Foundation, either version 3 of the License, or (at your option) any later
-! version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
+! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
+! Public License as published by the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
 !
 ! Copyright 1997-2020 Leibniz Universitaet Hannover
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 !
 ! Current revisions:
@@ -25 +24 @@
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+! 4536 2020-05-17 17:24:13Z raasch
 ! unneccessary query for restart data format removed
-! 
+!
 ! 4534 2020-05-14 18:35:22Z raasch
 ! adjustments for I/O on reduced number of cores using shared memory MPI
-! 
+!
 ! 4514 2020-04-30 16:29:59Z suehring
 ! Enable output of qsurf and ssurf
-! 
+!
 ! 4513 2020-04-30 13:45:47Z raasch
 ! unused modules removed
-! 
+!
 ! 4511 2020-04-30 12:20:40Z raasch
 ! call of chem_boundary_conds removed (respective settings are now done in the chemistry module)
-! 
+!
 ! 4495 2020-04-13 20:11:20Z raasch
 ! check new restart_data_format parameters
-! 
+!
 ! 4493 2020-04-10 09:49:43Z pavelkrc
 ! Kolmogorov length scale eta added to profile output
@@ -134 +136 @@
 ! ------------
 !> Check control parameters and deduce further quantities.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE check_parameters
 
@@ -142 +144 @@
     USE basic_constants_and_equations_mod
 
-    USE bulk_cloud_model_mod,                                                  &
+    USE bulk_cloud_model_mod,                                                                      &
         ONLY:  bulk_cloud_model
 
@@ -153 +155 @@
     USE indices
 
-    USE model_1d_mod,                                                          &
+    USE model_1d_mod,                                                                              &
         ONLY:  damp_level_1d, damp_level_ind_1d
 
-    USE module_interface,                                                      &
-        ONLY:  module_interface_check_parameters,                              &
-               module_interface_check_data_output_ts,                          &
-               module_interface_check_data_output_pr,                          &
-               module_interface_check_data_output
-
-    USE netcdf_data_input_mod,                                                 &
-        ONLY:  init_model, input_pids_static, netcdf_data_input_check_dynamic, &
+    USE module_interface,                                                                          &
+        ONLY:  module_interface_check_data_output,                                                 &
+               module_interface_check_data_output_pr,                                              &
+               module_interface_check_data_output_ts,                                              &
+               module_interface_check_parameters
+
+    USE netcdf_data_input_mod,                                                                     &
+        ONLY:  init_model, input_pids_static, netcdf_data_input_check_dynamic,                     &
                netcdf_data_input_check_static
 
-    USE netcdf_interface,                                                      &
-        ONLY:  dopr_unit, do2d_unit, do3d_unit, netcdf_data_format,            &
-               netcdf_data_format_string, dots_unit, heatflux_output_unit,     &
-               waterflux_output_unit, momentumflux_output_unit,                &
-               dots_max, dots_num, dots_label
-
-    USE particle_attributes,                                                   &
+    USE netcdf_interface,                                                                          &
+        ONLY:  do2d_unit, do3d_unit, dopr_unit, dots_label, dots_max, dots_num, dots_unit,         &
+               heatflux_output_unit, momentumflux_output_unit, netcdf_data_format,                 &
+               netcdf_data_format_string, waterflux_output_unit
+
+    USE particle_attributes,                                                                       &
         ONLY:  particle_advection, use_sgs_for_particles
 
     USE pegrid
 
-    USE pmc_interface,                                                         &
+    USE pmc_interface,                                                                             &
         ONLY:  cpl_id, nested_run
 
@@ -189 +190 @@
 
 #if defined( __parallel )
-    USE vertical_nesting_mod,                                                  &
-        ONLY:  vnested,                                                        &
+    USE vertical_nesting_mod,                                                                      &
+        ONLY:  vnested,                                                                            &
                vnest_check_parameters
 #endif
@@ -228 +229 @@
     CALL location_message( 'checking parameters', 'start' )
 !
-!-- At first, check static and dynamic input for consistency
+!-- At first, check static and dynamic input for consistency.
     CALL netcdf_data_input_check_dynamic
     CALL netcdf_data_input_check_static
@@ -240 +241 @@
 !
 !-- Check the coupling mode
-    IF ( coupling_mode /= 'uncoupled'            .AND.                         &
-         coupling_mode /= 'precursor_atmos'      .AND.                         &
-         coupling_mode /= 'precursor_ocean'      .AND.                         &
-         coupling_mode /= 'vnested_crse'         .AND.                         &
-         coupling_mode /= 'vnested_fine'         .AND.                         &
-         coupling_mode /= 'atmosphere_to_ocean'  .AND.                         &
+    IF ( coupling_mode /= 'uncoupled'            .AND.                                             &
+         coupling_mode /= 'precursor_atmos'      .AND.                                             &
+         coupling_mode /= 'precursor_ocean'      .AND.                                             &
+         coupling_mode /= 'vnested_crse'         .AND.                                             &
+         coupling_mode /= 'vnested_fine'         .AND.                                             &
+         coupling_mode /= 'atmosphere_to_ocean'  .AND.                                             &
          coupling_mode /= 'ocean_to_atmosphere' )  THEN
        message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
@@ -252 +253 @@
 
 !
-!-- Check if humidity is set to TRUE in case of the atmospheric run (for coupled runs)
+!-- Check if humidity is set to .TRUE. in case of the atmospheric run (for coupled runs)
     IF ( coupling_mode == 'atmosphere_to_ocean' .AND. .NOT. humidity) THEN
-       message_string = ' Humidity has to be set to .T. in the _p3d file ' //  &
+       message_string = ' Humidity has to be set to .T. in the _p3d file ' //                      &
                         'for coupled runs between ocean and atmosphere.'
        CALL message( 'check_parameters', 'PA0476', 1, 2, 0, 6, 0 )
@@ -297 +298 @@
 !
 !-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
-    IF ( coupling_mode /= 'uncoupled'       .AND.                              &
-         coupling_mode(1:8) /= 'vnested_'   .AND.                              &
-         coupling_mode /= 'precursor_atmos' .AND.                              &
+    IF ( coupling_mode /= 'uncoupled'        .AND.                                                 &
+         coupling_mode(1:8) /= 'vnested_'    .AND.                                                 &
+         coupling_mode /= 'precursor_atmos'  .AND.                                                 &
          coupling_mode /= 'precursor_ocean' )  THEN
 
        IF ( dt_coupling == 9999999.9_wp )  THEN
-          message_string = 'dt_coupling is not set but required for coup' //   &
-                           'ling mode "' //  TRIM( coupling_mode ) // '"'
+          message_string = 'dt_coupling is not set but required for coupling mode "' //            &
+                           TRIM( coupling_mode ) // '"'
           CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
        ENDIF
@@ -312 +313 @@
 
        IF ( myid == 0 ) THEN
-          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter,  &
-                         ierr )
-          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter,       &
-                         status, ierr )
+          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, status, ierr )
        ENDIF
        CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
 
        IF ( dt_coupling /= remote )  THEN
-          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
-                 '": dt_coupling = ', dt_coupling, '& is not equal to ',       &
-                 'dt_coupling_remote = ', remote
+          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), '": dt_coupling = ',&
+                 dt_coupling, '& is not equal to ', 'dt_coupling_remote = ', remote
           CALL message( 'check_parameters', 'PA0004', 1, 2, 0, 6, 0 )
        ENDIF
@@ -329 +327 @@
           IF ( myid == 0  ) THEN
              CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
-             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter,    &
-                            status, ierr )
+             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, status, ierr )
           ENDIF
           CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
 
           dt_coupling = MAX( dt_max, remote )
-          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
-                 '": dt_coupling <= 0.0 & is not allowed and is reset to ',    &
-                 'MAX(dt_max(A,O)) = ', dt_coupling
+          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                     &
+                 '": dt_coupling <= 0.0 & is not allowed and is reset to ', 'MAX(dt_max(A,O)) = ', &
+                 dt_coupling
           CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
        ENDIF
 
        IF ( myid == 0 ) THEN
-          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
-                         ierr )
-          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter,       &
-                         status, ierr )
+          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, status, ierr )
        ENDIF
        CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
 
        IF ( restart_time /= remote )  THEN
-          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
-                 '": restart_time = ', restart_time, '& is not equal to ',     &
+          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                     &
+                 '": restart_time = ', restart_time, '& is not equal to ',                         &
                  'restart_time_remote = ', remote
           CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
@@ -357 +352 @@
 
        IF ( myid == 0 ) THEN
-          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter,   &
-                         ierr )
-          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter,       &
-                         status, ierr )
+          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, status, ierr )
        ENDIF
        CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
 
        IF ( dt_restart /= remote )  THEN
-          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
-                 '": dt_restart = ', dt_restart, '& is not equal to ',         &
-                 'dt_restart_remote = ', remote
+          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), '": dt_restart = ', &
+                 dt_restart, '& is not equal to ', 'dt_restart_remote = ', remote
           CALL message( 'check_parameters', 'PA0007', 1, 2, 0, 6, 0 )
        ENDIF
@@ -373 +365 @@
        time_to_be_simulated_from_reference_point = end_time-coupling_start_time
 
-       IF  ( myid == 0 ) THEN
-          CALL MPI_SEND( time_to_be_simulated_from_reference_point, 1,         &
-                         MPI_REAL, target_id, 14, comm_inter, ierr )
-          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter,       &
-                         status, ierr )
+       IF ( myid == 0 ) THEN
+          CALL MPI_SEND( time_to_be_simulated_from_reference_point, 1, MPI_REAL, target_id, 14,    &
+                         comm_inter, ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, status, ierr )
        ENDIF
        CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
 
        IF ( time_to_be_simulated_from_reference_point /= remote )  THEN
-          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
-                 '": time_to_be_simulated_from_reference_point = ',            &
-                 time_to_be_simulated_from_reference_point, '& is not equal ', &
-                 'to time_to_be_simulated_from_reference_point_remote = ',     &
+          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                     &
+                 '": time_to_be_simulated_from_reference_point = ',                                &
+                 time_to_be_simulated_from_reference_point, '& is not equal ',                     &
+                 'to time_to_be_simulated_from_reference_point_remote = ',                         &
                  remote
           CALL message( 'check_parameters', 'PA0008', 1, 2, 0, 6, 0 )
@@ -392 +383 @@
        IF ( myid == 0 ) THEN
           CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
-          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter,       &
-                                                             status, ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, status, ierr )
        ENDIF
        CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
@@ -401 +391 @@
 
           IF ( dx < remote ) THEN
-             WRITE( message_string, * ) 'coupling mode "',                     &
-                   TRIM( coupling_mode ),                                      &
-           '": dx in Atmosphere is not equal to or not larger than dx in ocean'
+             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
+                    '": dx in Atmosphere is not equal to or not larger than dx in ocean'
              CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
           ENDIF
 
           IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
-             WRITE( message_string, * ) 'coupling mode "',                     &
-                    TRIM( coupling_mode ),                                     &
-             '": Domain size in x-direction is not equal in ocean and atmosphere'
+             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
+                    '": Domain size in x-direction is not equal in ocean and atmosphere'
              CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
           ENDIF
@@ -418 +406 @@
        IF ( myid == 0) THEN
           CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
-          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter,       &
-                         status, ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, status, ierr )
        ENDIF
        CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
@@ -426 +413 @@
 
           IF ( dy < remote )  THEN
-             WRITE( message_string, * ) 'coupling mode "',                     &
-                    TRIM( coupling_mode ),                                     &
-                 '": dy in Atmosphere is not equal to or not larger than dy in ocean'
+             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
+                    '": dy in Atmosphere is not equal to or not larger than dy in ocean'
              CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
           ENDIF
 
           IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
-             WRITE( message_string, * ) 'coupling mode "',                     &
-                   TRIM( coupling_mode ),                                      &
-             '": Domain size in y-direction is not equal in ocean and atmosphere'
+             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
+                    '": Domain size in y-direction is not equal in ocean and atmosphere'
              CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
           ENDIF
 
-          IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
-             WRITE( message_string, * ) 'coupling mode "',                     &
-                   TRIM( coupling_mode ),                                      &
-             '": nx+1 in ocean is not divisible by nx+1 in',                   &
-             ' atmosphere without remainder'
+          IF ( MOD( nx_o+1, nx_a+1 ) /= 0 )  THEN
+             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
+                    '": nx+1 in ocean is not divisible by nx+1 in', ' atmosphere without remainder'
              CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
           ENDIF
 
           IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
-             WRITE( message_string, * ) 'coupling mode "',                     &
-                   TRIM( coupling_mode ),                                      &
-             '": ny+1 in ocean is not divisible by ny+1 in', &
-             ' atmosphere without remainder'
+             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
+                    '": ny+1 in ocean is not divisible by ny+1 in', ' atmosphere without remainder'
 
              CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )
@@ -458 +439 @@
        ENDIF
 #else
-       WRITE( message_string, * ) 'coupling requires PALM to be compiled with',&
-            ' cpp-option "-D__parallel"'
+       WRITE( message_string, * ) 'coupling requires PALM to be compiled with',                    &
+              ' cpp-option "-D__parallel"'
        CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
 #endif
@@ -467 +448 @@
 !
 !-- Exchange via intercommunicator
-    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 )  THEN
-       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter,     &
-                      ierr )
-    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0)  THEN
-       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19,          &
-                      comm_inter, status, ierr )
+    IF ( coupling_mode == 'atmosphere_to_ocean'  .AND.  myid == 0 )  THEN
+       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, ierr )
+    ELSEIF ( coupling_mode == 'ocean_to_atmosphere'  .AND.  myid == 0)  THEN
+       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, comm_inter, status, ierr )
     ENDIF
     CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
@@ -479 +458 @@
 
 !
-!-- User settings for restart times requires that "restart" has been given as
-!-- file activation string. Otherwise, binary output would not be saved by
-!-- palmrun.
-    IF (  ( restart_time /= 9999999.9_wp  .OR.  dt_restart /= 9999999.9_wp )   &
+!-- User settings for restart times requires that "restart" has been given as file activation
+!-- string. Otherwise, binary output would not be saved by palmrun.
+    IF ( ( restart_time /= 9999999.9_wp  .OR.  dt_restart /= 9999999.9_wp )                        &
          .AND.  .NOT. write_binary )  THEN
-       WRITE( message_string, * ) 'manual restart settings requires file ',    &
+       WRITE( message_string, * ) 'manual restart settings requires file ',                        &
                                   'activation string "restart"'
        CALL message( 'check_parameters', 'PA0001', 1, 2, 0, 6, 0 )
@@ -491 +469 @@
 
 !
-!-- Generate the file header which is used as a header for most of PALM's
-!-- output files
+!-- Generate the file header which is used as a header for most of PALM's output files
     CALL DATE_AND_TIME( date, time, run_zone )
-    run_date = date(1:4)//'-'//date(5:6)//'-'//date(7:8)
-    run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
+    run_date = date(1:4) // '-' // date(5:6) // '-' // date(7:8)
+    run_time = time(1:2) // ':' // time(3:4) // ':' // time(5:6)
     IF ( coupling_mode == 'uncoupled' )  THEN
        coupling_string = ''
@@ -518 +495 @@
     ENDIF
 
-    WRITE ( run_description_header,                                            &
-            '(A,2X,A,2X,A,A,A,I2.2,A,A,A,2X,A,A,2X,A,1X,A)' )                  &
-          TRIM( version ), TRIM( revision ), 'run: ',                          &
-          TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ),         &
-          TRIM( nest_string ), TRIM( ensemble_string), 'host: ', TRIM( host ), &
-          run_date, run_time
+    WRITE ( run_description_header, '(A,2X,A,2X,A,A,A,I2.2,A,A,A,2X,A,A,2X,A,1X,A)' )              &
+          TRIM( version ), TRIM( revision ), 'run: ', TRIM( run_identifier ), '.', runnr,          &
+          TRIM( coupling_string ), TRIM( nest_string ), TRIM( ensemble_string), 'host: ',          &
+          TRIM( host ), run_date, run_time
 
 !
@@ -533 +508 @@
 
        CASE DEFAULT
-          message_string = 'illegal value given for loop_optimization: "' //   &
+          message_string = 'illegal value given for loop_optimization: "' //                       &
                            TRIM( loop_optimization ) // '"'
           CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
@@ -543 +518 @@
     IF ( topography /= 'flat' )  THEN
        action = ' '
-       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme'      &
-          )  THEN
+       IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' )  THEN
           WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
        ENDIF
-       IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' )&
-       THEN
+       IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' )  THEN
           WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
        ENDIF
@@ -563 +536 @@
           WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
        ENDIF
-       IF ( .NOT. constant_flux_layer .AND. topography /= 'closed_channel' )   &
-       THEN
+       IF ( .NOT. constant_flux_layer  .AND.  topography /= 'closed_channel' )  THEN
           WRITE( action, '(A)' )  'constant_flux_layer = .FALSE.'
        ENDIF
        IF ( action /= ' ' )  THEN
-          message_string = 'The specified topography does not allow ' //       &
-                           TRIM( action )
+          message_string = 'The specified topography does not allow ' // TRIM( action )
           CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
        ENDIF
 !
 !--    Check illegal/untested parameter combinations for closed channel
-       If ( topography == 'closed_channel' ) THEN
+       If ( topography == 'closed_channel' )  THEN
           symmetry_flag = 1
           message_string = 'Bottom and top boundary are treated equal'
           CALL message( 'check_parameters', 'PA0699', 0, 0, 0, 6, 0 )
 
-          IF ( dz(1) /= dz(COUNT( dz /= -1.0_wp )) .OR.                        &
-               dz_stretch_level /= -9999999.9_wp) THEN
-             WRITE( message_string, * )  'dz should be equal close to the ' // &
+          IF ( dz(1) /= dz(COUNT( dz /= -1.0_wp ))  .OR.  dz_stretch_level /= -9999999.9_wp)  THEN
+             WRITE( message_string, * )  'dz should be equal close to the ' //                     &
                                          'boundaries due to symmetrical problem'
              CALL message( 'check_parameters', 'PA0700', 1, 2, 0, 6, 0 )
           ENDIF
 
-          IF ( constant_flux_layer ) THEN
-             WRITE( message_string, * )  'A constant flux layer is not '//     &
-                                         'allowed if a closed channel '//      &
-                                         'shall be used'
+          IF ( constant_flux_layer )  THEN
+             WRITE( message_string, * )  'A constant flux layer is not ' //                        &
+                                         'allowed if a closed channel shall be used'
              CALL message( 'check_parameters', 'PA0701', 1, 2, 0, 6, 0 )
           ENDIF
 
-          IF ( ocean_mode ) THEN
-             WRITE( message_string, * )  'The ocean mode is not allowed if '// &
+          IF ( ocean_mode )  THEN
+             WRITE( message_string, * )  'The ocean mode is not allowed if ' //                    &
                                          'a closed channel shall be used'
              CALL message( 'check_parameters', 'PA0702', 1, 2, 0, 6, 0 )
           ENDIF
 
-          IF ( momentum_advec /= 'ws-scheme' .OR.                              &
-               scalar_advec /= 'ws-scheme' ) THEN
-             WRITE( message_string, * )  'A closed channel require the '//     &
-                                         'upwind scheme of Wicker and ' //     &
-                                         'Skamarock as the advection scheme'
+          IF ( momentum_advec /= 'ws-scheme'  .OR.                                                 &
+               scalar_advec /= 'ws-scheme' )  THEN
+             WRITE( message_string, * )  'A closed channel require the upwind scheme of Wicker' // &
+                                         ' and Skamarock as the advection scheme'
              CALL message( 'check_parameters', 'PA0703', 1, 2, 0, 6, 0 )
           ENDIF
@@ -611 +579 @@
 !
 !-- Check approximation
-    IF ( TRIM( approximation ) /= 'boussinesq'   .AND.                         &
-         TRIM( approximation ) /= 'anelastic' )  THEN
-       message_string = 'unknown approximation: approximation = "' //          &
-                        TRIM( approximation ) // '"'
+    IF ( TRIM( approximation ) /= 'boussinesq'  .AND.  TRIM( approximation ) /= 'anelastic' )  THEN
+       message_string = 'unknown approximation: approximation = "' // TRIM( approximation ) // '"'
        CALL message( 'check_parameters', 'PA0446', 1, 2, 0, 6, 0 )
     ENDIF
@@ -620 +586 @@
 !
 !-- Check approximation requirements
-    IF ( TRIM( approximation ) == 'anelastic'   .AND.                          &
-         TRIM( momentum_advec ) /= 'ws-scheme' )  THEN
-       message_string = 'Anelastic approximation requires ' //                 &
-                        'momentum_advec = "ws-scheme"'
+    IF ( TRIM( approximation ) == 'anelastic'  .AND.  TRIM( momentum_advec ) /= 'ws-scheme' )  THEN
+       message_string = 'Anelastic approximation requires momentum_advec = "ws-scheme"'
        CALL message( 'check_parameters', 'PA0447', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( TRIM( approximation ) == 'anelastic'   .AND.                          &
-         TRIM( psolver ) == 'multigrid' )  THEN
-       message_string = 'Anelastic approximation currently only supports ' //  &
-                        'psolver = "poisfft", ' //                             &
-                        'psolver = "sor" and ' //                              &
-                        'psolver = "multigrid_noopt"'
+    IF ( TRIM( approximation ) == 'anelastic'  .AND.  TRIM( psolver ) == 'multigrid' )  THEN
+       message_string = 'Anelastic approximation currently only supports psolver = "poisfft", ' // &
+                        'psolver = "sor" and psolver = "multigrid_noopt"'
        CALL message( 'check_parameters', 'PA0448', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( TRIM( approximation ) == 'anelastic'   .AND.                          &
-         conserve_volume_flow )  THEN
-       message_string = 'Anelastic approximation is not allowed with ' //      &
+    IF ( TRIM( approximation ) == 'anelastic'  .AND.  conserve_volume_flow )  THEN
+       message_string = 'Anelastic approximation is not allowed with ' //                          &
                         'conserve_volume_flow = .TRUE.'
        CALL message( 'check_parameters', 'PA0449', 1, 2, 0, 6, 0 )
@@ -643 +603 @@
 !
 !-- Check flux input mode
-    IF ( TRIM( flux_input_mode ) /= 'dynamic'    .AND.                         &
-         TRIM( flux_input_mode ) /= 'kinematic'  .AND.                         &
-         TRIM( flux_input_mode ) /= 'approximation-specific' )  THEN
-       message_string = 'unknown flux input mode: flux_input_mode = "' //      &
+    IF ( TRIM( flux_input_mode ) /= 'dynamic'  .AND.  TRIM( flux_input_mode ) /= 'kinematic'       &
+         .AND.  TRIM( flux_input_mode ) /= 'approximation-specific' )  THEN
+       message_string = 'unknown flux input mode: flux_input_mode = "' //                          &
                         TRIM( flux_input_mode ) // '"'
        CALL message( 'check_parameters', 'PA0450', 1, 2, 0, 6, 0 )
@@ -662 +621 @@
 !
 !-- Check flux output mode
-    IF ( TRIM( flux_output_mode ) /= 'dynamic'    .AND.                        &
-         TRIM( flux_output_mode ) /= 'kinematic'  .AND.                        &
-         TRIM( flux_output_mode ) /= 'approximation-specific' )  THEN
-       message_string = 'unknown flux output mode: flux_output_mode = "' //    &
+    IF ( TRIM( flux_output_mode ) /= 'dynamic'  .AND.  TRIM( flux_output_mode ) /= 'kinematic'     &
+         .AND.  TRIM( flux_output_mode ) /= 'approximation-specific' )  THEN
+       message_string = 'unknown flux output mode: flux_output_mode = "' //                        &
                         TRIM( flux_output_mode ) // '"'
        CALL message( 'check_parameters', 'PA0451', 1, 2, 0, 6, 0 )
@@ -681 +639 @@
 
 !
-!-- When the land- or urban-surface model is used, the flux output must be
-!-- dynamic.
+!-- When the land- or urban-surface model is used, the flux output must be dynamic.
     IF ( land_surface  .OR.  urban_surface )  THEN
        flux_output_mode = 'dynamic'
@@ -689 +646 @@
 !
 !-- Set the flux output units according to flux_output_mode
-    IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN
+    IF ( TRIM( flux_output_mode ) == 'kinematic' )  THEN
         heatflux_output_unit              = 'K m/s'
         waterflux_output_unit             = 'kg/kg m/s'
         momentumflux_output_unit          = 'm2/s2'
-    ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
+    ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' )  THEN
         heatflux_output_unit              = 'W/m2'
         waterflux_output_unit             = 'W/m2'
@@ -712 +669 @@
 !-- Check if maximum number of allowed timeseries is exceeded
     IF ( dots_num > dots_max )  THEN
-       WRITE( message_string, * ) 'number of time series quantities exceeds',  &
-                                  ' its maximum of dots_max = ', dots_max,     &
+       WRITE( message_string, * ) 'number of time series quantities exceeds',                      &
+                                  ' its maximum of dots_max = ', dots_max,                         &
                                   '&Please increase dots_max in modules.f90.'
        CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )
@@ -721 +678 @@
 !-- Check whether there are any illegal values
 !-- Pressure solver:
-    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'sor'  .AND.                  &
-         psolver /= 'multigrid'  .AND.  psolver /= 'multigrid_noopt' )  THEN
-       message_string = 'unknown solver for perturbation pressure: psolver' // &
+    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'sor'  .AND.  psolver /= 'multigrid'  .AND.       &
+         psolver /= 'multigrid_noopt' )  THEN
+       message_string = 'unknown solver for perturbation pressure: psolver' //                     &
                         ' = "' // TRIM( psolver ) // '"'
        CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
@@ -734 +691 @@
           gamma_mg = 1
        ELSE
-          message_string = 'unknown multigrid cycle: cycle_mg = "' //          &
-                           TRIM( cycle_mg ) // '"'
+          message_string = 'unknown multigrid cycle: cycle_mg = "' //  TRIM( cycle_mg ) // '"'
           CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
        ENDIF
     ENDIF
 
-    IF ( fft_method /= 'singleton-algorithm'  .AND.                            &
-         fft_method /= 'temperton-algorithm'  .AND.                            &
-         fft_method /= 'fftw'                 .AND.                            &
-         fft_method /= 'system-specific' )  THEN
-       message_string = 'unknown fft-algorithm: fft_method = "' //             &
-                        TRIM( fft_method ) // '"'
+    IF ( fft_method /= 'singleton-algorithm'  .AND.  fft_method /= 'temperton-algorithm'  .AND.    &
+         fft_method /= 'fftw'                 .AND.  fft_method /= 'system-specific' )  THEN
+       message_string = 'unknown fft-algorithm: fft_method = "' // TRIM( fft_method ) // '"'
        CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
     ENDIF
 
-    IF( momentum_advec == 'ws-scheme' .AND.                                    &
-        .NOT. call_psolver_at_all_substeps  ) THEN
-        message_string = 'psolver must be called at each RK3 substep when "'// &
-                      TRIM(momentum_advec) // ' "is used for momentum_advec'
+    IF( momentum_advec == 'ws-scheme' .AND.  .NOT. call_psolver_at_all_substeps  ) THEN
+        message_string = 'psolver must be called at each RK3 substep when "'//                     &
+                         TRIM(momentum_advec) // ' "is used for momentum_advec'
         CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 )
     END IF
 !
 !-- Advection schemes:
-    IF ( momentum_advec /= 'pw-scheme'  .AND.                                  &
-         momentum_advec /= 'ws-scheme'  .AND.                                  &
-         momentum_advec /= 'up-scheme' )                                       &
-    THEN
-       message_string = 'unknown advection scheme: momentum_advec = "' //      &
+    IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme'  .AND.                &
+         momentum_advec /= 'up-scheme' )  THEN
+       message_string = 'unknown advection scheme: momentum_advec = "' //                          &
                         TRIM( momentum_advec ) // '"'
        CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( ( momentum_advec == 'ws-scheme' .OR.  scalar_advec == 'ws-scheme' )   &
-           .AND. ( timestep_scheme == 'euler' .OR.                             &
-                   timestep_scheme == 'runge-kutta-2' ) )                      &
-    THEN
-       message_string = 'momentum_advec or scalar_advec = "'                   &
-         // TRIM( momentum_advec ) // '" is not allowed with ' //              &
-         'timestep_scheme = "' // TRIM( timestep_scheme ) // '"'
+    IF ( ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' )                        &
+         .AND. ( timestep_scheme == 'euler' .OR.  timestep_scheme == 'runge-kutta-2' ) )  THEN
+       message_string = 'momentum_advec or scalar_advec = "' // TRIM( momentum_advec ) //          &
+                        '" is not allowed with timestep_scheme = "' //                             &
+                        TRIM( timestep_scheme ) // '"'
        CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
-         scalar_advec /= 'bc-scheme' .AND. scalar_advec /= 'up-scheme' )       &
-    THEN
-       message_string = 'unknown advection scheme: scalar_advec = "' //        &
-                        TRIM( scalar_advec ) // '"'
+    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme'  .AND.                    &
+         scalar_advec /= 'bc-scheme'  .AND.  scalar_advec /= 'up-scheme' )  THEN
+       message_string = 'unknown advection scheme: scalar_advec = "' // TRIM( scalar_advec ) // '"'
        CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )    &
-    THEN
-       message_string = 'advection_scheme scalar_advec = "'                    &
-         // TRIM( scalar_advec ) // '" not implemented for ' //                &
-         'loop_optimization = "' // TRIM( loop_optimization ) // '"'
+    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )  THEN
+       message_string = 'advection_scheme scalar_advec = "' // TRIM( scalar_advec ) //             &
+                        '" not implemented for loop_optimization = "' //                           &
+                        TRIM( loop_optimization ) // '"'
        CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )
     ENDIF
 
-    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.             &
-         .NOT. use_upstream_for_tke  .AND.                                     &
-         scalar_advec /= 'ws-scheme'                                           &
-       )  THEN
+    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.  .NOT. use_upstream_for_tke     &
+         .AND.  scalar_advec /= 'ws-scheme' )  THEN
        use_upstream_for_tke = .TRUE.
-       message_string = 'use_upstream_for_tke is set to .TRUE. because ' //    &
-                        'use_sgs_for_particles = .TRUE. '          //          &
-                        'and scalar_advec /= ws-scheme'
+       message_string = 'use_upstream_for_tke is set to .TRUE. because ' //                        &
+                        'use_sgs_for_particles = .TRUE. and scalar_advec /= ws-scheme'
        CALL message( 'check_parameters', 'PA0025', 0, 1, 0, 6, 0 )
     ENDIF
@@ -820 +762 @@
 
        CASE DEFAULT
-          message_string = 'unknown timestep scheme: timestep_scheme = "' //   &
+          message_string = 'unknown timestep scheme: timestep_scheme = "' //                       &
                            TRIM( timestep_scheme ) // '"'
           CALL message( 'check_parameters', 'PA0027', 1, 2, 0, 6, 0 )
@@ -826 +768 @@
     END SELECT
 
-    IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme')   &
-         .AND. timestep_scheme(1:5) == 'runge' ) THEN
-       message_string = 'momentum advection scheme "' // &
-                        TRIM( momentum_advec ) // '" & does not work with ' // &
-                        'timestep_scheme "' // TRIM( timestep_scheme ) // '"'
+    IF ( ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' )                     &
+         .AND.  timestep_scheme(1:5) == 'runge' ) THEN
+       message_string = 'momentum advection scheme "' // TRIM( momentum_advec ) //                 &
+                        '" & does not work with timestep_scheme "' // TRIM( timestep_scheme )      &
+                        // '"'
        CALL message( 'check_parameters', 'PA0029', 1, 2, 0, 6, 0 )
     ENDIF
@@ -836 +778 @@
 !-- Check for proper settings for microphysics
     IF ( bulk_cloud_model  .AND.  cloud_droplets )  THEN
-       message_string = 'bulk_cloud_model = .TRUE. is not allowed with ' //    &
-                        'cloud_droplets = .TRUE.'
+       message_string = 'bulk_cloud_model = .TRUE. is not allowed with cloud_droplets = .TRUE.'
        CALL message( 'check_parameters', 'PA0442', 1, 2, 0, 6, 0 )
     ENDIF
@@ -848 +789 @@
     ENDIF
 
-    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
          TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
 !
@@ -857 +798 @@
           SELECT CASE ( action(1:position-1) )
 
-             CASE ( 'set_constant_profiles', 'set_1d-model_profiles',          &
-                    'by_user', 'initialize_vortex', 'initialize_ptanom',       &
-                    'initialize_bubble', 'inifor' )
+             CASE ( 'set_constant_profiles', 'set_1d-model_profiles', 'by_user',                   &
+                    'initialize_vortex', 'initialize_ptanom', 'initialize_bubble', 'inifor' )
                 action = action(position+1:)
 
              CASE DEFAULT
-                message_string = 'initializing_action = "' //                  &
+                message_string = 'initializing_action = "' //                                      &
                                  TRIM( action ) // '" unknown or not allowed'
                 CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
@@ -871 +811 @@
     ENDIF
 
-    IF ( TRIM( initializing_actions ) == 'initialize_vortex'  .AND.            &
-         conserve_volume_flow ) THEN
-         message_string = 'initializing_actions = "initialize_vortex"' //      &
-                        ' is not allowed with conserve_volume_flow = .T.'
+    IF ( TRIM( initializing_actions ) == 'initialize_vortex'  .AND.  conserve_volume_flow ) THEN
+         message_string = 'initializing_actions = "initialize_vortex"' //                          &
+                          ' is not allowed with conserve_volume_flow = .T.'
        CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
     ENDIF
 
 
-    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
+    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.                        &
          INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
-       message_string = 'initializing_actions = "set_constant_profiles"' //    &
-                        ' and "set_1d-model_profiles" are not allowed ' //     &
-                        'simultaneously'
+       message_string = 'initializing_actions = "set_constant_profiles"' //                        &
+                        ' and "set_1d-model_profiles" are not allowed simultaneously'
        CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
     ENDIF
 
-    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
+    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.                        &
          INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
-       message_string = 'initializing_actions = "set_constant_profiles"' //    &
+       message_string = 'initializing_actions = "set_constant_profiles"' //                        &
                         ' and "by_user" are not allowed simultaneously'
        CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
     ENDIF
 
-    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND.                  &
+    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND.                                      &
          INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
-       message_string = 'initializing_actions = "by_user" and ' //             &
+       message_string = 'initializing_actions = "by_user" and ' //                                 &
                         '"set_1d-model_profiles" are not allowed simultaneously'
        CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
     ENDIF
 !
-!-- In case of spinup and nested run, spinup end time must be identical
-!-- in order to have synchronously running simulations.
+!-- In case of spinup and nested run, spinup end time must be identical in order to have
+!-- synchronously running simulations.
     IF ( nested_run )  THEN
 #if defined( __parallel )
-       CALL MPI_ALLREDUCE( spinup_time, spinup_time_max, 1, MPI_REAL,          &
-                           MPI_MAX, MPI_COMM_WORLD, ierr )
-       CALL MPI_ALLREDUCE( dt_spinup,   dt_spinup_max,   1, MPI_REAL,          &
-                           MPI_MAX, MPI_COMM_WORLD, ierr )
-
-       IF ( spinup_time /= spinup_time_max  .OR.  dt_spinup /= dt_spinup_max ) &
-       THEN
-          message_string = 'In case of nesting, spinup_time and ' //           &
-                           'dt_spinup must be identical in all parent ' //     &
-                           'and child domains.'
+       CALL MPI_ALLREDUCE( spinup_time, spinup_time_max, 1, MPI_REAL, MPI_MAX, MPI_COMM_WORLD,     &
+                           ierr )
+       CALL MPI_ALLREDUCE( dt_spinup,   dt_spinup_max,   1, MPI_REAL, MPI_MAX, MPI_COMM_WORLD,     &
+                           ierr )
+
+       IF ( spinup_time /= spinup_time_max  .OR.  dt_spinup /= dt_spinup_max )  THEN
+          message_string = 'In case of nesting, spinup_time and ' //                               &
+                           'dt_spinup must be identical in all parent and child domains.'
           CALL message( 'check_parameters', 'PA0489', 3, 2, 0, 6, 0 )
        ENDIF
@@ -920 +856 @@
     ENDIF
 
-    IF ( bulk_cloud_model  .AND.  .NOT.  humidity )  THEN
-       WRITE( message_string, * ) 'bulk_cloud_model = ', bulk_cloud_model,     &
+    IF ( bulk_cloud_model  .AND.  .NOT. humidity )  THEN
+       WRITE( message_string, * ) 'bulk_cloud_model = ', bulk_cloud_model,                         &
               ' is not allowed with humidity = ', humidity
        CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
@@ -927 +863 @@
 
     IF ( humidity  .AND.  sloping_surface )  THEN
-       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' //   &
+       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' //                       &
                         'are not allowed simultaneously'
        CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
@@ -936 +872 @@
 
 !
-!-- In case of no restart run, check initialising parameters and calculate
-!-- further quantities
+!-- In case of no restart run, check initialising parameters and calculate further quantities
     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
 
@@ -947 +882 @@
 
 !--
-!--    If required, compute initial profile of the geostrophic wind
-!--    (component ug)
+!--    If required, compute initial profile of the geostrophic wind (component ug)
        i = 1
        gradient = 0.0_wp
@@ -958 +892 @@
           DO  k = 1, nzt+1
              IF ( i < 11 )  THEN
-                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND.              &
+                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND.                                  &
                      ug_vertical_gradient_level(i) >= 0.0_wp )  THEN
                    gradient = ug_vertical_gradient(i) / 100.0_wp
@@ -982 +916 @@
           DO  k = nzt, nzb, -1
              IF ( i < 11 )  THEN
-                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND.              &
+                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND.                                  &
                      ug_vertical_gradient_level(i) <= 0.0_wp )  THEN
                    gradient = ug_vertical_gradient(i) / 100.0_wp
@@ -1011 +945 @@
 !
 !--
-!--    If required, compute initial profile of the geostrophic wind
-!--    (component vg)
+!--    If required, compute initial profile of the geostrophic wind (component vg)
        i = 1
        gradient = 0.0_wp
@@ -1022 +955 @@
           DO  k = 1, nzt+1
              IF ( i < 11 )  THEN
-                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND.              &
+                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND.                                  &
                      vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
                    gradient = vg_vertical_gradient(i) / 100.0_wp
@@ -1046 +979 @@
           DO  k = nzt, nzb, -1
              IF ( i < 11 )  THEN
-                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND.              &
+                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND.                                  &
                      vg_vertical_gradient_level(i) <= 0.0_wp )  THEN
                    gradient = vg_vertical_gradient(i) / 100.0_wp
@@ -1074 +1007 @@
 
 !
-!--    Let the initial wind profiles be the calculated ug/vg profiles or
-!--    interpolate them from wind profile data (if given)
+!--    Let the initial wind profiles be the calculated ug/vg profiles or interpolate them from wind
+!--    profile data (if given)
        IF ( u_profile(1) == 9999999.9_wp  .AND.  v_profile(1) == 9999999.9_wp )  THEN
 
@@ -1089 +1022 @@
 
           IF ( omega /= 0.0_wp )  THEN
-             message_string = 'Coriolis force must be switched off (by setting omega=0.0)' //  &
+             message_string = 'Coriolis force must be switched off (by setting omega=0.0)' //      &
                               ' when prescribing the forcing by u_profile and v_profile'
              CALL message( 'check_parameters', 'PA0347', 1, 2, 0, 6, 0 )
@@ -1110 +1043 @@
 
              IF ( kk < 200  .AND.  uv_heights(kk+1) /= 9999999.9_wp )  THEN
-                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
-                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
-                                       ( u_profile(kk+1) - u_profile(kk) )
-                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
-                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
-                                       ( v_profile(kk+1) - v_profile(kk) )
+                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /                           &
+                                            ( uv_heights(kk+1) - uv_heights(kk) ) *                &
+                                            ( u_profile(kk+1) - u_profile(kk) )
+                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /                           &
+                                            ( uv_heights(kk+1) - uv_heights(kk) ) *                &
+                                            ( v_profile(kk+1) - v_profile(kk) )
              ELSE
                 u_init(k) = u_profile(kk)
@@ -1133 +1066 @@
 !--    Compute initial temperature profile using the given temperature gradients
        IF (  .NOT.  neutral )  THEN
-          CALL init_vertical_profiles( pt_vertical_gradient_level_ind,          &
-                                       pt_vertical_gradient_level,              &
-                                       pt_vertical_gradient, pt_init,           &
-                                       pt_surface, bc_pt_t_val )
+          CALL init_vertical_profiles( pt_vertical_gradient_level_ind, pt_vertical_gradient_level, &
+                                       pt_vertical_gradient, pt_init, pt_surface, bc_pt_t_val )
        ENDIF
 !
 !--    Compute initial humidity profile using the given humidity gradients
        IF ( humidity )  THEN
-          CALL init_vertical_profiles( q_vertical_gradient_level_ind,          &
-                                       q_vertical_gradient_level,              &
-                                       q_vertical_gradient, q_init,            &
-                                       q_surface, bc_q_t_val )
+          CALL init_vertical_profiles( q_vertical_gradient_level_ind, q_vertical_gradient_level,   &
+                                       q_vertical_gradient, q_init, q_surface, bc_q_t_val )
        ENDIF
 !
 !--    Compute initial scalar profile using the given scalar gradients
        IF ( passive_scalar )  THEN
-          CALL init_vertical_profiles( s_vertical_gradient_level_ind,          &
-                                       s_vertical_gradient_level,              &
-                                       s_vertical_gradient, s_init,            &
-                                       s_surface, bc_s_t_val )
+          CALL init_vertical_profiles( s_vertical_gradient_level_ind, s_vertical_gradient_level,   &
+                                       s_vertical_gradient, s_init, s_surface, bc_s_t_val )
        ENDIF
 !
@@ -1164 +1091 @@
 !-- Check if the control parameter use_subsidence_tendencies is used correctly
     IF ( use_subsidence_tendencies  .AND.  .NOT.  large_scale_subsidence )  THEN
-       message_string = 'The usage of use_subsidence_tendencies ' //           &
-                            'requires large_scale_subsidence = .T..'
+       message_string = 'The usage of use_subsidence_tendencies ' //                               &
+                        'requires large_scale_subsidence = .T..'
        CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
     ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
        message_string = 'The usage of use_subsidence_tendencies ' //           &
-                            'requires large_scale_forcing = .T..'
+                        'requires large_scale_forcing = .T..'
        CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1176 +1103 @@
 !-- Initialize large scale subsidence if required
     If ( large_scale_subsidence )  THEN
-       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp  .AND.            &
-                                     .NOT.  large_scale_forcing )  THEN
+       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp  .AND. .NOT. large_scale_forcing )    &
+       THEN
           CALL init_w_subsidence
        ENDIF
 !
-!--    In case large_scale_forcing is used, profiles for subsidence velocity
-!--    are read in from file LSF_DATA
-
-       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp  .AND.            &
-            .NOT.  large_scale_forcing )  THEN
-          message_string = 'There is no default large scale vertical ' //      &
-                           'velocity profile set. Specify the subsidence ' //  &
-                           'velocity profile via subs_vertical_gradient ' //   &
-                           'and subs_vertical_gradient_level.'
+!--    In case large_scale_forcing is used, profiles for subsidence velocity are read in from file
+!--    LSF_DATA
+
+       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp  .AND. .NOT. large_scale_forcing )    &
+       THEN
+          message_string = 'There is no default large scale vertical velocity profile set. ' //    &
+                           'Specify the subsidence velocity profile via subs_vertical_gradient' // &
+                           ' and subs_vertical_gradient_level.'
           CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 )
        ENDIF
     ELSE
         IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
-           message_string = 'Enable usage of large scale subsidence by ' //    &
+           message_string = 'Enable usage of large scale subsidence by ' //                        &
                             'setting large_scale_subsidence = .T..'
           CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
@@ -1223 +1149 @@
        vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
     ELSE
-       message_string = 'illegal value for reference_state: "' //              &
-                        TRIM( reference_state ) // '"'
+       message_string = 'illegal value for reference_state: "' // TRIM( reference_state ) // '"'
        CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1232 +1157 @@
     IF ( alpha_surface /= 0.0_wp )  THEN
        IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
-          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface,   &
-                                     ' ) must be < 90.0'
+          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface, ' ) must be < 90.0'
           CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
        ENDIF
@@ -1262 +1186 @@
           ENDIF
        ELSE
-          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor,              &
+          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor,                                  &
                  ' out of range &0.0 < cfl_factor <= 1.0 is required'
           CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
@@ -1273 +1197 @@
 
 !
-!-- Store reference time for coupled runs and change the coupling flag,
-!-- if ...
+!-- Store reference time for coupled runs and change the coupling flag, if ...
     IF ( simulated_time == 0.0_wp )  THEN
        IF ( coupling_start_time == 0.0_wp )  THEN
@@ -1286 +1209 @@
 !-- Set wind speed in the Galilei-transformed system
     IF ( galilei_transformation )  THEN
-       IF ( use_ug_for_galilei_tr                    .AND.                     &
-            ug_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
-            ug_vertical_gradient(1) == 0.0_wp        .AND.                     &
-            vg_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
-            vg_vertical_gradient(1) == 0.0_wp )  THEN
+       IF ( use_ug_for_galilei_tr                    .AND.                                         &
+            ug_vertical_gradient_level(1) == 0.0_wp  .AND.                                         &
+            ug_vertical_gradient(1)       == 0.0_wp  .AND.                                         &
+            vg_vertical_gradient_level(1) == 0.0_wp  .AND.                                         &
+            vg_vertical_gradient(1)       == 0.0_wp )  THEN
           u_gtrans = ug_surface * 0.6_wp
           v_gtrans = vg_surface * 0.6_wp
-       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
-                ( ug_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
-                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
-          message_string = 'baroclinity (ug) not allowed simultaneously' //    &
+       ELSEIF ( use_ug_for_galilei_tr  .AND.  ( ug_vertical_gradient_level(1) /= 0.0_wp .OR.       &
+                                                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
+          message_string = 'baroclinity (ug) not allowed simultaneously' //                        &
                            ' with galilei transformation'
           CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
-       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
-                ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
-                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
-          message_string = 'baroclinity (vg) not allowed simultaneously' //    &
+       ELSEIF ( use_ug_for_galilei_tr  .AND.  ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.      &
+                                                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
+          message_string = 'baroclinity (vg) not allowed simultaneously' //                        &
                            ' with galilei transformation'
           CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
        ELSE
-          message_string = 'variable translation speed used for Galilei-' //   &
-             'transformation, which may cause & instabilities in stably ' //   &
-             'stratified regions'
+          message_string = 'variable translation speed used for Galilei-transformation, which ' // &
+                           'may cause & instabilities in stably stratified regions'
           CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
        ENDIF
@@ -1314 +1234 @@
 
 !
-!-- In case of using a prandtl-layer, calculated (or prescribed) surface
-!-- fluxes have to be used in the diffusion-terms
+!-- In case of using a prandtl-layer, calculated (or prescribed) surface fluxes have to be used in
+!-- the diffusion-terms
     IF ( constant_flux_layer )  use_surface_fluxes = .TRUE.
 
 !
 !-- Check boundary conditions and set internal variables:
-!-- Attention: the lateral boundary conditions have been already checked in
-!-- parin
-!
-!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
-!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
-!-- and tools do not work with non-cyclic boundary conditions.
+!-- Attention: the lateral boundary conditions have been already checked in parin
+!
+!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-Willimas or
+!-- Wicker - Skamarock advection scheme. Several schemes and tools do not work with non-cyclic
+!-- boundary conditions.
     IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
        IF ( psolver(1:9) /= 'multigrid' )  THEN
-          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
-                           'psolver = "' // TRIM( psolver ) // '"'
+          message_string = 'non-cyclic lateral boundaries do not allow ' // 'psolver = "' //       &
+                           TRIM( psolver ) // '"'
           CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
        ENDIF
-       IF ( momentum_advec /= 'pw-scheme'  .AND.                               &
-            momentum_advec /= 'ws-scheme' )  THEN
-
-          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
-                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
+       IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' )  THEN
+
+          message_string = 'non-cyclic lateral boundaries do not allow momentum_advec = "' //      &
+                           TRIM( momentum_advec ) // '"'
           CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
        ENDIF
-       IF ( scalar_advec /= 'pw-scheme'  .AND.                                 &
-            scalar_advec /= 'ws-scheme' )  THEN
-          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
-                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
+       IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' )  THEN
+          message_string = 'non-cyclic lateral boundaries do not allow scalar_advec = "' //        &
+                           TRIM( scalar_advec ) // '"'
           CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
        ENDIF
        IF ( galilei_transformation )  THEN
-          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
-                           'galilei_transformation = .T.'
+          message_string = 'non-cyclic lateral boundaries do not allow galilei_transformation = .T.'
           CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
        ENDIF
@@ -1361 +1277 @@
           bc_e_b = 'neumann'
           ibc_e_b = 1
-          message_string = 'boundary condition bc_e_b changed to "' //         &
-                           TRIM( bc_e_b ) // '"'
+          message_string = 'boundary condition bc_e_b changed to "' // TRIM( bc_e_b ) // '"'
           CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
        ENDIF
     ELSE
-       message_string = 'unknown boundary condition: bc_e_b = "' //            &
-                        TRIM( bc_e_b ) // '"'
+       message_string = 'unknown boundary condition: bc_e_b = "' // TRIM( bc_e_b ) // '"'
        CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1378 +1292 @@
        ibc_p_b = 1
     ELSE
-       message_string = 'unknown boundary condition: bc_p_b = "' //            &
-                        TRIM( bc_p_b ) // '"'
+       message_string = 'unknown boundary condition: bc_p_b = "' // TRIM( bc_p_b ) // '"'
        CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1389 +1302 @@
        ibc_p_t = 1
     ELSE
-       message_string = 'unknown boundary condition: bc_p_t = "' //            &
-                        TRIM( bc_p_t ) // '"'
+       message_string = 'unknown boundary condition: bc_p_t = "' // TRIM( bc_p_t ) // '"'
        CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1404 +1316 @@
           ibc_pt_b = 1
        ELSE
-          message_string = 'unknown boundary condition: bc_pt_b = "' //        &
-                           TRIM( bc_pt_b ) // '"'
+          message_string = 'unknown boundary condition: bc_pt_b = "' // TRIM( bc_pt_b ) // '"'
           CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
        ENDIF
@@ -1419 +1330 @@
        ibc_pt_t = 3
     ELSE
-       message_string = 'unknown boundary condition: bc_pt_t = "' //           &
-                        TRIM( bc_pt_t ) // '"'
+       message_string = 'unknown boundary condition: bc_pt_t = "' // TRIM( bc_pt_t ) // '"'
        CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
     ENDIF
 
-    IF ( ANY( wall_heatflux /= 0.0_wp )  .AND.                        &
-         surface_heatflux == 9999999.9_wp )  THEN
-       message_string = 'wall_heatflux additionally requires ' //     &
-                        'setting of surface_heatflux'
+    IF ( ANY( wall_heatflux /= 0.0_wp )  .AND.  surface_heatflux == 9999999.9_wp )  THEN
+       message_string = 'wall_heatflux additionally requires setting of surface_heatflux'
        CALL message( 'check_parameters', 'PA0443', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1451 +1359 @@
     IF ( neutral )  THEN
 
-       IF ( surface_heatflux /= 0.0_wp  .AND.                                  &
-            surface_heatflux /= 9999999.9_wp )  THEN
+       IF ( surface_heatflux /= 0.0_wp  .AND.  surface_heatflux /= 9999999.9_wp )  THEN
           message_string = 'heatflux must not be set for pure neutral flow'
           CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
        ENDIF
 
-       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp )      &
-       THEN
+       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp )  THEN
           message_string = 'heatflux must not be set for pure neutral flow'
           CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
@@ -1465 +1371 @@
     ENDIF
 
-    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.                             &
-         top_momentumflux_v /= 9999999.9_wp )  THEN
+    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.  top_momentumflux_v /= 9999999.9_wp )  THEN
        constant_top_momentumflux = .TRUE.
-    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.                &
+    ELSEIF ( .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.                                     &
            top_momentumflux_v == 9999999.9_wp ) )  THEN
-       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' //  &
-                        'must be set'
+       message_string = 'both, top_momentumflux_u AND top_momentumflux_v must be set'
        CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
     ENDIF
 
 !
-!-- A given surface temperature implies Dirichlet boundary condition for
-!-- temperature. In this case specification of a constant heat flux is
-!-- forbidden.
-    IF ( ibc_pt_b == 0  .AND.  constant_heatflux  .AND.                        &
-         surface_heatflux /= 0.0_wp )  THEN
-       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
+!-- A given surface temperature implies Dirichlet boundary condition for temperature. In this case
+!-- specification of a constant heat flux is forbidden.
+    IF ( ibc_pt_b == 0  .AND.  constant_heatflux  .AND.  surface_heatflux /= 0.0_wp )  THEN
+       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //                    &
                         '& is not allowed with constant_heatflux = .TRUE.'
        CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
     ENDIF
     IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0_wp )  THEN
-       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo',  &
-               'wed with pt_surface_initial_change (/=0) = ',                  &
-               pt_surface_initial_change
+       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo',                      &
+               'wed with pt_surface_initial_change (/=0) = ', pt_surface_initial_change
        CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
     ENDIF
 
 !
-!-- A given temperature at the top implies Dirichlet boundary condition for
-!-- temperature. In this case specification of a constant heat flux is
-!-- forbidden.
-    IF ( ibc_pt_t == 0  .AND.  constant_top_heatflux  .AND.                    &
-         top_heatflux /= 0.0_wp )  THEN
-       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
+!-- A given temperature at the top implies Dirichlet boundary condition for temperature. In this
+!-- case specification of a constant heat flux is forbidden.
+    IF ( ibc_pt_t == 0  .AND.  constant_top_heatflux  .AND.  top_heatflux /= 0.0_wp )  THEN
+       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //                    &
                         '" is not allowed with constant_top_heatflux = .TRUE.'
        CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
@@ -1507 +1406 @@
     IF ( humidity )  THEN
 
-       IF ( ANY( wall_humidityflux /= 0.0_wp )  .AND.                        &
-            surface_waterflux == 9999999.9_wp )  THEN
-          message_string = 'wall_humidityflux additionally requires ' //     &
-                           'setting of surface_waterflux'
+       IF ( ANY( wall_humidityflux /= 0.0_wp )  .AND.  surface_waterflux == 9999999.9_wp )  THEN
+          message_string = 'wall_humidityflux additionally requires setting of surface_waterflux'
           CALL message( 'check_parameters', 'PA0444', 1, 2, 0, 6, 0 )
        ENDIF
 
-       CALL set_bc_scalars( 'q', bc_q_b, bc_q_t, ibc_q_b, ibc_q_t,           &
-                            'PA0071', 'PA0072' )
+       CALL set_bc_scalars( 'q', bc_q_b, bc_q_t, ibc_q_b, ibc_q_t, 'PA0071', 'PA0072' )
 
        IF ( surface_waterflux == 9999999.9_wp  )  THEN
@@ -1530 +1426 @@
        ENDIF
 
-       CALL check_bc_scalars( 'q', bc_q_b, ibc_q_b, 'PA0073', 'PA0074',        &
-                              constant_waterflux, q_surface_initial_change )
+       CALL check_bc_scalars( 'q', bc_q_b, ibc_q_b, 'PA0073', 'PA0074', constant_waterflux,        &
+                              q_surface_initial_change )
 
     ENDIF
@@ -1537 +1433 @@
     IF ( passive_scalar )  THEN
 
-       IF ( ANY( wall_scalarflux /= 0.0_wp )  .AND.                            &
-            surface_scalarflux == 9999999.9_wp )  THEN
-          message_string = 'wall_scalarflux additionally requires ' //         &
-                           'setting of surface_scalarflux'
+       IF ( ANY( wall_scalarflux /= 0.0_wp )  .AND.  surface_scalarflux == 9999999.9_wp )  THEN
+          message_string = 'wall_scalarflux additionally requires setting of surface_scalarflux'
           CALL message( 'check_parameters', 'PA0445', 1, 2, 0, 6, 0 )
        ENDIF
@@ -1546 +1440 @@
        IF ( surface_scalarflux == 9999999.9_wp )  constant_scalarflux = .FALSE.
 
-       CALL set_bc_scalars( 's', bc_s_b, bc_s_t, ibc_s_b, ibc_s_t,             &
-                            'PA0071', 'PA0072' )
-
-       CALL check_bc_scalars( 's', bc_s_b, ibc_s_b, 'PA0073', 'PA0074',        &
-                              constant_scalarflux, s_surface_initial_change )
+       CALL set_bc_scalars( 's', bc_s_b, bc_s_t, ibc_s_b, ibc_s_t, 'PA0071', 'PA0072' )
+
+       CALL check_bc_scalars( 's', bc_s_b, ibc_s_b, 'PA0073', 'PA0074', constant_scalarflux,       &
+                              s_surface_initial_change )
 
        IF ( top_scalarflux == 9999999.9_wp )  constant_top_scalarflux = .FALSE.
 !
-!--    A fixed scalar concentration at the top implies Dirichlet boundary
-!--    condition for scalar. Hence, in this case specification of a constant
-!--    scalar flux is forbidden.
-       IF ( ( ibc_s_t == 0 .OR. ibc_s_t == 2 )  .AND.  constant_top_scalarflux &
-               .AND.  top_scalarflux /= 0.0_wp )  THEN
-          message_string = 'boundary condition: bc_s_t = "' //                 &
-                           TRIM( bc_s_t ) // '" is not allowed with ' //       &
-                           'top_scalarflux /= 0.0'
+!--    A fixed scalar concentration at the top implies Dirichlet boundary condition for scalar.
+!--    Hence, in this case specification of a constant scalar flux is forbidden.
+       IF ( ( ibc_s_t == 0 .OR. ibc_s_t == 2 )  .AND.  constant_top_scalarflux  .AND.              &
+              top_scalarflux /= 0.0_wp )  THEN
+          message_string = 'boundary condition: bc_s_t = "' // TRIM( bc_s_t ) //                   &
+                           '" is not allowed with top_scalarflux /= 0.0'
           CALL message( 'check_parameters', 'PA0441', 1, 2, 0, 6, 0 )
        ENDIF
@@ -1573 +1464 @@
        ibc_uv_b = 1
        IF ( constant_flux_layer )  THEN
-          message_string = 'boundary condition: bc_uv_b = "' //                &
-               TRIM( bc_uv_b ) // '" is not allowed with constant_flux_layer'  &
-               // ' = .TRUE.'
+          message_string = 'boundary condition: bc_uv_b = "' // TRIM( bc_uv_b ) //                 &
+                           '" is not allowed with constant_flux_layer = .TRUE.'
           CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
        ENDIF
     ELSE
-       message_string = 'unknown boundary condition: bc_uv_b = "' //           &
-                        TRIM( bc_uv_b ) // '"'
+       message_string = 'unknown boundary condition: bc_uv_b = "' // TRIM( bc_uv_b ) // '"'
        CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
     ENDIF
 !
-!-- In case of coupled simulations u and v at the ground in atmosphere will be
-!-- assigned with the u and v values of the ocean surface
+!-- In case of coupled simulations u and v at the ground in atmosphere will be assigned with the u
+!-- and v values of the ocean surface
     IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
        ibc_uv_b = 2
@@ -1598 +1487 @@
           IF ( bc_uv_t == 'dirichlet_0' )  THEN
 !
-!--          Velocities for the initial u,v-profiles are set zero at the top
-!--          in case of dirichlet_0 conditions
+!--          Velocities for the initial u,v-profiles are set zero at the top in case of dirichlet_0
+!--          conditions
              u_init(nzt+1)    = 0.0_wp
              v_init(nzt+1)    = 0.0_wp
@@ -1608 +1497 @@
           ibc_uv_t = 3
        ELSE
-          message_string = 'unknown boundary condition: bc_uv_t = "' //        &
-                           TRIM( bc_uv_t ) // '"'
+          message_string = 'unknown boundary condition: bc_uv_t = "' // TRIM( bc_uv_t ) // '"'
           CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
        ENDIF
@@ -1619 +1507 @@
        rayleigh_damping_factor = 0.0_wp
     ELSE
-       IF ( rayleigh_damping_factor < 0.0_wp  .OR.                             &
-            rayleigh_damping_factor > 1.0_wp )  THEN
-          WRITE( message_string, * )  'rayleigh_damping_factor = ',            &
-                              rayleigh_damping_factor, ' out of range [0.0,1.0]'
+       IF ( rayleigh_damping_factor < 0.0_wp  .OR.  rayleigh_damping_factor > 1.0_wp )  THEN
+          WRITE( message_string, * )  'rayleigh_damping_factor = ', rayleigh_damping_factor,       &
+                 ' out of range [0.0,1.0]'
           CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
        ENDIF
@@ -1628 +1515 @@
 
     IF ( rayleigh_damping_height == -1.0_wp )  THEN
-       IF (  .NOT.  ocean_mode )  THEN
+       IF ( .NOT. ocean_mode )  THEN
           rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
        ELSE
@@ -1634 +1521 @@
        ENDIF
     ELSE
-       IF (  .NOT.  ocean_mode )  THEN
-          IF ( rayleigh_damping_height < 0.0_wp  .OR.                          &
-               rayleigh_damping_height > zu(nzt) )  THEN
-             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
-                   rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
+       IF ( .NOT. ocean_mode )  THEN
+          IF ( rayleigh_damping_height < 0.0_wp  .OR.  rayleigh_damping_height > zu(nzt) )  THEN
+             WRITE( message_string, * )  'rayleigh_damping_height = ',  rayleigh_damping_height,   &
+                    ' out of range [0.0,', zu(nzt), ']'
              CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
           ENDIF
        ELSE
-          IF ( rayleigh_damping_height > 0.0_wp  .OR.                          &
-               rayleigh_damping_height < zu(nzb) )  THEN
-             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
-                   rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
+          IF ( rayleigh_damping_height > 0.0_wp  .OR.  rayleigh_damping_height < zu(nzb) )  THEN
+             WRITE( message_string, * )  'rayleigh_damping_height = ', rayleigh_damping_height,    &
+                    ' out of range [0.0,', zu(nzb), ']'
              CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
           ENDIF
@@ -1654 +1539 @@
 !-- Check number of chosen statistic regions
     IF ( statistic_regions < 0 )  THEN
-       WRITE ( message_string, * ) 'number of statistic_regions = ',           &
-                   statistic_regions+1, ' is not allowed'
+       WRITE ( message_string, * ) 'number of statistic_regions = ', statistic_regions+1,          &
+               ' is not allowed'
        CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( normalizing_region > statistic_regions  .OR.                          &
-         normalizing_region < 0)  THEN
-       WRITE ( message_string, * ) 'normalizing_region = ',                    &
-                normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
-                ' (value of statistic_regions)'
+    IF ( normalizing_region > statistic_regions  .OR.  normalizing_region < 0)  THEN
+       WRITE ( message_string, * ) 'normalizing_region = ', normalizing_region,                    &
+               ' must be >= 0 and <= ',statistic_regions, ' (value of statistic_regions)'
        CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1684 +1567 @@
 !
 !-- Set the default skip time intervals for data output, if necessary
-    IF ( skip_time_dopr    == 9999999.9_wp )                                   &
-                                       skip_time_dopr    = skip_time_data_output
-    IF ( skip_time_do2d_xy == 9999999.9_wp )                                   &
-                                       skip_time_do2d_xy = skip_time_data_output
-    IF ( skip_time_do2d_xz == 9999999.9_wp )                                   &
-                                       skip_time_do2d_xz = skip_time_data_output
-    IF ( skip_time_do2d_yz == 9999999.9_wp )                                   &
-                                       skip_time_do2d_yz = skip_time_data_output
-    IF ( skip_time_do3d    == 9999999.9_wp )                                   &
-                                       skip_time_do3d    = skip_time_data_output
-    IF ( skip_time_data_output_av == 9999999.9_wp )                            &
-                                skip_time_data_output_av = skip_time_data_output
+    IF ( skip_time_dopr    == 9999999.9_wp )  skip_time_dopr    = skip_time_data_output
+    IF ( skip_time_do2d_xy == 9999999.9_wp )  skip_time_do2d_xy = skip_time_data_output
+    IF ( skip_time_do2d_xz == 9999999.9_wp )  skip_time_do2d_xz = skip_time_data_output
+    IF ( skip_time_do2d_yz == 9999999.9_wp )  skip_time_do2d_yz = skip_time_data_output
+    IF ( skip_time_do3d    == 9999999.9_wp )  skip_time_do3d    = skip_time_data_output
+    IF ( skip_time_data_output_av == 9999999.9_wp )                                                &
+                                       skip_time_data_output_av = skip_time_data_output
     DO  mid = 1, max_masks
-       IF ( skip_time_domask(mid) == 9999999.9_wp )                            &
-                                skip_time_domask(mid)    = skip_time_data_output
+       IF ( skip_time_domask(mid) == 9999999.9_wp )                                                &
+                                       skip_time_domask(mid)    = skip_time_data_output
     ENDDO
 
@@ -1704 +1582 @@
 !-- Check the average intervals (first for 3d-data, then for profiles)
     IF ( averaging_interval > dt_data_output_av )  THEN
-       WRITE( message_string, * )  'averaging_interval = ',                    &
-             averaging_interval, ' must be <= dt_data_output_av = ',           &
-             dt_data_output_av
+       WRITE( message_string, * )  'averaging_interval = ', averaging_interval,                    &
+              ' must be <= dt_data_output_av = ', dt_data_output_av
        CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1715 +1592 @@
 
     IF ( averaging_interval_pr > dt_dopr )  THEN
-       WRITE( message_string, * )  'averaging_interval_pr = ',                 &
-             averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
+       WRITE( message_string, * )  'averaging_interval_pr = ', averaging_interval_pr,              &
+              ' must be <= dt_dopr = ', dt_dopr
        CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
     ENDIF
@@ -1727 +1604 @@
 
 !
-!-- Set the default interval for the output of timeseries to a reasonable
-!-- value (tries to minimize the number of calls of flow_statistics)
+!-- Set the default interval for the output of timeseries to a reasonable value (tries to minimize
+!-- the number of calls of flow_statistics)
     IF ( dt_dots == 9999999.9_wp )  THEN
        IF ( averaging_interval_pr == 0.0_wp )  THEN
@@ -1740 +1617 @@
 !-- Check the sample rate for averaging (first for 3d-data, then for profiles)
     IF ( dt_averaging_input > averaging_interval )  THEN
-       WRITE( message_string, * )  'dt_averaging_input = ',                    &
-                dt_averaging_input, ' must be <= averaging_interval = ',       &
-                averaging_interval
+       WRITE( message_string, * )  'dt_averaging_input = ', dt_averaging_input,                    &
+              ' must be <= averaging_interval = ', averaging_interval
        CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
     ENDIF
 
     IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
-       WRITE( message_string, * )  'dt_averaging_input_pr = ',                 &
-                dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
-                averaging_interval_pr
+       WRITE( message_string, * )  'dt_averaging_input_pr = ', dt_averaging_input_pr,              &
+              ' must be <= averaging_interval_pr = ', averaging_interval_pr
        CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
     ENDIF
 
 !
-!-- Determine the number of output profiles and check whether they are
-!-- permissible
+!-- Determine the number of output profiles and check whether they are permissible
     DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
 
@@ -1762 +1636 @@
 
 !
-!--    Determine internal profile number (for hom, homs)
-!--    and store height levels
+!--    Determine internal profile number (for hom, homs) and store height levels
        SELECT CASE ( TRIM( data_output_pr(i) ) )
 
@@ -1792 +1665 @@
 
           CASE ( 'theta', '#theta' )
-             IF ( .NOT. bulk_cloud_model ) THEN
+             IF ( .NOT. bulk_cloud_model )  THEN
                 dopr_index(i) = 4
                 dopr_unit(i)  = 'K'
@@ -1974 +1847 @@
           CASE ( 'q', '#q' )
              IF ( .NOT. humidity )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for humidity = .FALSE.'
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for humidity = .FALSE.'
                 CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
              ELSE
@@ -1992 +1864 @@
           CASE ( 's', '#s' )
              IF ( .NOT. passive_scalar )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for passive_scalar = .FALSE.'
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for passive_scalar = .FALSE.'
                 CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
              ELSE
@@ -2033 +1904 @@
           CASE ( 'thetal', '#thetal' )
              IF ( .NOT. bulk_cloud_model ) THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for bulk_cloud_model = .FALSE.'
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for bulk_cloud_model = .FALSE.'
                 CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
              ELSE
@@ -2076 +1946 @@
 
           CASE ( 'w"q"' )
-             IF (  .NOT.  humidity )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for humidity = .FALSE.'
+             IF ( .NOT. humidity )  THEN
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for humidity = .FALSE.'
                 CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
              ELSE
@@ -2088 +1957 @@
 
           CASE ( 'w*q*' )
-             IF (  .NOT.  humidity )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for humidity = .FALSE.'
+             IF ( .NOT. humidity )  THEN
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for humidity = .FALSE.'
                 CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
              ELSE
@@ -2100 +1968 @@
 
           CASE ( 'wq' )
-             IF (  .NOT.  humidity )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for humidity = .FALSE.'
+             IF ( .NOT. humidity )  THEN
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for humidity = .FALSE.'
                 CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
              ELSE
@@ -2112 +1979 @@
 
           CASE ( 'w"s"' )
-             IF (  .NOT.  passive_scalar )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for passive_scalar = .FALSE.'
+             IF ( .NOT. passive_scalar )  THEN
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for passive_scalar = .FALSE.'
                 CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
              ELSE
@@ -2124 +1990 @@
 
           CASE ( 'w*s*' )
-             IF (  .NOT.  passive_scalar )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for passive_scalar = .FALSE.'
+             IF ( .NOT. passive_scalar )  THEN
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for passive_scalar = .FALSE.'
                 CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
              ELSE
@@ -2136 +2001 @@
 
           CASE ( 'ws' )
-             IF (  .NOT.  passive_scalar )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for passive_scalar = .FALSE.'
+             IF ( .NOT. passive_scalar )  THEN
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for passive_scalar = .FALSE.'
                 CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
              ELSE
@@ -2148 +2012 @@
 
           CASE ( 'w"qv"' )
-             IF ( humidity  .AND.  .NOT.  bulk_cloud_model )  THEN
+             IF ( humidity  .AND.  .NOT. bulk_cloud_model )  THEN
                 dopr_index(i) = 48
                 dopr_unit(i)  = TRIM ( waterflux_output_unit )
@@ -2157 +2021 @@
                 hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for bulk_cloud_model = .FALSE. ' // &
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for bulk_cloud_model = .FALSE. ' //          &
                                  'and humidity = .FALSE.'
                 CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
@@ -2169 +2032 @@
                 dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
-             ELSEIF( humidity .AND. bulk_cloud_model ) THEN
+             ELSEIF( humidity  .AND.  bulk_cloud_model )  THEN
                 dopr_index(i) = 52
                 dopr_unit(i)  = TRIM ( waterflux_output_unit )
                 hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for bulk_cloud_model = .FALSE. ' // &
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for bulk_cloud_model = .FALSE. ' //          &
                                  'and humidity = .FALSE.'
                 CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
@@ -2182 +2044 @@
 
           CASE ( 'wqv' )
-             IF ( humidity  .AND.  .NOT.  bulk_cloud_model )  THEN
+             IF ( humidity  .AND.  .NOT. bulk_cloud_model )  THEN
                 dopr_index(i) = 50
                 dopr_unit(i)  = TRIM ( waterflux_output_unit )
@@ -2191 +2053 @@
                 hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for bulk_cloud_model = .FALSE. ' // &
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for bulk_cloud_model = .FALSE. ' //          &
                                  'and humidity = .FALSE.'
                 CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
@@ -2199 +2060 @@
 
           CASE ( 'ql' )
-             IF (  .NOT.  bulk_cloud_model  .AND.  .NOT.  cloud_droplets )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for bulk_cloud_model = .FALSE. ' // &
+             IF ( .NOT. bulk_cloud_model  .AND.  .NOT. cloud_droplets )  THEN
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for bulk_cloud_model = .FALSE. ' //          &
                                  'and cloud_droplets = .FALSE.'
                 CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
@@ -2268 +2128 @@
           CASE ( 'q*2' )
              IF (  .NOT.  humidity )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for humidity = .FALSE.'
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for humidity = .FALSE.'
                 CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
              ELSE
@@ -2305 +2164 @@
           CASE ( 'w_subs' )
              IF (  .NOT.  large_scale_subsidence )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for large_scale_subsidence = .FALSE.'
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for large_scale_subsidence = .FALSE.'
                 CALL message( 'check_parameters', 'PA0382', 1, 2, 0, 6, 0 )
              ELSE
@@ -2317 +2175 @@
           CASE ( 's*2' )
              IF (  .NOT.  passive_scalar )  THEN
-                message_string = 'data_output_pr = ' //                        &
-                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
-                                 'lemented for passive_scalar = .FALSE.'
+                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
+                                 ' is not implemented for passive_scalar = .FALSE.'
                 CALL message( 'check_parameters', 'PA0185', 1, 2, 0, 6, 0 )
              ELSE
@@ -2338 +2195 @@
 !
 !--          Check for other modules
-             CALL module_interface_check_data_output_pr( data_output_pr(i), i, &
-                                                         unit, dopr_unit(i) )
+             CALL module_interface_check_data_output_pr( data_output_pr(i), i, unit, dopr_unit(i) )
 
 !
@@ -2345 +2201 @@
              IF ( unit == 'illegal' )  THEN
                 IF ( data_output_pr_user(1) /= ' ' )  THEN
-                   message_string = 'illegal value for data_output_pr or ' //  &
-                                    'data_output_pr_user = "' //               &
-                                    TRIM( data_output_pr(i) ) // '"'
+                   message_string = 'illegal value for data_output_pr or ' //                      &
+                                    'data_output_pr_user = "' // TRIM( data_output_pr(i) ) // '"'
                    CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
                 ELSE
-                   message_string = 'illegal value for data_output_pr = "' //  &
+                   message_string = 'illegal value for data_output_pr = "' //                      &
                                     TRIM( data_output_pr(i) ) // '"'
                    CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
@@ -2371 +2226 @@
        DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 500 )
           IF ( i > 500 )  THEN
-             message_string = 'number of output quantitities given by data' // &
-                '_output and data_output_user exceeds the limit of 500'
+             message_string = 'number of output quantitities given by data' //                     &
+                              '_output and data_output_user exceeds the limit of 500'
              CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
           ENDIF
@@ -2401 +2256 @@
        var = data_output(i)(1:ilen)
        IF ( ilen > 3 )  THEN
-          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR.                      &
-               data_output(i)(ilen-2:ilen) == '_xz'  .OR.                      &
-               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
+          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR.  data_output(i)(ilen-2:ilen) == '_xz'    &
+               .OR.  data_output(i)(ilen-2:ilen) == '_yz' )  THEN
              k = 1                                             ! 2d data
              var = data_output(i)(1:ilen-3)
@@ -2415 +2269 @@
           CASE ( 'e' )
              IF ( constant_diffusion )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res constant_diffusion = .FALSE.'
+                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
+                                 'constant_diffusion = .FALSE.'
                 CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
              ENDIF
@@ -2422 +2276 @@
 
           CASE ( 'thetal' )
-             IF (  .NOT.  bulk_cloud_model )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                         'res bulk_cloud_model = .TRUE.'
+             IF ( .NOT. bulk_cloud_model )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
+                                 'bulk_cloud_model = .TRUE.'
                 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
              ENDIF
@@ -2430 +2284 @@
 
           CASE ( 'pc', 'pr' )
-             IF (  .NOT.  particle_advection )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requir' // &
-                   'es a "particle_parameters"-NAMELIST in the parameter ' //  &
-                   'file (PARIN)'
+             IF ( .NOT. particle_advection )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
+                                 'a "particle_parameters"-NAMELIST in the parameter file (PARIN)'
                 CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
              ENDIF
@@ -2440 +2293 @@
 
           CASE ( 'q', 'thetav' )
-             IF (  .NOT.  humidity )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res humidity = .TRUE.'
+             IF ( .NOT. humidity )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requires humidity = .TRUE.'
                 CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
              ENDIF
@@ -2450 +2302 @@
           CASE ( 'ql' )
              IF ( .NOT.  ( bulk_cloud_model  .OR.  cloud_droplets ) )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                      'res bulk_cloud_model = .TRUE. or cloud_droplets = .TRUE.'
+                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
+                                 'bulk_cloud_model = .TRUE. or cloud_droplets = .TRUE.'
                 CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
              ENDIF
@@ -2457 +2309 @@
 
           CASE ( 'ql_c', 'ql_v', 'ql_vp' )
-             IF (  .NOT.  cloud_droplets )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res cloud_droplets = .TRUE.'
+             IF ( .NOT. cloud_droplets )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
+                                 'cloud_droplets = .TRUE.'
                 CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
              ENDIF
@@ -2467 +2319 @@
 
           CASE ( 'qv' )
-             IF (  .NOT.  bulk_cloud_model )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res bulk_cloud_model = .TRUE.'
+             IF ( .NOT. bulk_cloud_model )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
+                                 'bulk_cloud_model = .TRUE.'
                 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
              ENDIF
@@ -2475 +2327 @@
 
           CASE ( 's' )
-             IF (  .NOT.  passive_scalar )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res passive_scalar = .TRUE.'
+             IF ( .NOT. passive_scalar )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
+                                 'passive_scalar = .TRUE.'
                 CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
              ENDIF
@@ -2490 +2342 @@
              CONTINUE
 
-          CASE ( 'ghf*', 'lwp*', 'ol*', 'qsurf*', 'qsws*', 'r_a*',             &
-                 'shf*', 'ssurf*', 'ssws*', 't*', 'tsurf*', 'us*',             &
-                 'z0*', 'z0h*', 'z0q*'  )
+          CASE ( 'ghf*', 'lwp*', 'ol*', 'qsurf*', 'qsws*', 'r_a*', 'shf*', 'ssurf*', 'ssws*', 't*',&
+                 'tsurf*', 'us*', 'z0*', 'z0h*', 'z0q*' )
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
-                message_string = 'illegal value for data_output: "' //         &
-                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
-                                 'cross sections are allowed for this value'
+                message_string = 'illegal value for data_output: "' // TRIM( var ) //              &
+                                 '" & only 2d-horizontal cross sections are allowed for this value'
                 CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
              ENDIF
 
              IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. bulk_cloud_model )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res bulk_cloud_model = .TRUE.'
+                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
+                                 'bulk_cloud_model = .TRUE.'
                 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
              ENDIF
              IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT.  humidity )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res humidity = .TRUE.'
+                message_string = 'output of "' // TRIM( var ) // '" requires humidity = .TRUE.'
                 CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
              ENDIF
 
              IF ( TRIM( var ) == 'ghf*'  .AND.  .NOT.  land_surface )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE.'
+                message_string = 'output of "' // TRIM( var ) // '" requires land_surface = .TRUE.'
                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
              ENDIF
 
-             IF ( ( TRIM( var ) == 'r_a*' .OR.  TRIM( var ) == 'ghf*' )        &
-                 .AND.  .NOT.  land_surface  .AND.  .NOT.  urban_surface )     &
-             THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res land_surface = .TRUE. or ' //            &
-                                 'urban_surface = .TRUE.'
+             IF ( (  TRIM( var ) == 'r_a*' .OR. TRIM( var ) == 'ghf*' )  .AND.  .NOT. land_surface &
+                   .AND.  .NOT. urban_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
+                                 'land_surface = .TRUE. or ' // 'urban_surface = .TRUE.'
                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
              ENDIF
 
-             IF ( TRIM( var ) == 'ssws*'  .AND.  .NOT.  passive_scalar )  THEN
-                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
-                                 'res passive_scalar = .TRUE.'
+             IF ( TRIM( var ) == 'ssws*'  .AND.  .NOT. passive_scalar )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
+                                 'passive_scalar = .TRUE.'
                 CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
              ENDIF
@@ -2547 +2393 @@
              IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
 !
-!--          Output of surface latent and sensible heat flux will be in W/m2
-!--          in case of natural- and urban-type surfaces, even if
-!--          flux_output_mode is set to kinematic units.
+!--          Output of surface latent and sensible heat flux will be in W/m2 in case of natural- and
+!--          urban-type surfaces, even if flux_output_mode is set to kinematic units.
              IF ( land_surface  .OR.  urban_surface )  THEN
                 IF ( TRIM( var ) == 'shf*'   )  unit = 'W/m2'
@@ -2562 +2407 @@
              IF ( unit == 'illegal' )  THEN
                 IF ( data_output_user(1) /= ' ' )  THEN
-                   message_string = 'illegal value for data_output or ' //     &
-                         'data_output_user = "' // TRIM( data_output(i) ) // '"'
+                   message_string = 'illegal value for data_output or ' //                         &
+                                    'data_output_user = "' // TRIM( data_output(i) ) // '"'
                    CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
                 ELSE
-                   message_string = 'illegal value for data_output = "' //     &
+                   message_string = 'illegal value for data_output = "' //                         &
                                     TRIM( data_output(i) ) // '"'
                    CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
@@ -2614 +2459 @@
 !-- Averaged 2d or 3d output requires that an averaging interval has been set
     IF ( doav_n > 0  .AND.  averaging_interval == 0.0_wp )  THEN
-       WRITE( message_string, * )  'output of averaged quantity "',            &
-                                   TRIM( doav(1) ), '_av" requires to set a ', &
-                                   'non-zero averaging interval'
+       WRITE( message_string, * )  'output of averaged quantity "', TRIM( doav(1) ),               &
+                                   '_av" requires to set a ', 'non-zero averaging interval'
        CALL message( 'check_parameters', 'PA0323', 1, 2, 0, 6, 0 )
     ENDIF
@@ -2638 +2482 @@
     section(:,3) = section_yz
 
-    IF ( ANY( data_output_xy ) .AND. .NOT. ANY( section(:,1) /= -9999 ) )  THEN
-       WRITE( message_string, * )  'section_xy not defined for requested '  // &
-                                   'xy-cross section output.&At least one ' // &
-                                   'cross section must be given.'
+    IF ( ANY( data_output_xy )  .AND.  .NOT. ANY( section(:,1) /= -9999 ) )  THEN
+       WRITE( message_string, * )  'section_xy not defined for requested xy-cross section ' //     &
+                                   'output.&At least one cross section must be given.'
        CALL message( 'check_parameters', 'PA0681', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( ANY( data_output_xz ) .AND. .NOT. ANY( section(:,2) /= -9999 ) )  THEN
-       WRITE( message_string, * )  'section_xz not defined for requested '  // &
-                                   'xz-cross section output.&At least one ' // &
-                                   'cross section must be given.'
+    IF ( ANY( data_output_xz )  .AND.  .NOT. ANY( section(:,2) /= -9999 ) )  THEN
+       WRITE( message_string, * )  'section_xz not defined for requested xz-cross section ' //     &
+                                   'output.&At least one cross section must be given.'
        CALL message( 'check_parameters', 'PA0681', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( ANY( data_output_yz ) .AND. .NOT. ANY( section(:,3) /= -9999 ) )  THEN
-       WRITE( message_string, * )  'section_yz not defined for requested '  // &
-                                   'yz-cross section output.&At least one ' // &
-                                   'cross section must be given.'
+    IF ( ANY( data_output_yz )  .AND.  .NOT. ANY( section(:,3) /= -9999 ) )  THEN
+       WRITE( message_string, * )  'section_yz not defined for requested yz-cross section ' //     &
+                                   'output.&At least one cross section must be given.'
        CALL message( 'check_parameters', 'PA0681', 1, 2, 0, 6, 0 )
     ENDIF
@@ -2681 +2522 @@
 !-- Check mask conditions
     DO mid = 1, max_masks
-       IF ( data_output_masks(mid,1) /= ' '  .OR.                              &
-            data_output_masks_user(mid,1) /= ' ' ) THEN
+       IF ( data_output_masks(mid,1) /= ' '  .OR.  data_output_masks_user(mid,1) /= ' ' )  THEN
           masks = masks + 1
        ENDIF
@@ -2688 +2528 @@
 
     IF ( masks < 0  .OR.  masks > max_masks )  THEN
-       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
-            '<= ', max_masks, ' (=max_masks)'
+       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ', '<= ', max_masks,     &
+              ' (=max_masks)'
        CALL message( 'check_parameters', 'PA0325', 1, 2, 0, 6, 0 )
     ENDIF
@@ -2697 +2537 @@
        mask_scale(3) = mask_scale_z
        IF ( ANY( mask_scale <= 0.0_wp ) )  THEN
-          WRITE( message_string, * )                                           &
-               'illegal value: mask_scale_x, mask_scale_y and mask_scale_z',   &
-               'must be > 0.0'
+          WRITE( message_string, * )  'illegal value: mask_scale_x, mask_scale_y and mask_scale_z',&
+                 'must be > 0.0'
           CALL message( 'check_parameters', 'PA0326', 1, 2, 0, 6, 0 )
        ENDIF
 !
-!--    Generate masks for masked data output
-!--    Parallel netcdf output is not tested so far for masked data, hence
-!--    netcdf_data_format is switched back to non-parallel output.
+!--    Generate masks for masked data output.
+!--    Parallel netcdf output is not tested so far for masked data, hence netcdf_data_format is
+!--    switched back to non-parallel output.
        netcdf_data_format_save = netcdf_data_format
        IF ( netcdf_data_format > 4 )  THEN
           IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
           IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
-          message_string = 'netCDF file formats '//                            &
-                           '5 (parallel netCDF 4) and ' //                     &
-                           '6 (parallel netCDF 4 Classic model) '//            &
-                           '& are currently not supported (not yet tested) ' //&
-                           'for masked data. &Using respective non-parallel' //&
+          message_string = 'netCDF file formats '// '5 (parallel netCDF 4) and 6 (parallel ' //    &
+                           'netCDF 4 Classic model) & are currently not supported (not yet ' //    &
+                           'tested) for masked data. &Using respective non-parallel' //            &
                            ' output for masked data.'
           CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
@@ -2728 +2565 @@
        CONTINUE
 #else
-       message_string = 'netCDF: netCDF4 format requested but no ' //          &
-                        'cpp-directive __netcdf4 given & switch '  //          &
-                        'back to 64-bit offset format'
+       message_string = 'netCDF: netCDF4 format requested but no ' //                              &
+                        'cpp-directive __netcdf4 given & switch back to 64-bit offset format'
        CALL message( 'check_parameters', 'PA0171', 0, 1, 0, 6, 0 )
        netcdf_data_format = 2
@@ -2739 +2575 @@
        CONTINUE
 #else
-       message_string = 'netCDF: netCDF4 parallel output requested but no ' // &
-                        'cpp-directive __netcdf4_parallel given & switch '   //&
+       message_string = 'netCDF: netCDF4 parallel output requested but no ' //                     &
+                        'cpp-directive __netcdf4_parallel given & switch ' //                      &
                         'back to netCDF4 non-parallel output'
        CALL message( 'check_parameters', 'PA0099', 0, 1, 0, 6, 0 )
@@ -2749 +2585 @@
 !
 !-- Calculate fixed number of output time levels for parallel netcdf output.
-!-- The time dimension has to be defined as limited for parallel output,
-!-- because otherwise the I/O performance drops significantly.
+!-- The time dimension has to be defined as limited for parallel output, because otherwise the I/O
+!-- performance drops significantly.
     IF ( netcdf_data_format > 4 )  THEN
 
 !
-!--    Check if any of the follwoing data output interval is 0.0s, which is
-!--    not allowed for parallel output.
+!--    Check if any of the follwoing data output interval is 0.0s, which is not allowed for parallel
+!--    output.
        CALL check_dt_do( dt_do3d,           'dt_do3d'           )
        CALL check_dt_do( dt_do2d_xy,        'dt_do2d_xy'        )
@@ -2762 +2598 @@
        CALL check_dt_do( dt_data_output_av, 'dt_data_output_av' )
 
-!--    Set needed time levels (ntdim) to
-!--    saved time levels + to be saved time levels.
-       ntdim_3d(0) = do3d_time_count(0) + CEILING(                                    &
-                     ( end_time - MAX(                                                &
-                         MERGE(skip_time_do3d, skip_time_do3d + spinup_time,          &
-                               data_output_during_spinup ),                           &
-                         simulated_time_at_begin )                                    &
+!--    Set needed time levels (ntdim) to saved time levels + to be saved time levels.
+       ntdim_3d(0) = do3d_time_count(0) + CEILING(                                                 &
+                     ( end_time - MAX(                                                             &
+                         MERGE( skip_time_do3d, skip_time_do3d + spinup_time,                      &
+                                data_output_during_spinup ),                                       &
+                         simulated_time_at_begin )                                                 &
                      ) / dt_do3d )
        IF ( do3d_at_begin ) ntdim_3d(0) = ntdim_3d(0) + 1
 
-       ntdim_3d(1) = do3d_time_count(1) + CEILING(                                    &
-                     ( end_time - MAX(                                                &
-                         MERGE(   skip_time_data_output_av, skip_time_data_output_av  &
-                                + spinup_time, data_output_during_spinup ),           &
-                         simulated_time_at_begin )                                    &
+       ntdim_3d(1) = do3d_time_count(1) + CEILING(                                                 &
+                     ( end_time - MAX(                                                             &
+                         MERGE( skip_time_data_output_av, skip_time_data_output_av + spinup_time,  &
+                                data_output_during_spinup ),                                       &
+                         simulated_time_at_begin )                                                 &
                      ) / dt_data_output_av )
 
-       ntdim_2d_xy(0) = do2d_xy_time_count(0) + CEILING(                              &
-                        ( end_time - MAX(                                             &
-                           MERGE(skip_time_do2d_xy, skip_time_do2d_xy + spinup_time,  &
-                                 data_output_during_spinup ),                         &
-                           simulated_time_at_begin )                                  &
+       ntdim_2d_xy(0) = do2d_xy_time_count(0) + CEILING(                                           &
+                        ( end_time - MAX(                                                          &
+                            MERGE( skip_time_do2d_xy, skip_time_do2d_xy + spinup_time,             &
+                                   data_output_during_spinup ),                                    &
+                            simulated_time_at_begin )                                              &
                         ) / dt_do2d_xy )
 
-       ntdim_2d_xz(0) = do2d_xz_time_count(0) + CEILING(                              &
-                        ( end_time - MAX(                                             &
-                         MERGE(skip_time_do2d_xz, skip_time_do2d_xz + spinup_time,    &
-                               data_output_during_spinup ),                           &
-                         simulated_time_at_begin )                                    &
+       ntdim_2d_xz(0) = do2d_xz_time_count(0) + CEILING(                                           &
+                        ( end_time - MAX(                                                          &
+                            MERGE( skip_time_do2d_xz, skip_time_do2d_xz + spinup_time,             &
+                                   data_output_during_spinup ),                                    &
+                            simulated_time_at_begin )                                              &
                         ) / dt_do2d_xz )
 
-       ntdim_2d_yz(0) = do2d_yz_time_count(0) + CEILING(                              &
-                        ( end_time - MAX(                                             &
-                         MERGE(skip_time_do2d_yz, skip_time_do2d_yz + spinup_time,    &
-                               data_output_during_spinup ),                           &
-                         simulated_time_at_begin )                                    &
+       ntdim_2d_yz(0) = do2d_yz_time_count(0) + CEILING(                                           &
+                        ( end_time - MAX(                                                          &
+                            MERGE( skip_time_do2d_yz, skip_time_do2d_yz + spinup_time,             &
+                                   data_output_during_spinup ),                                    &
+                            simulated_time_at_begin )                                              &
                         ) / dt_do2d_yz )
 
@@ -2806 +2641 @@
        ENDIF
 !
-!--    Please note, for averaged 2D data skip_time_data_output_av is the relavant
-!--    output control parameter.
-       ntdim_2d_xy(1) = do2d_xy_time_count(1) + CEILING(                              &
-                     ( end_time - MAX( MERGE( skip_time_data_output_av,               &
-                                              skip_time_data_output_av + spinup_time, &
-                                              data_output_during_spinup ),            &
-                                       simulated_time_at_begin )                      &
-                     ) / dt_data_output_av )
-
-       ntdim_2d_xz(1) = do2d_xz_time_count(1) + CEILING(                              &
-                     ( end_time - MAX( MERGE( skip_time_data_output_av,               &
-                                              skip_time_data_output_av + spinup_time, &
-                                              data_output_during_spinup ),            &
-                                       simulated_time_at_begin )                      &
-                     ) / dt_data_output_av )
-
-       ntdim_2d_yz(1) = do2d_yz_time_count(1) + CEILING(                              &
-                     ( end_time - MAX( MERGE( skip_time_data_output_av,               &
-                                              skip_time_data_output_av + spinup_time, &
-                                              data_output_during_spinup ),            &
-                                       simulated_time_at_begin )                      &
-                     ) / dt_data_output_av )
+!--    Please note, for averaged 2D data skip_time_data_output_av is the relavant output control
+!--    parameter.
+       ntdim_2d_xy(1) = do2d_xy_time_count(1) + CEILING(                                           &
+                        ( end_time - MAX( MERGE( skip_time_data_output_av,                         &
+                                                 skip_time_data_output_av + spinup_time,           &
+                                                 data_output_during_spinup ),                      &
+                                          simulated_time_at_begin )                                &
+                        ) / dt_data_output_av )
+
+       ntdim_2d_xz(1) = do2d_xz_time_count(1) + CEILING(                                           &
+                        ( end_time - MAX( MERGE( skip_time_data_output_av,                         &
+                                                 skip_time_data_output_av + spinup_time,           &
+                                                 data_output_during_spinup ),                      &
+                                          simulated_time_at_begin )                                &
+                        ) / dt_data_output_av )
+
+       ntdim_2d_yz(1) = do2d_yz_time_count(1) + CEILING(                                           &
+                        ( end_time - MAX( MERGE( skip_time_data_output_av,                         &
+                                                 skip_time_data_output_av + spinup_time,           &
+                                                 data_output_during_spinup ),                      &
+                                          simulated_time_at_begin )                                &
+                        ) / dt_data_output_av )
 
     ENDIF
@@ -2839 +2674 @@
        ELSE
           IF ( prandtl_number < 0.0_wp )  THEN
-             WRITE( message_string, * )  'prandtl_number = ', prandtl_number,  &
-                                         ' < 0.0'
+             WRITE( message_string, * )  'prandtl_number = ', prandtl_number, ' < 0.0'
              CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
           ENDIF
@@ -2846 +2680 @@
 
           IF ( constant_flux_layer )  THEN
-             message_string = 'constant_flux_layer is not allowed with fixed ' &
-                              // 'value of km'
+             message_string = 'constant_flux_layer is not allowed with fixed value of km'
              CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
           ENDIF
@@ -2854 +2687 @@
 
 !
-!-- In case of non-cyclic lateral boundaries and a damping layer for the
-!-- potential temperature, check the width of the damping layer
+!-- In case of non-cyclic lateral boundaries and a damping layer for the potential temperature,
+!-- check the width of the damping layer
     IF ( bc_lr /= 'cyclic' ) THEN
-       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
-            pt_damping_width > REAL( (nx+1) * dx ) )  THEN
+       IF ( pt_damping_width < 0.0_wp  .OR.  pt_damping_width > REAL( (nx+1) * dx ) )  THEN
           message_string = 'pt_damping_width out of range'
           CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
@@ -2865 +2697 @@
 
     IF ( bc_ns /= 'cyclic' )  THEN
-       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
-            pt_damping_width > REAL( (ny+1) * dy ) )  THEN
+       IF ( pt_damping_width < 0.0_wp  .OR.  pt_damping_width > REAL( (ny+1) * dy ) )  THEN
           message_string = 'pt_damping_width out of range'
           CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
@@ -2875 +2706 @@
 !-- Check value range for zeta = z/L
     IF ( zeta_min >= zeta_max )  THEN
-       WRITE( message_string, * )  'zeta_min = ', zeta_min, ' must be less ',  &
-                                   'than zeta_max = ', zeta_max
+       WRITE( message_string, * )  'zeta_min = ', zeta_min, ' must be less ', 'than zeta_max = ',  &
+              zeta_max
        CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
     ENDIF
@@ -2882 +2713 @@
 !
 !-- Check random generator
-    IF ( (random_generator /= 'system-specific'      .AND.                     &
-          random_generator /= 'random-parallel'   )  .AND.                     &
+    IF ( (random_generator /= 'system-specific'      .AND.                                         &
+          random_generator /= 'random-parallel'   )  .AND.                                         &
           random_generator /= 'numerical-recipes' )  THEN
-       message_string = 'unknown random generator: random_generator = "' //    &
+       message_string = 'unknown random generator: random_generator = "' //                        &
                         TRIM( random_generator ) // '"'
        CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
@@ -2901 +2732 @@
        ENDIF
     ELSEIF ( disturbance_level_b < zu(3) )  THEN
-       WRITE( message_string, * )  'disturbance_level_b = ',                   &
-                           disturbance_level_b, ' must be >= ', zu(3), '(zu(3))'
+       WRITE( message_string, * )  'disturbance_level_b = ', disturbance_level_b, ' must be >= ',  &
+              zu(3), '(zu(3))'
        CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
     ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
-       WRITE( message_string, * )  'disturbance_level_b = ',                   &
-                   disturbance_level_b, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
+       WRITE( message_string, * )  'disturbance_level_b = ', disturbance_level_b, ' must be <= ',  &
+              zu(nzt-2), '(zu(nzt-2))'
        CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
     ELSE
@@ -2926 +2757 @@
        ENDIF
     ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
-       WRITE( message_string, * )  'disturbance_level_t = ',                   &
-                   disturbance_level_t, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
+       WRITE( message_string, * )  'disturbance_level_t = ', disturbance_level_t, ' must be <= ',  &
+              zu(nzt-2), '(zu(nzt-2))'
        CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
     ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
-       WRITE( message_string, * )  'disturbance_level_t = ',                   &
-                   disturbance_level_t, ' must be >= disturbance_level_b = ',  &
-                   disturbance_level_b
+       WRITE( message_string, * )  'disturbance_level_t = ', disturbance_level_t,                  &
+             ' must be >= disturbance_level_b = ', disturbance_level_b
        CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
     ELSE
@@ -2945 +2775 @@
 !
 !-- Check again whether the levels determined this way are ok.
-!-- Error may occur at automatic determination and too few grid points in
-!-- z-direction.
+!-- Error may occur at automatic determination and too few grid points in z-direction.
     IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
-       WRITE( message_string, * )  'disturbance_level_ind_t = ',               &
-                disturbance_level_ind_t, ' must be >= ',                       &
-                'disturbance_level_ind_b = ', disturbance_level_ind_b
+       WRITE( message_string, * )  'disturbance_level_ind_t = ', disturbance_level_ind_t,          &
+              ' must be >= ', 'disturbance_level_ind_b = ', disturbance_level_ind_b
        CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
     ENDIF
@@ -2956 +2784 @@
 !
 !-- Determine the horizontal index range for random perturbations.
-!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
-!-- near the inflow and the perturbation area is further limited to ...(1)
-!-- after the initial phase of the flow.
+!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed near the inflow and
+!-- the perturbation area is further limited to ...(1) after the initial phase of the flow.
 
     IF ( bc_lr /= 'cyclic' )  THEN
@@ -2964 +2791 @@
           inflow_disturbance_begin = MIN( 10, nx/2 )
        ENDIF
-       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
-       THEN
+       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )  THEN
           message_string = 'inflow_disturbance_begin out of range'
           CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
@@ -2972 +2798 @@
           inflow_disturbance_end = MIN( 100, 3*nx/4 )
        ENDIF
-       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
-       THEN
+       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )  THEN
           message_string = 'inflow_disturbance_end out of range'
           CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
@@ -2981 +2806 @@
           inflow_disturbance_begin = MIN( 10, ny/2 )
        ENDIF
-       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
-       THEN
+       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )  THEN
           message_string = 'inflow_disturbance_begin out of range'
           CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
@@ -2989 +2813 @@
           inflow_disturbance_end = MIN( 100, 3*ny/4 )
        ENDIF
-       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
-       THEN
+       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )  THEN
           message_string = 'inflow_disturbance_end out of range'
           CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
@@ -3045 +2868 @@
 
 !
-!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
-!-- boundary (so far, a turbulent inflow is realized from the left side only)
+!-- A turbulent inflow requires Dirichlet conditions at the respective inflow boundary (so far, a
+!-- turbulent inflow is realized from the left side only).
     IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
-       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' //      &
+       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' //                          &
                         'condition at the inflow boundary'
        CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
@@ -3054 +2877 @@
 
 !
-!-- Turbulent inflow requires that 3d arrays have been cyclically filled with
-!-- data from prerun in the first main run
-    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND. &
+!-- Turbulent inflow requires that 3d arrays have been cyclically filled with data from prerun in
+!-- the first main run
+    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND.                     &
          initializing_actions /= 'read_restart_data' )  THEN
-       message_string = 'turbulent_inflow = .T. requires ' //                  &
-                        'initializing_actions = ''cyclic_fill'' or ' //        &
+       message_string = 'turbulent_inflow = .T. requires ' //                                      &
+                        'initializing_actions = ''cyclic_fill'' or ' //                            &
                         'initializing_actions = ''read_restart_data'' '
        CALL message( 'check_parameters', 'PA0055', 1, 2, 0, 6, 0 )
@@ -3071 +2894 @@
 !--    Calculate the index of the recycling plane
        IF ( recycling_width <= dx  .OR.  recycling_width >= nx * dx )  THEN
-          WRITE( message_string, * )  'illegal value for recycling_width: ',   &
-                                      recycling_width
+          WRITE( message_string, * )  'illegal value for recycling_width: ', recycling_width
           CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
        ENDIF
@@ -3080 +2902 @@
 !
 !--   Check for correct input of recycling method for thermodynamic quantities
-       IF ( TRIM( recycling_method_for_thermodynamic_quantities ) /= 'turbulent_fluctuation' .AND. &
+       IF ( TRIM( recycling_method_for_thermodynamic_quantities ) /= 'turbulent_fluctuation'  .AND.&
             TRIM( recycling_method_for_thermodynamic_quantities ) /= 'absolute_value' )  THEN
           WRITE( message_string, * )  'unknown recycling method for thermodynamic quantities: ',   &
@@ -3092 +2914 @@
     IF ( turbulent_outflow )  THEN
 !
-!--    Turbulent outflow requires Dirichlet conditions at the respective inflow
-!--    boundary (so far, a turbulent outflow is realized at the right side only)
+!--    Turbulent outflow requires Dirichlet conditions at the respective inflow boundary (so far, a
+!--    turbulent outflow is realized at the right side only).
        IF ( bc_lr /= 'dirichlet/radiation' )  THEN
-          message_string = 'turbulent_outflow = .T. requires ' //              &
-                           'bc_lr = "dirichlet/radiation"'
+          message_string = 'turbulent_outflow = .T. requires bc_lr = "dirichlet/radiation"'
           CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
        ENDIF
 !
 !--    The ouflow-source plane must lay inside the model domain
-       IF ( outflow_source_plane < dx  .OR.  &
-            outflow_source_plane > nx * dx )  THEN
-          WRITE( message_string, * )  'illegal value for outflow_source'//     &
-                                      '_plane: ', outflow_source_plane
+       IF ( outflow_source_plane < dx  .OR.  outflow_source_plane > nx * dx )  THEN
+          WRITE( message_string, * )  'illegal value for outflow_source'// '_plane: ',             &
+                 outflow_source_plane
           CALL message( 'check_parameters', 'PA0145', 1, 2, 0, 6, 0 )
        ENDIF
@@ -3116 +2936 @@
           damp_level_ind_1d = nzt + 1
        ELSEIF ( damp_level_1d < 0.0_wp  .OR.  damp_level_1d > zu(nzt+1) )  THEN
-          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d,       &
-                 ' must be >= 0.0 and <= ', zu(nzt+1), '(zu(nzt+1))'
+          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d, ' must be >= 0.0 and <= ',&
+                 zu(nzt+1), '(zu(nzt+1))'
           CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
        ELSE
@@ -3131 +2951 @@
 !
 !-- Check some other 1d-model parameters
-    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
+    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                                       &
          TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
-       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) //  &
-                        '" is unknown'
+       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) // '" is unknown'
        CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND.                     &
-         TRIM( dissipation_1d ) /= 'detering'  .AND.                           &
+    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND.                                         &
+         TRIM( dissipation_1d ) /= 'detering'        .AND.                                         &
          TRIM( dissipation_1d ) /= 'prognostic' )  THEN
-       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
-                        '" is unknown'
+       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) // '" is unknown'
        CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
+    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                                       &
          TRIM( dissipation_1d ) == 'as_in_3d_model' )  THEN
-       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
+       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //                          &
                         '" requires mixing_length_1d = "as_in_3d_model"'
        CALL message( 'check_parameters', 'PA0485', 1, 2, 0, 6, 0 )
@@ -3152 +2970 @@
 
 !
-!-- Set time for the next user defined restart (time_restart is the
-!-- internal parameter for steering restart events)
+!-- Set time for the next user defined restart (time_restart is the internal parameter for steering
+!-- restart events)
     IF ( restart_time /= 9999999.9_wp )  THEN
        IF ( restart_time > time_since_reference_point )  THEN
@@ -3160 +2978 @@
     ELSE
 !
-!--    In case of a restart run, set internal parameter to default (no restart)
-!--    if the NAMELIST-parameter restart_time is omitted
+!--    In case of a restart run, set internal parameter to default (no restart) if the
+!--    NAMELIST-parameter restart_time is omitted
        time_restart = 9999999.9_wp
     ENDIF
@@ -3168 +2986 @@
 !-- Check pressure gradient conditions
     IF ( dp_external  .AND.  conserve_volume_flow )  THEN
-       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
-            'w are .TRUE. but one of them must be .FALSE.'
+       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo',                     &
+              'w are .TRUE. but one of them must be .FALSE.'
        CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
     ENDIF
     IF ( dp_external )  THEN
        IF ( dp_level_b < zu(nzb)  .OR.  dp_level_b > zu(nzt) )  THEN
-          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ', &
-               ' of range [zu(nzb), zu(nzt)]'
+          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ',                     &
+                 ' of range [zu(nzb), zu(nzt)]'
           CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
        ENDIF
        IF ( .NOT. ANY( dpdxy /= 0.0_wp ) )  THEN
-          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
-               'ro, i.e. the external pressure gradient will not be applied'
+          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze',                     &
+                 'ro, i.e. the external pressure gradient will not be applied'
           CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
        ENDIF
     ENDIF
     IF ( ANY( dpdxy /= 0.0_wp )  .AND.  .NOT.  dp_external )  THEN
-       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',     &
-            '.FALSE., i.e. the external pressure gradient & will not be applied'
+       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',                         &
+              '.FALSE., i.e. the external pressure gradient & will not be applied'
        CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
     ENDIF
@@ -3194 +3012 @@
           conserve_volume_flow_mode = 'initial_profiles'
 
-       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
-            TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
-          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ',   &
-               conserve_volume_flow_mode
+       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles'  .AND.                     &
+                TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
+          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ',                       &
+                 conserve_volume_flow_mode
           CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
        ENDIF
-       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND.                &
-          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
-          WRITE( message_string, * )  'non-cyclic boundary conditions ',       &
-               'require  conserve_volume_flow_mode = ''initial_profiles'''
+       IF ( ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND.                                   &
+            TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
+          WRITE( message_string, * )  'non-cyclic boundary conditions ',                           &
+                 'require  conserve_volume_flow_mode = ''initial_profiles'''
           CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
        ENDIF
     ENDIF
-    IF ( ( u_bulk /= 0.0_wp  .OR.  v_bulk /= 0.0_wp )  .AND.                   &
-         ( .NOT. conserve_volume_flow  .OR.                                    &
-         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
-       WRITE( message_string, * )  'nonzero bulk velocity requires ',          &
-            'conserve_volume_flow = .T. and ',                                 &
-            'conserve_volume_flow_mode = ''bulk_velocity'''
+    IF ( ( u_bulk /= 0.0_wp .OR. v_bulk /= 0.0_wp )  .AND.                                         &
+         ( .NOT. conserve_volume_flow .OR. TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )&
+    THEN
+       WRITE( message_string, * )  'nonzero bulk velocity requires ',                              &
+              'conserve_volume_flow = .T. and ', 'conserve_volume_flow_mode = ''bulk_velocity'''
        CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
     ENDIF
 
 !
-!-- Prevent empty time records in volume, cross-section and masked data in case
-!-- of non-parallel netcdf-output in restart runs
+!-- Prevent empty time records in volume, cross-section and masked data in case of non-parallel
+!-- netcdf-output in restart runs
     IF ( netcdf_data_format < 5 )  THEN
        IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
@@ -3232 +3049 @@
 !
 !-- Check roughness length, which has to be smaller than dz/2
-    IF ( ( constant_flux_layer .OR.  &
-           INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )       &
-         .AND. roughness_length >= 0.5 * dz(1) )  THEN
+    IF ( ( constant_flux_layer .OR. INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  &
+         .AND.  roughness_length >= 0.5 * dz(1) )  THEN
        message_string = 'roughness_length must be smaller than dz/2'
        CALL message( 'check_parameters', 'PA0424', 1, 2, 0, 6, 0 )
@@ -3248 +3064 @@
 !-- Check if topography is read from file in case of complex terrain simulations
     IF ( complex_terrain  .AND.  TRIM( topography ) /= 'read_from_file' )  THEN
-       message_string = 'complex_terrain requires topography' //               &
-                        ' = ''read_from_file'''
+       message_string = 'complex_terrain requires topography = ''read_from_file'''
        CALL message( 'check_parameters', 'PA0295', 1, 2, 0, 6, 0 )
     ENDIF
 
 !
-!-- Check if vertical grid stretching is switched off in case of complex
-!-- terrain simulations
-    IF ( complex_terrain  .AND.                                                &
-         ANY( dz_stretch_level_start /= -9999999.9_wp ) )  THEN
-       message_string = 'Vertical grid stretching is not allowed for ' //      &
-                        'complex_terrain = .T.'
+!-- Check if vertical grid stretching is switched off in case of complex terrain simulations
+    IF ( complex_terrain  .AND.  ANY( dz_stretch_level_start /= -9999999.9_wp ) )  THEN
+       message_string = 'Vertical grid stretching is not allowed for complex_terrain = .TRUE.'
        CALL message( 'check_parameters', 'PA0473', 1, 2, 0, 6, 0 )
     ENDIF
@@ -3267 +3079 @@
  CONTAINS
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Check the length of data output intervals. In case of parallel NetCDF output
-!> the time levels of the output files need to be fixed. Therefore setting the
-!> output interval to 0.0s (usually used to output each timestep) is not
-!> possible as long as a non-fixed timestep is used.
-!------------------------------------------------------------------------------!
+!> Check the length of data output intervals. In case of parallel NetCDF output the time levels of
+!> the output files need to be fixed. Therefore setting the output interval to 0.0s (usually used to
+!> output each timestep) is not possible as long as a non-fixed timestep is used.
+!--------------------------------------------------------------------------------------------------!
 
     SUBROUTINE check_dt_do( dt_do, dt_do_name )
@@ -3286 +3097 @@
        IF ( dt_do == 0.0_wp )  THEN
           IF ( dt_fixed )  THEN
-             WRITE( message_string, '(A,F9.4,A)' )  'Output at every '  //     &
-                    'timestep is wanted (' // dt_do_name // ' = 0.0).&'//      &
-                    'The output interval is set to the fixed timestep dt '//   &
-                    '= ', dt, 's.'
+             WRITE( message_string, '(A,F9.4,A)' )  'Output at every timestep is wanted (' //      &
+                    dt_do_name // ' = 0.0).&'//                                                    &
+                    'The output interval is set to the fixed timestep dt '// '= ', dt, 's.'
              CALL message( 'check_parameters', 'PA0060', 0, 0, 0, 6, 0 )
              dt_do = dt
           ELSE
-             message_string = dt_do_name // ' = 0.0 while using a ' //         &
-                              'variable timestep and parallel netCDF4 ' //     &
-                              'is not allowed.'
+             message_string = dt_do_name // ' = 0.0 while using a ' //                             &
+                              'variable timestep and parallel netCDF4 is not allowed.'
              CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
           ENDIF
@@ -3304 +3113 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Set the bottom and top boundary conditions for humidity and scalars.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 
     SUBROUTINE set_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, err_nr_b, err_nr_t )
@@ -3315 +3124 @@
        IMPLICIT NONE
 
-       CHARACTER (LEN=1)   ::  sq         !< name of scalar quantity
        CHARACTER (LEN=*)   ::  bc_b       !< bottom boundary condition
        CHARACTER (LEN=*)   ::  bc_t       !< top boundary condition
        CHARACTER (LEN=*)   ::  err_nr_b   !< error number if bottom bc is unknown
        CHARACTER (LEN=*)   ::  err_nr_t   !< error number if top bc is unknown
+       CHARACTER (LEN=1)   ::  sq         !< name of scalar quantity
+
 
        INTEGER(iwp)        ::  ibc_b      !< index for bottom boundary condition
@@ -3325 +3135 @@
 
 !
-!--    Set Integer flags and check for possilbe errorneous settings for bottom
-!--    boundary condition
+!--    Set Integer flags and check for possilbe errorneous settings for bottom boundary condition
        IF ( bc_b == 'dirichlet' )  THEN
           ibc_b = 0
@@ -3332 +3141 @@
           ibc_b = 1
        ELSE
-          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
-                           '_b ="' // TRIM( bc_b ) // '"'
+          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // '_b ="' //           &
+                           TRIM( bc_b ) // '"'
           CALL message( 'check_parameters', err_nr_b, 1, 2, 0, 6, 0 )
        ENDIF
 !
-!--    Set Integer flags and check for possilbe errorneous settings for top
-!--    boundary condition
+!--    Set Integer flags and check for possilbe errorneous settings for top boundary condition
        IF ( bc_t == 'dirichlet' )  THEN
           ibc_t = 0
@@ -3348 +3156 @@
           ibc_t = 3
        ELSE
-          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
-                           '_t ="' // TRIM( bc_t ) // '"'
+          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // '_t ="' //           &
+                           TRIM( bc_t ) // '"'
           CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 )
        ENDIF
@@ -3358 +3166 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Check for consistent settings of bottom boundary conditions for humidity
-!> and scalars.
-!------------------------------------------------------------------------------!
-
-    SUBROUTINE check_bc_scalars( sq, bc_b, ibc_b,                      &
-                                 err_nr_1, err_nr_2,                   &
-                                 constant_flux, surface_initial_change )
+!> Check for consistent settings of bottom boundary conditions for humidity and scalars.
+!--------------------------------------------------------------------------------------------------!
+
+    SUBROUTINE check_bc_scalars( sq, bc_b, ibc_b, err_nr_1, err_nr_2, constant_flux,               &
+                                 surface_initial_change )
 
 
        IMPLICIT NONE
 
-       CHARACTER (LEN=1)   ::  sq                       !< name of scalar quantity
        CHARACTER (LEN=*)   ::  bc_b                     !< bottom boundary condition
        CHARACTER (LEN=*)   ::  err_nr_1                 !< error number of first error
        CHARACTER (LEN=*)   ::  err_nr_2                 !< error number of second error
+       CHARACTER (LEN=1)   ::  sq                       !< name of scalar quantity
+
 
        INTEGER(iwp)        ::  ibc_b                    !< index of bottom boundary condition
@@ -3384 +3191 @@
 
 !
-!--    A given surface value implies Dirichlet boundary condition for
-!--    the respective quantity. In this case specification of a constant flux is
-!--    forbidden. However, an exception is made for large-scale forcing as well
-!--    as land-surface model.
+!--    A given surface value implies Dirichlet boundary condition for the respective quantity. In
+!--    this case specification of a constant flux is forbidden. However, an exception is made for
+!--    large-scale forcing as well as land-surface model.
        IF ( .NOT. land_surface  .AND.  .NOT. large_scale_forcing )  THEN
           IF ( ibc_b == 0  .AND.  constant_flux )  THEN
-             message_string = 'boundary condition: bc_' // TRIM( sq ) //       &
-                              '_b ' // '= "' // TRIM( bc_b ) //                &
-                              '" is not allowed with prescribed surface flux'
+             message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b = "' //               &
+                              TRIM( bc_b ) // '" is not allowed with prescribed surface flux'
              CALL message( 'check_parameters', err_nr_1, 1, 2, 0, 6, 0 )
           ENDIF
        ENDIF
        IF ( constant_flux  .AND.  surface_initial_change /= 0.0_wp )  THEN
-          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
-                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
-                 surface_initial_change
+          WRITE( message_string, * )  'a prescribed surface flux is not allo', 'wed with ', sq,    &
+                '_surface_initial_change (/=0) = ', surface_initial_change
           CALL message( 'check_parameters', err_nr_2, 1, 2, 0, 6, 0 )
        ENDIF

palm/trunk/SOURCE/chem_emissions_mod.f90

@@ -1 +1 @@
 !> @file chem_emissions_mod.f90
-!--------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! This file is part of 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 2018-2020 Leibniz Universitaet Hannover
 ! Copyright 2018-2020 Freie Universitaet Berlin
 ! Copyright 2018-2020 Karlsruhe Institute of Technology
-!--------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-!
+! 
 ! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+! file re-formatted to follow the PALM coding standard
+!
+! 4481 2020-03-31 18:55:54Z maronga
 ! Implemented on-demand read mode for LOD 2 (ECC)
 !  - added following module global variables
@@ -55 +57 @@
 !
 ! 4242 2019-09-27 12:59:10Z suehring
-! Adjust index_hh access to new definition accompanied with new 
+! Adjust index_hh access to new definition accompanied with new
 ! palm_date_time_mod. Note, this is just a preleminary fix. (E Chan)
-! 
+!
 ! 4230 2019-09-11 13:58:14Z suehring
-! Bugfix, consider that time_since_reference_point can be also negative when 
+! Bugfix, consider that time_since_reference_point can be also negative when
 ! time indices are determined.
-! 
+!
 ! 4227 2019-09-10 18:04:34Z gronemeier
 ! implement new palm_date_time_mod
-! 
+!
 ! 4223 2019-09-10 09:20:47Z gronemeier
 ! Unused routine chem_emissions_check_parameters commented out due to uninitialized content
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 4157 2019-08-14 09:19:12Z suehring
 ! Replace global arrays also in mode_emis branch
-! 
+!
 ! 4154 2019-08-13 13:35:59Z suehring
 ! Replace global arrays for emissions by local ones.
-! 
+!
 ! 4144 2019-08-06 09:11:47Z raasch
 ! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
-! 
+!
 ! 4055 2019-06-27 09:47:29Z suehring
-! - replaced local thermo. constants w/ module definitions in
-!   basic_constants_and_equations_mod (rgas_univ, p_0, r_d_cp) 
+! - replaced local thermo. constants w/ module definitions in basic_constants_and_equations_mod
+!   (rgas_univ, p_0, r_d_cp)
 ! - initialize array emis_distribution immediately following allocation
-! - lots of minor formatting changes based on review sesson in 20190325 
-!   (E.C. Chan)
-! 
+! - lots of minor formatting changes based on review sesson in 20190325 (E.C. Chan)
+!
 ! 3968 2019-05-13 11:04:01Z suehring
-! - in subroutine chem_emissions_match replace all decision structures relating to
-!   mode_emis to emiss_lod
+! - in subroutine chem_emissions_match replace all decision structures relating to mode_emis to
+!   emiss_lod
 ! - in subroutine chem_check_parameters replace emt%nspec with emt%n_emiss_species
 ! - spring cleaning (E.C. Chan)
-! 
+!
 ! 3885 2019-04-11 11:29:34Z kanani
-! Changes related to global restructuring of location messages and introduction 
-! of additional debug messages
-! 
+! Changes related to global restructuring of location messages and introduction of additional debug
+! messages
+!
 ! 3831 2019-03-28 09:11:22Z forkel
-! added nvar to USE chem_gasphase_mod (chem_modules will not include nvar anymore) 
-! 
+! added nvar to USE chem_gasphase_mod (chem_modules will not include nvar anymore)
+!
 ! 3788 2019-03-07 11:40:09Z banzhafs
 ! Removed unused variables from chem_emissions_mod
-! 
-! 3772 2019-02-28 15:51:57Z suehring 
-! - In case of parametrized emissions, assure that emissions are only on natural
-!   surfaces (i.e. streets) and not on urban surfaces. 
+!
+! 3772 2019-02-28 15:51:57Z suehring
+! - In case of parametrized emissions, assure that emissions are only on natural surfaces
+!   (i.e. streets) and not on urban surfaces.
 ! - some unnecessary if clauses removed
 !
@@ -117 +118 @@
 ! @author Sabine Banzhaf  (FU-Berlin)
 ! @author Martijn Schaap  (FU-Berlin, TNO Utrecht)
-! 
+!
 ! Description:
 ! ------------
@@ -131 +132 @@
 !> @note <Enter notes on the module>
 !> @bug  <Enter known bugs here>
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 
  MODULE chem_emissions_mod
 
-    USE arrays_3d,                                                          &
+    USE arrays_3d,                                                                                 &
         ONLY:  rho_air
 
-    USE basic_constants_and_equations_mod,                                  &
-        ONLY:  rgas_univ, p_0, rd_d_cp
-
-    USE control_parameters,                                                 &
-        ONLY:  debug_output,                                                &
-               end_time, message_string, initializing_actions,              &
-               intermediate_timestep_count, dt_3d
- 
+    USE basic_constants_and_equations_mod,                                                         &
+        ONLY:  p_0, rd_d_cp, rgas_univ
+
+    USE control_parameters,                                                                        &
+        ONLY:  debug_output, end_time, initializing_actions, intermediate_timestep_count,          &
+               message_string, dt_3d
+
     USE indices
 
@@ -154 +154 @@
 #endif
 
-    USE netcdf_data_input_mod,                                               &
+    USE netcdf_data_input_mod,                                                                     &
         ONLY: chem_emis_att_type, chem_emis_val_type
 
-    USE chem_gasphase_mod,                                                   &
+    USE chem_gasphase_mod,                                                                         &
         ONLY: nvar, spc_names
- 
+
     USE chem_modules
 
-    USE statistics,                                                          &   
+    USE statistics,                                                                                &
         ONLY:  weight_pres
 
@@ -169 +169 @@
 !-- Added new palm_date_time_mod for on-demand emission reading
 
-    USE palm_date_time_mod,                                                  &
+    USE palm_date_time_mod,                                                                        &
         ONLY:  get_date_time
 
@@ -175 +175 @@
 
 !
-!-- Declare all global variables within the module 
+!-- Declare all global variables within the module
 
 !
 !-- 20200203 (ECC) - variable unused
-!    CHARACTER (LEN=80) ::  filename_emis             !< Variable for the name of the netcdf input file
+!    CHARACTER (LEN=80) ::  filename_emis           !< Variable for the name of the netcdf input file
 
 !
 !-- 20200203 (ECC) new variables for on-demand read mode
 
-    CHARACTER(LEN=512), ALLOCATABLE, DIMENSION(:) :: timestamps   !< timestamps in chemistry file
-    
     CHARACTER(LEN=*),   PARAMETER   :: input_file_chem = 'PIDS_CHEM' !< chemistry file
 
-    INTEGER(iwp) ::  dt_emis         !< Time Step Emissions 
+    CHARACTER(LEN=512), ALLOCATABLE, DIMENSION(:) :: timestamps      !< timestamps in chemistry file
+
+
+    INTEGER(iwp) ::  dt_emis         !< Time Step Emissions
     INTEGER(iwp) ::  i               !< index 1st selected dimension (some dims are not spatial)
-    INTEGER(iwp) ::  j               !< index 2nd selected dimension 
+    INTEGER(iwp) ::  j               !< index 2nd selected dimension
+    INTEGER(iwp) ::  i_end           !< Index to end read variable from netcdf in one dims
     INTEGER(iwp) ::  i_start         !< Index to start read variable from netcdf along one dims
-    INTEGER(iwp) ::  i_end           !< Index to end read variable from netcdf in one dims
+    INTEGER(iwp) ::  j_end           !< Index to end read variable from netcdf in additional dims
     INTEGER(iwp) ::  j_start         !< Index to start read variable from netcdf in additional dims
-    INTEGER(iwp) ::  j_end           !< Index to end read variable from netcdf in additional dims
     INTEGER(iwp) ::  len_index       !< length of index (used for several indices)
     INTEGER(iwp) ::  len_index_pm    !< length of PMs index
     INTEGER(iwp) ::  len_index_voc   !< length of voc index
     INTEGER(iwp) ::  previous_timestamp_index !< index for current timestamp (20200203 ECC)
+    INTEGER(iwp) ::  z_end           !< Index to end read variable from netcdf in additional dims
     INTEGER(iwp) ::  z_start         !< Index to start read variable from netcdf in additional dims
-    INTEGER(iwp) ::  z_end           !< Index to end read variable from netcdf in additional dims
 
     REAL(wp) ::  conversion_factor                   !< Units Conversion Factor
@@ -212 +213 @@
 !     END INTERFACE chem_emissions_check_parameters
 !
-!-- Matching Emissions actions  
+!-- Matching Emissions actions
     INTERFACE chem_emissions_match
        MODULE PROCEDURE chem_emissions_match
     END INTERFACE chem_emissions_match
 !
-!-- Initialization actions  
+!-- Initialization actions
     INTERFACE chem_emissions_init
        MODULE PROCEDURE chem_emissions_init
@@ -231 +232 @@
 
 !
-!-- initialization actions for on-demand mode 
+!-- initialization actions for on-demand mode
     INTERFACE chem_emissions_header_init
        MODULE PROCEDURE chem_emissions_header_init
@@ -246 +247 @@
 !    PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
 
-    PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup,    &
+    PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup,                        &
            chem_emissions_header_init, chem_emissions_update_on_demand
 !
@@ -254 +255 @@
  CONTAINS
 
-! !------------------------------------------------------------------------------!
+! !------------------------------------------------------------------------------------------------!
 ! ! Description:
 ! ! ------------
 ! !> Routine for checking input parameters
-! !------------------------------------------------------------------------------!
+! !------------------------------------------------------------------------------------------------!
 !  SUBROUTINE chem_emissions_check_parameters
 !
@@ -281 +282 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Matching the chemical species indices. The routine checks what are the
-!> indices of the emission input species and the corresponding ones of the
-!> model species. The routine gives as output a vector containing the number 
-!> of common species: it is important to note that while the model species
-!> are distinct, their values could be given to a single species in input.
-!> For example, in the case of NO2 and NO, values may be passed in input as
-!> NOX values.
-!------------------------------------------------------------------------------!
-
-SUBROUTINE chem_emissions_match( emt_att,len_index )   
-
-    INTEGER(iwp)  ::  ind_inp                    !< Parameters for cycling through chemical input species 
-    INTEGER(iwp)  ::  ind_mod                    !< Parameters for cycling through chemical model species 
+!> Matching the chemical species indices. The routine checks what are the indices of the emission
+!> input species and the corresponding ones of the model species. The routine gives as output a
+!> vector containing the number of common species: it is important to note that while the model
+!> species are distinct, their values could be given to a single species in input.
+!> For example, in the case of NO2 and NO, values may be passed in input as NOX values.
+!--------------------------------------------------------------------------------------------------!
+
+SUBROUTINE chem_emissions_match( emt_att,len_index )
+
+    INTEGER(iwp)  ::  ind_inp                    !< Parameters for cycling through chemical input species
+    INTEGER(iwp)  ::  ind_mod                    !< Parameters for cycling through chemical model species
     INTEGER(iwp)  ::  ind_voc                    !< Indices to check whether a split for voc should be done
     INTEGER(iwp)  ::  ispec                      !< index for cycle over effective number of emission species
@@ -312 +311 @@
 
     nspec_emis_inp = emt_att%n_emiss_species
-!    nspec_emis_inp=emt_att%nspec 
+!    nspec_emis_inp=emt_att%nspec
 
 !
 !-- Check the emission LOD: 0 (PARAMETERIZED), 1 (DEFAULT), 2 (PRE-PROCESSED)
-!
     SELECT CASE (emiss_lod)
 
 !
-!-- LOD 0 (PARAMETERIZED mode)
-
+!--    LOD 0 (PARAMETERIZED mode)
        CASE (0)
 
           len_index = 0
-
-! number of species and number of matched species can be different
-! but call is only made if both are greater than zero
-
+!
+!--       Number of species and number of matched species can be different but call is only made if
+!--       both are greater than zero.
           IF  ( nvar > 0  .AND.  nspec_emis_inp > 0 )  THEN
 
 !
-!-- Cycle over model species
-
+!--          Cycle over model species
              DO  ind_mod = 1, nvar
                 ind_inp = 1
-                DO  WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )    !< 'novalue' is the default  
+                DO  WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )    !< 'novalue' is the default
                    IF  ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
                       len_index = len_index + 1
@@ -347 +342 @@
 
 !
-!-- Allocation of Arrays of the matched species
-
-                ALLOCATE ( match_spec_input(len_index) ) 
+!--             Allocation of Arrays of the matched species
+                ALLOCATE ( match_spec_input(len_index) )
                 ALLOCATE ( match_spec_model(len_index) )
 
 !
-!-- Pass species indices to declared arrays
-
+!--             Pass species indices to declared arrays
                 len_index = 0
 
-                DO  ind_mod = 1, nvar  
+                DO  ind_mod = 1, nvar
                    ind_inp = 1
                    DO  WHILE  ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )
-                      IF  ( TRIM(surface_csflux_name(ind_inp)) ==          &
-                            TRIM(spc_names(ind_mod))            )  THEN
+                      IF  ( TRIM( surface_csflux_name(ind_inp) ) == TRIM(spc_names(ind_mod) ) )    &
+                      THEN
                          len_index = len_index + 1
                          match_spec_input(len_index) = ind_inp
@@ -371 +364 @@
 
 !
-!-- Check
-
+!--             Check
                 DO  ispec = 1, len_index
 
-                   IF  ( emiss_factor_main(match_spec_input(ispec) ) < 0    .AND.     &
+                   IF  ( emiss_factor_main(match_spec_input(ispec) ) < 0    .AND.                  &
                          emiss_factor_side(match_spec_input(ispec) ) < 0 )  THEN
 
-                      message_string = 'PARAMETERIZED emissions mode selected:'             //          &
-                                       ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES'        //          &
-                                       ' but values of scaling factors for street types'    //          &
-                                       ' emiss_factor_main AND emiss_factor_side'           //          &
-                                       ' not provided for each of the emissions species'    //          &
-                                       ' or not provided at all: PLEASE set a finite value' //          &
-                                       ' for these parameters in the chemistry namelist'          
+                      message_string = 'PARAMETERIZED emissions mode selected:'             //     &
+                                       ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES'        //     &
+                                       ' but values of scaling factors for street types'    //     &
+                                       ' emiss_factor_main AND emiss_factor_side'           //     &
+                                       ' not provided for each of the emissions species'    //     &
+                                       ' or not provided at all: PLEASE set a finite value' //     &
+                                       ' for these parameters in the chemistry namelist'
                       CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 )
- 
+
                    ENDIF
 
@@ -393 +385 @@
 
              ELSE
-                
-                message_string = 'Non of given Emission Species'           //           &
-                                 ' matches'                                //           &
-                                 ' model chemical species'                 //           &
-                                 ' Emission routine is not called'         
+
+                message_string = 'Non of given Emission Species'           //                      &
+                                 ' matches'                                //                      &
+                                 ' model chemical species'                 //                      &
+                                 ' Emission routine is not called'
                 CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 )
- 
+
              ENDIF
 
           ELSE
-     
-             message_string = 'Array of Emission species not allocated: '             //          &
-                              ' Either no emission species are provided as input or'  //          &
-                              ' no chemical species are used by PALM.'                //          &
-                              ' Emission routine is not called'                                   
-             CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 ) 
- 
+
+             message_string = 'Array of Emission species not allocated: '             //           &
+                              ' Either no emission species are provided as input or'  //           &
+                              ' no chemical species are used by PALM.'                //           &
+                              ' Emission routine is not called'
+             CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 )
+
           ENDIF
 
 !
-!-- LOD 1 (DEFAULT mode)
-
+!--    LOD 1 (DEFAULT mode)
        CASE (1)
 
-          len_index = 0          ! total number of species (to be accumulated)  
+          len_index = 0          ! total number of species (to be accumulated)
           len_index_voc = 0      ! total number of VOCs (to be accumulated)
           len_index_pm = 3       ! total number of PMs: PM1, PM2.5, PM10.
 
 !
-!-- number of model species and input species could be different
-!-- but process this only when both are non-zero
- 
+!--       Number of model species and input species could be different but process this only when both are
+!--       non-zero
           IF  ( nvar > 0  .AND.  nspec_emis_inp > 0 )  THEN
 
 !
-!-- Cycle over model species
+!--          Cycle over model species
              DO  ind_mod = 1, nvar
 
 !
-!-- Cycle over input species
-
+!--             Cycle over input species
                 DO  ind_inp = 1, nspec_emis_inp
 
 !
-!-- Check for VOC Species
-
+!--                Check for VOC Species
                    IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN
                       DO  ind_voc= 1, emt_att%nvoc
-                        
-                         IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
+
+                         IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )   &
+                         THEN
                             len_index = len_index + 1
                             len_index_voc = len_index_voc + 1
                          ENDIF
-                         
+
                       END DO
                    ENDIF
 
 !
-!-- PMs: There is one input species name for all PM 
-!-- This variable has 3 dimensions, one for PM1, PM2.5 and PM10
-
+!-- PMs: There is one input species name for all PM. This variable has 3 dimensions, one for PM1,
+!--      PM2.5 and PM10
                    IF  ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
 
@@ -467 +455 @@
 
 !
-!-- NOX: NO2 and NO
-
-                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" )  THEN 
+!--                NOX: NO2 and NO
+                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" )  THEN
 
                       IF     ( TRIM( spc_names(ind_mod) ) == "NO"  )  THEN
@@ -480 +467 @@
 
 !
-!-- SOX: SO2 and SO4
-
+!--                SOX: SO2 and SO4
                    IF  ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" )  THEN
 
@@ -493 +479 @@
 
 !
-!-- Other Species
-
-                   IF  ( TRIM( emt_att%species_name(ind_inp) ) ==             &
-                        TRIM( spc_names(ind_mod) ) )  THEN
+!--                Other Species
+                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
                       len_index = len_index + 1
                    ENDIF
@@ -506 +490 @@
 
 !
-!-- Allocate arrays
-
+!--          Allocate arrays
              IF  ( len_index > 0 )  THEN
 
@@ -514 +497 @@
 
                 IF  ( len_index_voc > 0 )  THEN
-
-!
-!-- Contains indices of the VOC model species
-
+!
+!--                Contains indices of the VOC model species
                    ALLOCATE( match_spec_voc_model(len_index_voc) )
-
-!
-!-- Contains the indices of different values of VOC composition
-!-- of input variable VOC_composition
-
+!
+!--                Contains the indices of different values of VOC composition of input variable
+!--                VOC_composition
                    ALLOCATE( match_spec_voc_input(len_index_voc) )
 
@@ -529 +508 @@
 
 !
-!-- Pass the species indices to declared arrays
-
+!--             Pass the species indices to declared arrays
                 len_index = 0
                 len_index_voc = 0
-                
+
                 DO  ind_mod = 1, nvar
-                   DO  ind_inp = 1, nspec_emis_inp 
-
-!
-!-- VOCs
-
-                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC"  .AND.        &
-                            ALLOCATED (match_spec_voc_input) )  THEN
+                   DO  ind_inp = 1, nspec_emis_inp
+
+!
+!--                   VOCs
+                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC"  .AND.                  &
+                            ALLOCATED( match_spec_voc_input ) )  THEN
 
                          DO  ind_voc = 1, emt_att%nvoc
 
-                            IF  ( TRIM( emt_att%voc_name(ind_voc) ) ==                  &
-                                 TRIM( spc_names(ind_mod) ) )  THEN
+                            IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )&
+                            THEN
 
                                len_index     = len_index + 1
                                len_index_voc = len_index_voc + 1
-                        
+
                                match_spec_input(len_index) = ind_inp
                                match_spec_model(len_index) = ind_mod
@@ -564 +541 @@
 
 !
-!-- PMs
-
+!--                   PMs
                       IF  ( TRIM( emt_att%species_name(ind_inp) ) == "PM" )  THEN
 
@@ -585 +561 @@
 
 !
-!-- NOX 
+!--                   NOX
                       IF  ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" )  THEN
 
@@ -599 +575 @@
                             match_spec_input(len_index) = ind_inp
                             match_spec_model(len_index) = ind_mod
-  
+
                          ENDIF
 
@@ -606 +582 @@
 
 !
-!-- SOX 
-
+!--                   SOX
                       IF  ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" ) THEN
 
@@ -623 +598 @@
 
 !
-!-- Other Species
-
-                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
+!--                   Other Species
+                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  &
+                      THEN
                          len_index = len_index + 1
                          match_spec_input(len_index) = ind_inp
@@ -636 +611 @@
 
 !
-!-- Error reporting (no matching)
-
+!--          Error reporting (no matching)
              ELSE
 
-                message_string = 'None of given Emission Species matches'           //   &
-                                 ' model chemical species'                          //   &
-                                 ' Emission routine is not called'         
-                CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 ) 
+                message_string = 'None of given Emission Species matches'           //             &
+                                 ' model chemical species'                          //             &
+                                 ' Emission routine is not called'
+                CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 )
 
              ENDIF
 
 !
-!-- Error reporting (no species)
-
+!--       Error reporting (no species)
           ELSE
 
-             message_string = 'Array of Emission species not allocated: '             // &
-                              ' Either no emission species are provided as input or'  // &
-                              ' no chemical species are used by PALM:'                // &
-                              ' Emission routine is not called'                                   
-             CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 ) 
- 
+             message_string = 'Array of Emission species not allocated: '             //           &
+                              ' Either no emission species are provided as input or'  //           &
+                              ' no chemical species are used by PALM:'                //           &
+                              ' Emission routine is not called'
+             CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 )
+
           ENDIF
 
 !
-!-- LOD 2 (PRE-PROCESSED mode)
-
+!--    LOD 2 (PRE-PROCESSED mode)
        CASE (2)
 
@@ -670 +642 @@
           IF  ( nvar > 0  .AND.  nspec_emis_inp > 0 )  THEN
 !
-!-- Cycle over model species
-             DO  ind_mod = 1, nvar 
-
-!
-!-- Cycle over input species  
+!--          Cycle over model species
+             DO  ind_mod = 1, nvar
+
+!
+!--             Cycle over input species
                 DO  ind_inp = 1, nspec_emis_inp
 
 !
-!-- Check for VOC Species 
-
-                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN        
+!--                Check for VOC Species
+                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" )  THEN
                       DO  ind_voc = 1, emt_att%nvoc
-                         IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
+                         IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )   &
+                         THEN
                             len_index     = len_index + 1
                             len_index_voc = len_index_voc + 1
@@ -690 +662 @@
 
 !
-!-- Other Species
-
-                   IF  ( TRIM(emt_att%species_name(ind_inp)) == TRIM(spc_names(ind_mod)) )  THEN
+!--                Other Species
+                   IF  ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
                       len_index = len_index + 1
                    ENDIF
@@ -699 +670 @@
 
 !
-!-- Allocate array for storing the indices of the matched species
-
+!--          Allocate array for storing the indices of the matched species
              IF  ( len_index > 0 )  THEN
- 
-                ALLOCATE ( match_spec_input(len_index) ) 
- 
+
+                ALLOCATE ( match_spec_input(len_index) )
+
                 ALLOCATE ( match_spec_model(len_index) )
 
                 IF  ( len_index_voc > 0 )  THEN
 !
-!-- contains indices of the VOC model species
+!--                Contains indices of the VOC model species
                    ALLOCATE( match_spec_voc_model(len_index_voc) )
 !
-!-- contains the indices of different values of VOC composition of input variable VOC_composition
+!--                Contains the indices of different values of VOC composition of input variable
+!--                VOC_composition
                    ALLOCATE( match_spec_voc_input(len_index_voc) )
 
@@ -718 +689 @@
 
 !
-!-- pass the species indices to declared arrays
-
+!--             Pass the species indices to declared arrays
                 len_index = 0
 
 !
-!-- Cycle over model species
-
+!--             Cycle over model species
                 DO  ind_mod = 1, nvar
- 
-!
-!-- Cycle over Input species  
-
+
+!
+!--                Cycle over Input species
                    DO  ind_inp = 1, nspec_emis_inp
 
 !
-!-- VOCs
-
-                      IF  ( TRIM(emt_att%species_name(ind_inp) ) == "VOC"  .AND.     &
-                           ALLOCATED(match_spec_voc_input) )  THEN
-                         
+!--                   VOCs
+                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == "VOC"  .AND.                  &
+                            ALLOCATED( match_spec_voc_input ) )  THEN
+
                          DO  ind_voc= 1, emt_att%nvoc
-                            IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )  THEN
+                            IF  ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )&
+                            THEN
                                len_index = len_index + 1
                                len_index_voc = len_index_voc + 1
-                        
+
                                match_spec_input(len_index) = ind_inp
                                match_spec_model(len_index) = ind_mod
 
                                match_spec_voc_input(len_index_voc) = ind_voc
-                               match_spec_voc_model(len_index_voc) = ind_mod                         
+                               match_spec_voc_model(len_index_voc) = ind_mod
                             ENDIF
                          END DO
@@ -753 +721 @@
 
 !
-!-- Other Species
-
-                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  THEN
+!--                   Other Species
+                      IF  ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) )  &
+                      THEN
                          len_index = len_index + 1
                          match_spec_input(len_index) = ind_inp
@@ -767 +735 @@
 
 !
-!-- in case there are no species matching (just informational message)
-
-                message_string = 'Non of given emission species'            //         &
-                                 ' matches'                                //          &
-                                 ' model chemical species:'                //          &
-                                 ' Emission routine is not called'  
+!--             In case there are no species matching (just informational message)
+                message_string = 'Non of given emission species'            //                     &
+                                 ' matches'                                //                      &
+                                 ' model chemical species:'                //                      &
+                                 ' Emission routine is not called'
                 CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
              ENDIF
 
 !
-!-- Error check (no matching)
- 
+!--       Error check (no matching)
           ELSE
 
 !
-!-- either spc_names is zero or nspec_emis_inp is not allocated
-             message_string = 'Array of Emission species not allocated:'              // &
-                              ' Either no emission species are provided as input or'  // &
-                              ' no chemical species are used by PALM:'                // &
-                              ' Emission routine is not called'                  
-             CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 ) 
-
-          ENDIF  
+!--          Either spc_names is zero or nspec_emis_inp is not allocated
+             message_string = 'Array of Emission species not allocated:'              //           &
+                              ' Either no emission species are provided as input or'  //           &
+                              ' no chemical species are used by PALM:'                //           &
+                              ' Emission routine is not called'
+             CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 )
+
+          ENDIF
 
 !
@@ -795 +761 @@
 
 !
-!-- Error check (no species)
-
+!--    Error check (no species)
        CASE DEFAULT
 
           message_string = 'Emission Module switched ON, but'                           //         &
                            ' either no emission mode specified or incorrectly given :'  //         &
-                           ' please, pass the correct value to the namelist parameter "mode_emis"'                  
-          CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )              
+                           ' please, pass the correct value to the namelist parameter "mode_emis"'
+          CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )
 
        END SELECT
@@ -810 +775 @@
  END SUBROUTINE chem_emissions_match
 
- 
-!------------------------------------------------------------------------------!
+
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Initialization:
-!> Netcdf reading, arrays allocation and first calculation of cssws
-!> fluxes at timestep 0 
-!------------------------------------------------------------------------------!
+!> Netcdf reading, arrays allocation and first calculation of cssws fluxes at timestep 0
+!--------------------------------------------------------------------------------------------------!
 
  SUBROUTINE chem_emissions_init
 
-    USE netcdf_data_input_mod,                                              &
+    USE netcdf_data_input_mod,                                                                     &
         ONLY:  chem_emis, chem_emis_att
-    
+
     IMPLICIT NONE
- 
+
     INTEGER(iwp) :: ispec                        !< running index
 
-!    
-!-- Actions for initial runs 
+!
+!-- Actions for initial runs
 !  IF  (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
-!--    ... 
-!    
+!--    ...
+!
 !
 !-- Actions for restart runs
@@ -845 +809 @@
 !
 !-- Matching
-
     CALL  chem_emissions_match( chem_emis_att, n_matched_vars )
- 
+
     IF  ( n_matched_vars == 0 )  THEN
-  
+
        emission_output_required = .FALSE.
 
@@ -858 +821 @@
 
 !
-!-- Set molecule masses  (in kg/mol)
-
+!--    Set molecule masses  (in kg/mol)
        ALLOCATE( chem_emis_att%xm(n_matched_vars) )
 
@@ -878 +840 @@
        ENDDO
 
-    
-!
-!-- Get emissions for the first time step base on LOD (if defined)
-!-- or emission mode (if no LOD defined)
-
-!
-!-- NOTE - I could use a combined if ( lod = xxx .or. mode = 'XXX' )
-!--        type of decision structure but I think it is much better
-!--        to implement it this way (i.e., conditional on lod if it
-!--        is defined, and mode if not) as we can easily take out
-!--        the case structure for mode_emis later on.
+
+!
+!-- Get emissions for the first time step base on LOD (if defined) or emission mode
+!-- (if no LOD defined)
+
+!
+!-- NOTE - I could use a combined if ( lod = xxx .or. mode = 'XXX' ) type of decision structure but
+!          I think it is much better to implement it this way (i.e., conditional on lod if it is
+!          defined, and mode if not) as we can easily take out the case structure for mode_emis
+!          later on.
 
        IF   ( emiss_lod < 0 )  THEN  !-- no LOD defined (not likely)
 
-          SELECT CASE ( TRIM( mode_emis ) )   
+          SELECT CASE ( TRIM( mode_emis ) )
 
              CASE ( 'PARAMETERIZED' )     ! LOD 0
 
-                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
                 ENDIF
@@ -904 +865 @@
              CASE ( 'DEFAULT' )           ! LOD 1
 
-                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
                 ENDIF
@@ -912 +873 @@
              CASE ( 'PRE-PROCESSED' )     ! LOD 2
 
-                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
-!
-!--                Note, at the moment emissions are considered only by surface fluxes
-!--                rather than by volume sources. Therefore, no vertical dimension is 
-!--                required and is thus allocated with 1. Later when volume sources
-!--                are considered, the vertical dimension will increase.                 
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
+!
+!--                Note, at the moment emissions are considered only by surface fluxes rather than
+!--                by volume sources. Therefore, no vertical dimension is required and is thus
+!--                allocated with 1. Later when volume sources are considered, the vertical
+!--                dimension will increase.
                    !ALLOCATE( emis_distribution(nzb:nzt+1,nys:nyn,nxl:nxr,n_matched_vars) )
                    ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
                 ENDIF
- 
+
                 CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
 
@@ -932 +893 @@
              CASE ( 0 )     ! parameterized mode
 
-                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
                 ENDIF
@@ -940 +901 @@
              CASE ( 1 )     ! default mode
 
-                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
                 ENDIF
@@ -948 +909 @@
              CASE ( 2 )     ! pre-processed mode
 
-                IF  (  .NOT.  ALLOCATED( emis_distribution) )  THEN
+                IF  ( .NOT. ALLOCATED( emis_distribution) )  THEN
                    ALLOCATE( emis_distribution(nzb:nzt+1,nys:nyn,nxl:nxr,n_matched_vars) )
                 ENDIF
- 
+
                 CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
 
@@ -959 +920 @@
 
 !
-! -- initialize
-
+! --   Initialize
        emis_distribution = 0.0_wp
 
@@ -971 +931 @@
 
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Routine for Update of Emission values at each timestep.
 !>
-!> @todo Clarify the correct usage of index_dd, index_hh and index_mm. Consider
-!>       renaming of these variables.
+!> @todo Clarify the correct usage of index_dd, index_hh and index_mm. Consider renaming of these
+!>       variables.
 !> @todo Clarify time used in emis_lod=2 mode. ATM, the used time seems strange.
-!-------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 
  SUBROUTINE chem_emissions_setup( emt_att, emt, n_matched_vars )
- 
-   USE surface_mod,                                               &
+
+   USE surface_mod,                                                                                &
        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
 
-   USE netcdf_data_input_mod,                                     &
+   USE netcdf_data_input_mod,                                                                      &
        ONLY:  street_type_f
 
-   USE arrays_3d,                                                 &        
-       ONLY: hyp, pt 
-
-    USE control_parameters, &
+   USE arrays_3d,                                                                                  &
+       ONLY: hyp, pt
+
+    USE control_parameters,                                                                        &
         ONLY: time_since_reference_point
 
-    USE palm_date_time_mod, &
+    USE palm_date_time_mod,                                                                        &
         ONLY: days_per_week, get_date_time, hours_per_day, months_per_year, seconds_per_day
-   
+
     IMPLICIT NONE
-  
-
-    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att                         !< variable to store emission information 
-
-    TYPE(chem_emis_val_type), INTENT(INOUT), ALLOCATABLE, DIMENSION(:) ::  emt  !< variable to store emission input values, 
+
+    INTEGER(iwp) ::  day_of_month                          !< day of the month
+    INTEGER(iwp) ::  day_of_week                           !< day of the week
+    INTEGER(iwp) ::  day_of_year                           !< day of the year
+    INTEGER(iwp) ::  days_since_reference_point            !< days since reference point
+    INTEGER(iwp) ::  i                                     !< running index for grid in x-direction
+    INTEGER(iwp) ::  i_pm_comp                             !< index for number of PM components
+    INTEGER(iwp) ::  icat                                  !< Index for number of categories
+    INTEGER(iwp) ::  index_dd                              !< index day
+    INTEGER(iwp) ::  index_hh                              !< index hour
+    INTEGER(iwp) ::  index_mm                              !< index month
+    INTEGER(iwp) ::  ispec                                 !< index for number of species
+    INTEGER(iwp) ::  ivoc                                  !< Index for number of VOCs
+    INTEGER(iwp) ::  hour_of_day                           !< hour of the day
+    INTEGER(iwp) ::  j                                     !< running index for grid in y-direction
+    INTEGER(iwp) ::  k                                     !< running index for grid in z-direction
+    INTEGER(iwp) ::  m                                     !< running index for horizontal surfaces
+    INTEGER(iwp) ::  month_of_year                         !< month of the year
+
+    INTEGER,INTENT(IN) ::  n_matched_vars                  !< Output of matching routine with number
+                                                           !< of matched species
+
+    REAL(wp) ::  time_utc_init                             !< second of day of initial date
+
+    TYPE(chem_emis_att_type), INTENT(INOUT) ::  emt_att    !< variable to store emission information
+
+    TYPE(chem_emis_val_type), INTENT(INOUT), ALLOCATABLE, DIMENSION(:) ::  emt  !< variable to store emission input values,
                                                                                 !< depending on the considered species
-
-    INTEGER,INTENT(IN) ::  n_matched_vars                                       !< Output of matching routine with number
-                                                                                !< of matched species
-
-    INTEGER(iwp) ::  i                                                          !< running index for grid in x-direction
-    INTEGER(iwp) ::  i_pm_comp                                                  !< index for number of PM components
-    INTEGER(iwp) ::  icat                                                       !< Index for number of categories
-    INTEGER(iwp) ::  ispec                                                      !< index for number of species
-    INTEGER(iwp) ::  ivoc                                                       !< Index for number of VOCs
-    INTEGER(iwp) ::  j                                                          !< running index for grid in y-direction
-    INTEGER(iwp) ::  k                                                          !< running index for grid in z-direction
-    INTEGER(iwp) ::  m                                                          !< running index for horizontal surfaces
-
-    INTEGER(iwp) ::  day_of_month                                               !< day of the month
-    INTEGER(iwp) ::  day_of_week                                                !< day of the week
-    INTEGER(iwp) ::  day_of_year                                                !< day of the year
-    INTEGER(iwp) ::  days_since_reference_point                                 !< days since reference point
-    INTEGER(iwp) ::  hour_of_day                                                !< hour of the day
-    INTEGER(iwp) ::  month_of_year                                              !< month of the year
-    INTEGER(iwp) ::  index_dd                                                   !< index day
-    INTEGER(iwp) ::  index_hh                                                   !< index hour
-    INTEGER(iwp) ::  index_mm                                                   !< index month
-
-    REAL(wp) ::  time_utc_init                                                  !< second of day of initial date
-
-    !
-    !-- CONVERSION FACTORS: TIME  
-    REAL(wp), PARAMETER ::  hour_per_year =  8760.0_wp  !< number of hours in a year of 365 days  
-    REAL(wp), PARAMETER ::  s_per_hour    =  3600.0_wp  !< number of sec per hour (s)/(hour)    
-    REAL(wp), PARAMETER ::  s_per_day     = 86400.0_wp  !< number of sec per day (s)/(day)  
+!
+!-- CONVERSION FACTORS: TIME
+    REAL(wp), PARAMETER ::  hour_per_year =  8760.0_wp                         !< number of hours in a year of 365 days
+    REAL(wp), PARAMETER ::  s_per_hour    =  3600.0_wp                         !< number of sec per hour (s)/(hour)
+    REAL(wp), PARAMETER ::  s_per_day     = 86400.0_wp                         !< number of sec per day (s)/(day)
 
     REAL(wp), PARAMETER ::  day_to_s      = 1.0_wp/s_per_day                   !< conversion day   -> sec
     REAL(wp), PARAMETER ::  hour_to_s     = 1.0_wp/s_per_hour                  !< conversion hours -> sec
     REAL(wp), PARAMETER ::  year_to_s     = 1.0_wp/(s_per_hour*hour_per_year)  !< conversion year  -> sec
-    !
-    !-- CONVERSION FACTORS: MASS
-    REAL(wp), PARAMETER ::  g_to_kg       = 1.0E-03_wp     !< Conversion from g to kg (kg/g) 
+
+
+!
+!-- CONVERSION FACTORS: MASS
+    REAL(wp), PARAMETER ::  g_to_kg       = 1.0E-03_wp     !< Conversion from g to kg (kg/g)
     REAL(wp), PARAMETER ::  miug_to_kg    = 1.0E-09_wp     !< Conversion from g to kg (kg/g)
-    REAL(wp), PARAMETER ::  tons_to_kg    = 100.0_wp       !< Conversion from tons to kg (kg/tons)    
-    !
-    !-- CONVERSION FACTORS: PPM
-    REAL(wp), PARAMETER ::  ratio2ppm     = 1.0E+06_wp 
- 
+    REAL(wp), PARAMETER ::  tons_to_kg    = 100.0_wp       !< Conversion from tons to kg (kg/tons)
+!
+!-- CONVERSION FACTORS: PPM
+    REAL(wp), PARAMETER ::  ratio2ppm     = 1.0E+06_wp
+
     REAL(wp), DIMENSION(24)                        ::  par_emis_time_factor    !< time factors for the parameterized mode:
                                                                                !< fixed houlry profile for example day
-    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  conv_to_ratio           !< factor used for converting input 
+    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  conv_to_ratio           !< factor used for converting input
                                                                                !< to concentration ratio
     REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  tmp_temp                !< temporary variable for abs. temperature
 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  delta_emis   !< incremental emission factor                     
     REAL(wp), DIMENSION(:),   ALLOCATABLE ::  time_factor  !< factor for time scaling of emissions
+
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  delta_emis   !< incremental emission factor
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  emis         !< emission factor
 
+
     IF  ( emission_output_required )  THEN
 
 !
-!-- Set emis_dt to be used - since chemistry ODEs can be stiff, the option
-!-- to solve them at every RK substep is present to help improve stability
-!-- should the need arises
- 
+!-- Set emis_dt to be used - since chemistry ODEs can be stiff, the option to solve them at every
+!-- RK substep is present to help improve stability should the need arises
+
        IF  ( call_chem_at_all_substeps )  THEN
 
@@ -1076 +1036 @@
 
  !
- !-- Conversion of units to the ones employed in PALM 
- !-- In PARAMETERIZED mode no conversion is performed: in this case input units are fixed
-
+ !--    Conversion of units to the ones employed in PALM
+ !--    In PARAMETERIZED mode no conversion is performed: in this case input units are fixed
         IF  ( TRIM( mode_emis ) == "DEFAULT"  .OR.  TRIM( mode_emis ) == "PRE-PROCESSED" )  THEN
 
@@ -1101 +1060 @@
 
 !
-!-- Error check (need units)
-
-             CASE DEFAULT  
-                message_string = 'The Units of the provided emission input'                 // &
-                                 ' are not the ones required by PALM-4U: please check '     // &
-                                 ' emission module documentation.'                                 
+!--          Error check (need units)
+             CASE DEFAULT
+                message_string = 'The Units of the provided emission input'                 //     &
+                                 ' are not the ones required by PALM-4U: please check '     //     &
+                                 ' emission module documentation.'
                 CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 )
 
@@ -1114 +1072 @@
 
 !
-!-- Conversion factor to convert  kg/m**2/s to ppm/s
-
+!--    Conversion factor to convert  kg/m**2/s to ppm/s
        DO  i = nxl, nxr
           DO  j = nys, nyn
 
 !
-!-- Derive Temperature from Potential Temperature
-
-             tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) *                     &
-                                       ( hyp(nzb:nzt+1) / p_0 )**rd_d_cp
- 
-!            
-!-- We need to pass to cssws <- (ppm/s) * dz 
-!-- Input is Nmole/(m^2*s)
-!-- To go to ppm*dz multiply the input by (m**2/N)*dz
-!-- (m**2/N)*dz == V/N 
-!-- V/N = RT/P 
-
-             conv_to_ratio(nzb:nzt+1,j,i) =  rgas_univ *                       &  ! J K-1 mol-1
-                                             tmp_temp(nzb:nzt+1,j,i) /         &  ! K
-                                             hyp(nzb:nzt+1)                       ! Pa
+!--          Derive Temperature from Potential Temperature
+             tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) * ( hyp(nzb:nzt+1) / p_0 )**rd_d_cp
+
+!
+!--          We need to pass to cssws <- (ppm/s) * dz
+!--          Input is Nmole/(m^2*s)
+!--          To go to ppm*dz multiply the input by (m**2/N)*dz
+!--          (m**2/N)*dz == V/N
+!--          V/N = RT/P
+             conv_to_ratio(nzb:nzt+1,j,i) =  rgas_univ *                            &  ! J K-1 mol-1
+                                             tmp_temp(nzb:nzt+1,j,i) /              &  ! K
+                                             hyp(nzb:nzt+1)                            ! Pa
 
 ! (ecc) for reference
-!                   m**3/Nmole               (J/mol)*K^-1           K                      Pa         
-!             conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) )  
+!                   m**3/Nmole               (J/mol)*K^-1           K                      Pa
+!             conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) )
 
           ENDDO
@@ -1150 +1104 @@
 !       emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp
 
-  
+
 !
 !-- LOD 2 (PRE-PROCESSED MODE)
@@ -1160 +1114 @@
 
 !
-!-- Update time indices
-
+!--       Update time indices
           CALL get_date_time( 0.0_wp, second_of_day=time_utc_init )
-          CALL get_date_time( MAX( 0.0_wp, time_since_reference_point ),       &
-                              hour=hour_of_day )
-
-          days_since_reference_point = INT( FLOOR( (                            &
-                    time_utc_init + MAX( 0.0_wp, time_since_reference_point ) ) &
-                                                   / seconds_per_day ) )
+          CALL get_date_time( MAX( 0.0_wp, time_since_reference_point ), hour=hour_of_day )
+
+          days_since_reference_point = INT( FLOOR( ( time_utc_init +                               &
+                                                     MAX( 0.0_wp, time_since_reference_point ) )   &
+                                            / seconds_per_day ) )
 
           index_hh = days_since_reference_point * hours_per_day + hour_of_day
 
 !
-!-- LOD 1 (DEFAULT MODE)
-
+!--     LOD 1 (DEFAULT MODE)
         ELSEIF  ( emiss_lod == 1 )  THEN
 
@@ -1181 +1132 @@
 
 !
-!-- Allocate array where to store temporary emission values
-
-          IF  (  .NOT.  ALLOCATED(emis) ) ALLOCATE( emis(nys:nyn,nxl:nxr) )
-
-!
-!-- Allocate time factor per category
-
+!--       Allocate array where to store temporary emission values
+          IF  ( .NOT. ALLOCATED(emis) ) ALLOCATE( emis(nys:nyn,nxl:nxr) )
+
+!
+!--       Allocate time factor per category
           ALLOCATE( time_factor(emt_att%ncat) )
 
 !
-!-- Read-in hourly emission time factor
-
-          IF  ( TRIM(time_fac_type) == "HOUR" )  THEN
-
-!
-!-- Update time indices
-
-             CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ),    &
+!--       Read-in hourly emission time factor
+          IF  ( TRIM( time_fac_type ) == "HOUR" )  THEN
+
+!
+!--          Update time indices
+             CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ),                        &
                                  day_of_year=day_of_year, hour=hour_of_day )
              index_hh = ( day_of_year - 1_iwp ) * hour_of_day
 
 !
-!-- Check if the index is less or equal to the temporal dimension of HOURLY emission files
-
+!--          Check if the index is less or equal to the temporal dimension of HOURLY emission files
              IF  ( index_hh <= SIZE( emt_att%hourly_emis_time_factor(1,:) ) )  THEN
 
 !
-!-- Read-in the correspondant time factor
-
-                time_factor(:) = emt_att%hourly_emis_time_factor(:,index_hh+1)      
-
-!
-!-- Error check (time out of range)
-
+!--             Read-in the correspondant time factor
+                time_factor(:) = emt_att%hourly_emis_time_factor(:,index_hh+1)
+
+!
+!--          Error check (time out of range)
              ELSE
 
-                message_string = 'The "HOUR" time-factors in the DEFAULT mode '           // &
-                              ' are not provided for each hour of the total simulation time'      
-                CALL message( 'chem_emissions_setup', 'CM0448', 2, 2, 0, 6, 0 ) 
+                message_string = 'The "HOUR" time-factors in the DEFAULT mode '           //       &
+                                 ' are not provided for each hour of the total simulation time'
+                CALL message( 'chem_emissions_setup', 'CM0448', 2, 2, 0, 6, 0 )
 
              ENDIF
 
 !
-!-- Read-in MDH emissions time factors
-
+!--       Read-in MDH emissions time factors
           ELSEIF  ( TRIM( time_fac_type ) == "MDH" )  THEN
 
 !
-!-- Update time indices
-             CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), &
-                                 month=month_of_year,                       &
-                                 day=day_of_month,                          &
-                                 hour=hour_of_day,                          &
-                                 day_of_week=day_of_week     )
+!--          Update time indices
+             CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ),                        &
+                                 month       = month_of_year,                                      &
+                                 day         = day_of_month,                                       &
+                                 hour        = hour_of_day,                                        &
+                                 day_of_week = day_of_week                                         &
+                               )
              index_mm = month_of_year
              index_dd = months_per_year + day_of_week
-             SELECT CASE(TRIM(daytype_mdh))
+             SELECT CASE( TRIM( daytype_mdh ) )
 
                 CASE ("workday")
@@ -1249 +1193 @@
 
              END SELECT
-!
-!-- Check if the index is less or equal to the temporal dimension of MDH emission files
-
-             IF  ( ( index_hh + index_dd + index_mm) <= SIZE( emt_att%mdh_emis_time_factor(1,:) ) )  THEN
-
-!
-!-- Read in corresponding time factor
-
-                time_factor(:) = emt_att%mdh_emis_time_factor(:,index_mm) *    &
-                                 emt_att%mdh_emis_time_factor(:,index_dd) *    &
+
+!
+!--          Check if the index is less or equal to the temporal dimension of MDH emission files
+             IF  ( ( index_hh + index_dd + index_mm) <= SIZE( emt_att%mdh_emis_time_factor(1,:) ) )&
+             THEN
+!
+!--             Read in corresponding time factor
+                time_factor(:) = emt_att%mdh_emis_time_factor(:,index_mm) *                        &
+                                 emt_att%mdh_emis_time_factor(:,index_dd) *                        &
                                  emt_att%mdh_emis_time_factor(:,index_hh+1)
 
 !
-!-- Error check (MDH time factor not provided)
-     
+!--          Error check (MDH time factor not provided)
              ELSE
 
-                message_string = 'The "MDH" time-factors in the DEFAULT mode '           //  &
-                              ' are not provided for each hour/day/month  of the total simulation time'      
+                message_string = 'The "MDH" time-factors in the DEFAULT mode '           //        &
+                           ' are not provided for each hour/day/month  of the total simulation time'
                 CALL message( 'chem_emissions_setup', 'CM0449', 2, 2, 0, 6, 0 )
 
-             ENDIF  
-
-!
-!-- Error check (no time factor defined)
-
+             ENDIF
+
+!
+!--       Error check (no time factor defined)
           ELSE
-                 
-             message_string = 'In the DEFAULT mode the time factor'           //  &
-                              ' has to be defined in the NAMELIST'      
+
+             message_string = 'In the DEFAULT mode the time factor'           //                   &
+                              ' has to be defined in the NAMELIST'
              CALL message( 'chem_emissions_setup', 'CM0450', 2, 2, 0, 6, 0 )
-         
+
           ENDIF
 
 !
-!-- PARAMETERIZED MODE
-
+!--    PARAMETERIZED MODE
        ELSEIF ( emiss_lod == 0 )  THEN
 
@@ -1291 +1231 @@
 ! for reference (ecc)
 !       ELSEIF  ( TRIM( mode_emis ) == "PARAMETERIZED" )  THEN
-        
-!
-!-- assign constant values of time factors, diurnal profile for traffic sector
-
-          par_emis_time_factor( : ) =  (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039,   &
-                                          0.056, 0.053, 0.051, 0.051, 0.052, 0.055,   &
-                                          0.059, 0.061, 0.064, 0.067, 0.069, 0.069,   &
+
+!
+!--       Assign constant values of time factors, diurnal profile for traffic sector
+          par_emis_time_factor(:) =    (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039,                &
+                                          0.056, 0.053, 0.051, 0.051, 0.052, 0.055,                &
+                                          0.059, 0.061, 0.064, 0.067, 0.069, 0.069,                &
                                           0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /)
-          
-          IF  (  .NOT.  ALLOCATED (time_factor) )  ALLOCATE (time_factor(1))
+
+          IF  ( .NOT. ALLOCATED( time_factor ) )  ALLOCATE( time_factor(1) )
 
 !
 !--       Get time-factor for specific hour
-          CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ),              &
-                              hour=hour_of_day )
+          CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), hour = hour_of_day )
 
           index_hh = hour_of_day
@@ -1317 +1255 @@
 
 !
-!-- LOD 0 (PARAMETERIZED mode)
-
+!--    LOD 0 (PARAMETERIZED mode)
        IF  ( emiss_lod == 0 )  THEN
 
@@ -1328 +1265 @@
 !
 !-- Units are micromoles/m**2*day (or kilograms/m**2*day for PMs)
-
-             emis_distribution(1,nys:nyn,nxl:nxr,ispec) =            &
-                       surface_csflux(match_spec_input(ispec)) *     &
-                       time_factor(1) * hour_to_s
-
-          ENDDO  
-
-          
-!
-!-- LOD 1 (DEFAULT mode)
-
-
+             emis_distribution(1,nys:nyn,nxl:nxr,ispec) = surface_csflux( match_spec_input(ispec) )&
+                                                          * time_factor(1) * hour_to_s
+
+          ENDDO
+
+
+!
+!--    LOD 1 (DEFAULT mode)
        ELSEIF  ( emiss_lod == 1 )  THEN
 
@@ -1346 +1279 @@
 
 !
-!-- Allocate array for the emission value corresponding to a specific category and time factor
-
+!--       Allocate array for the emission value corresponding to a specific category and time factor
           ALLOCATE (delta_emis(nys:nyn,nxl:nxr))
 
 !
-!-- Cycle over categories
-
+!--       Cycle over categories
           DO  icat = 1, emt_att%ncat
- 
-!
-!-- Cycle over Species:  n_matched_vars represents the number of species
-!-- in common between the emission input data and the chemistry mechanism used
-
+
+!
+!--          Cycle over Species: n_matched_vars represents the number of species in common between
+!--                              the emission input data and the chemistry mechanism used
              DO  ispec = 1, n_matched_vars
 
-                emis(nys:nyn,nxl:nxr) =                    &
-                       emt(match_spec_input(ispec))%       &
-                           default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1)
-
-!
-!-- NO 
-
-                IF  ( TRIM(spc_names(match_spec_model(ispec))) == "NO" )  THEN
-                
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
-                                                 time_factor(icat) *           &
-                                                 emt_att%nox_comp(icat,1) *    &
+                emis(nys:nyn,nxl:nxr) = emt( match_spec_input(ispec) )%                            &
+                                        default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1)
+
+!
+!--             NO
+                IF  ( TRIM( spc_names( match_spec_model(ispec) ) ) == "NO" )  THEN
+
+                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *         &         ! kg/m2/s
+                                                 time_factor(icat) *                               &
+                                                 emt_att%nox_comp(icat,1) *                        &
                                                  conversion_factor * hours_per_day
 
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
-                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
-                               delta_emis(nys:nyn,nxl:nxr)
-!
-!-- NO2
-
-                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "NO2" )  THEN
-
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
-                                                 time_factor(icat) *           &
-                                                 emt_att%nox_comp(icat,2) *    &
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                    &
+                                                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)&
+                                                         + delta_emis(nys:nyn,nxl:nxr)
+!
+!--             NO2
+                ELSEIF  ( TRIM( spc_names( match_spec_model(ispec) ) ) == "NO2" )  THEN
+
+                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *         &         ! kg/m2/s
+                                                 time_factor(icat) *                               &
+                                                 emt_att%nox_comp(icat,2) *                        &
                                                  conversion_factor * hours_per_day
 
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
-                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
-                               delta_emis(nys:nyn,nxl:nxr)
- 
-!
-!-- SO2
-                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "SO2" )  THEN
-                  
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
-                                                 time_factor(icat) *           &
-                                                 emt_att%sox_comp(icat,1) *    &
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                    &
+                                                        emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
+                                                        + delta_emis(nys:nyn,nxl:nxr)
+
+!
+!--             SO2
+                ELSEIF  ( TRIM( spc_names( match_spec_model(ispec) ) ) == "SO2" )  THEN
+
+                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *         &         ! kg/m2/s
+                                                 time_factor(icat) *                               &
+                                                 emt_att%sox_comp(icat,1) *                        &
                                                  conversion_factor * hours_per_day
 
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                &
-                               emis_distribution(1,nys:nyn,nxl:nxr,ispec) +    &
-                               delta_emis(nys:nyn,nxl:nxr)
-
-!
-!-- SO4
-
-                ELSEIF  ( TRIM(spc_names(match_spec_model(ispec))) == "SO4" )  THEN
-
-                  
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) *       &         ! kg/m2/s
-                                                 time_factor(icat) *           &
-                                                 emt_att%sox_comp(icat,2) *    &