Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 1787)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 1788)
@@ -19,6 +19,8 @@
 ! Current revisions:
 ! -----------------
-!
-!
+! Added check for use of most_method = 'lookup' in combination with water
+! surface presribed in the land surface model. Added output of z0q.
+! Syntax layout improved.
+! 
 ! Former revisions:
 ! -----------------
@@ -397,5 +399,5 @@
 
        IF ( dt_coupling == 9999999.9_wp )  THEN
-          message_string = 'dt_coupling is not set but required for coup' // &
+          message_string = 'dt_coupling is not set but required for coup' //   &
                            'ling mode "' //  TRIM( coupling_mode ) // '"'
           CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
@@ -419,7 +421,7 @@
 #if ! defined( __check )
        IF ( myid == 0 ) THEN
-          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &
+          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, &
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter,       &
                          status, ierr )
        ENDIF
@@ -436,5 +438,5 @@
           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, &
+             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter,    &
                             status, ierr )
           ENDIF   
@@ -451,5 +453,5 @@
           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, &
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter,       &
                          status, ierr )
        ENDIF 
@@ -464,7 +466,7 @@
 #if ! defined( __check )
        IF ( myid == 0 ) THEN
-          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &
+          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, &
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter,       &
                          status, ierr )
        ENDIF   
@@ -483,5 +485,5 @@
           CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
                          14, comm_inter, ierr )
-          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter,       &
                          status, ierr )    
        ENDIF 
@@ -499,5 +501,5 @@
        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, &
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter,       &
                                                              status, ierr )
        ENDIF
@@ -508,6 +510,6 @@
 
           IF ( dx < remote ) THEN
-             WRITE( message_string, * ) 'coupling mode "', &
-                   TRIM( coupling_mode ),                  &
+             WRITE( message_string, * ) 'coupling mode "',                     &
+                   TRIM( coupling_mode ),                                      &
            '": dx in Atmosphere is not equal to or not larger then dx in ocean'
              CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
@@ -515,6 +517,6 @@
 
           IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
-             WRITE( message_string, * ) 'coupling mode "', &
-                    TRIM( coupling_mode ), &
+             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 )
@@ -526,5 +528,5 @@
        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, &
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter,       &
                          status, ierr )
        ENDIF
@@ -534,6 +536,6 @@
 
           IF ( dy < remote )  THEN
-             WRITE( message_string, * ) 'coupling mode "', &
-                    TRIM( coupling_mode ), &
+             WRITE( message_string, * ) 'coupling mode "',                     &
+                    TRIM( coupling_mode ),                                     &
                  '": dy in Atmosphere is not equal to or not larger then dy in ocean'
              CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
@@ -541,6 +543,6 @@
 
           IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
-             WRITE( message_string, * ) 'coupling mode "', &
-                   TRIM( coupling_mode ), &
+             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 )
@@ -548,6 +550,6 @@
 
           IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
-             WRITE( message_string, * ) 'coupling mode "', &
-                   TRIM( coupling_mode ), &
+             WRITE( message_string, * ) 'coupling mode "',                     &
+                   TRIM( coupling_mode ),                                      &
              '": nx+1 in ocean is not divisible without remainder with nx+1 in', & 
              ' atmosphere'
@@ -556,6 +558,6 @@
 
           IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
-             WRITE( message_string, * ) 'coupling mode "', &
-                   TRIM( coupling_mode ), &
+             WRITE( message_string, * ) 'coupling mode "',                     &
+                   TRIM( coupling_mode ),                                      &
              '": ny+1 in ocean is not divisible without remainder with ny+1 in', & 
              ' atmosphere'
@@ -565,5 +567,5 @@
        ENDIF
 #else
-       WRITE( message_string, * ) 'coupling requires PALM to be called with', &
+       WRITE( message_string, * ) 'coupling requires PALM to be called with',  &
             ' ''mrun -K parallel'''
        CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
@@ -575,8 +577,8 @@
 !-- 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, &
+       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, &
+       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19,          &
                       comm_inter, status, ierr )
     ENDIF
@@ -633,5 +635,5 @@
 
        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 )
@@ -659,9 +661,9 @@
     IF ( topography /= 'flat' )  THEN
        action = ' '
-       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme'     &
+       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme'      &
       .AND. scalar_advec /= 'ws-scheme-mono' )  THEN
           WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
        ENDIF
-       IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' ) &
+       IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' )&
        THEN
           WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
@@ -686,5 +688,5 @@
        ENDIF
        IF ( action /= ' ' )  THEN
-          message_string = 'a non-flat topography does not allow ' // &
+          message_string = 'a non-flat topography does not allow ' //          &
                            TRIM( action )
           CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
@@ -695,15 +697,15 @@
 !--    is applicable. If this is not possible, abort.
        IF ( TRIM( topography_grid_convention ) == ' ' )  THEN
-          IF ( TRIM( topography ) /= 'single_building' .AND.  &
-               TRIM( topography ) /= 'single_street_canyon' .AND.  &
+          IF ( TRIM( topography ) /= 'single_building' .AND.                   &
+               TRIM( topography ) /= 'single_street_canyon' .AND.              &
                TRIM( topography ) /= 'read_from_file' )  THEN
 !--          The default value is not applicable here, because it is only valid
 !--          for the two standard cases 'single_building' and 'read_from_file'
 !--          defined in init_grid.
-             WRITE( message_string, * )  &
-                  'The value for "topography_grid_convention" ',  &
-                  'is not set. Its default value is & only valid for ',  &
-                  '"topography" = ''single_building'', ',  &
-                  '''single_street_canyon'' & or ''read_from_file''.',  &
+             WRITE( message_string, * )                                        &
+                  'The value for "topography_grid_convention" ',               &
+                  'is not set. Its default value is & only valid for ',        &
+                  '"topography" = ''single_building'', ',                      &
+                  '''single_street_canyon'' & or ''read_from_file''.',         &
                   ' & Choose ''cell_edge'' or ''cell_center''.'
              CALL message( 'user_check_parameters', 'PA0239', 1, 2, 0, 6, 0 )
@@ -711,5 +713,5 @@
 !--          The default value is applicable here.
 !--          Set convention according to topography.
-             IF ( TRIM( topography ) == 'single_building' .OR.  &
+             IF ( TRIM( topography ) == 'single_building' .OR.                 &
                   TRIM( topography ) == 'single_street_canyon' )  THEN
                 topography_grid_convention = 'cell_edge'
@@ -718,8 +720,8 @@
              ENDIF
           ENDIF
-       ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND.  &
+       ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND.        &
                 TRIM( topography_grid_convention ) /= 'cell_center' )  THEN
-          WRITE( message_string, * )  &
-               'The value for "topography_grid_convention" is ', &
+          WRITE( message_string, * )                                           &
+               'The value for "topography_grid_convention" is ',               &
                'not recognized. & Choose ''cell_edge'' or ''cell_center''.'
           CALL message( 'user_check_parameters', 'PA0240', 1, 2, 0, 6, 0 )
@@ -738,5 +740,5 @@
        ENDIF
 
-    ELSEIF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.  &
+    ELSEIF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.                       &
              TRIM( coupling_char ) == '_O' )  THEN
 
@@ -744,6 +746,5 @@
 !--    Check whether an (uncoupled) atmospheric run has been declared as an
 !--    ocean run (this setting is done via mrun-option -y)
-
-       message_string = 'ocean = .F. does not allow coupling_char = "' // &
+       message_string = 'ocean = .F. does not allow coupling_char = "' //      &
                         TRIM( coupling_char ) // '" set by mrun-option "-y"'
        CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
@@ -777,5 +778,5 @@
           gamma_mg = 1
        ELSE
-          message_string = 'unknown multigrid cycle: cycle_mg = "' // &
+          message_string = 'unknown multigrid cycle: cycle_mg = "' //          &
                            TRIM( cycle_mg ) // '"'
           CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
@@ -783,16 +784,16 @@
     ENDIF
 
-    IF ( fft_method /= 'singleton-algorithm'  .AND.  &
-         fft_method /= 'temperton-algorithm'  .AND.  &
-         fft_method /= 'fftw'                 .AND.  &
+    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 = "' // &
+       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. & 
+    IF( momentum_advec == 'ws-scheme' .AND.                                    & 
         .NOT. call_psolver_at_all_substeps  ) THEN
-        message_string = 'psolver must be called at each RK3 substep when "'//&
+        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 )
@@ -802,5 +803,5 @@
     IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' ) &
     THEN
-       message_string = 'unknown advection scheme: momentum_advec = "' // &
+       message_string = 'unknown advection scheme: momentum_advec = "' //      &
                         TRIM( momentum_advec ) // '"'
        CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
@@ -811,5 +812,5 @@
                    timestep_scheme == 'runge-kutta-2' ) )                      &
     THEN
-       message_string = 'momentum_advec or scalar_advec = "' &
+       message_string = 'momentum_advec or scalar_advec = "'                   &
          // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // &
          TRIM( timestep_scheme ) // '"'
@@ -819,5 +820,5 @@
          scalar_advec /= 'ws-scheme-mono' .AND. scalar_advec /= 'bc-scheme' )  &
     THEN
-       message_string = 'unknown advection scheme: scalar_advec = "' // &
+       message_string = 'unknown advection scheme: scalar_advec = "' //        &
                         TRIM( scalar_advec ) // '"'
        CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
@@ -825,5 +826,5 @@
     IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )    &
     THEN
-       message_string = 'advection_scheme scalar_advec = "' &
+       message_string = 'advection_scheme scalar_advec = "'                    &
          // TRIM( scalar_advec ) // '" not implemented for & loop_optimization = "' // &
          TRIM( loop_optimization ) // '"'
@@ -851,5 +852,5 @@
 !-- Set LOGICAL switches to enhance performance
     IF ( momentum_advec == 'ws-scheme' )       ws_scheme_mom = .TRUE.
-    IF ( scalar_advec   == 'ws-scheme' .OR.                                   &
+    IF ( scalar_advec   == 'ws-scheme' .OR.                                    &
          scalar_advec   == 'ws-scheme-mono' )  ws_scheme_sca = .TRUE.
     IF ( scalar_advec   == 'ws-scheme-mono' )  monotonic_adjustment = .TRUE.
@@ -908,19 +909,19 @@
     IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
 
-    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
          TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
 !
 !--    No restart run: several initialising actions are possible
        action = initializing_actions
-       DO WHILE ( TRIM( action ) /= '' )
+       DO  WHILE ( TRIM( action ) /= '' )
           position = INDEX( action, ' ' )
           SELECT CASE ( action(1:position-1) )
 
-             CASE ( 'set_constant_profiles', 'set_1d-model_profiles', &
+             CASE ( 'set_constant_profiles', 'set_1d-model_profiles',          &
                     'by_user', 'initialize_vortex',     'initialize_ptanom' )
                 action = action(position+1:)
 
              CASE DEFAULT
-                message_string = 'initializing_action = "' // &
+                message_string = 'initializing_action = "' //                  &
                                  TRIM( action ) // '" unkown or not allowed'
                 CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
@@ -930,7 +931,7 @@
     ENDIF
 
-    IF ( TRIM( initializing_actions ) == 'initialize_vortex' .AND. &
+    IF ( TRIM( initializing_actions ) == 'initialize_vortex'  .AND.            &
          conserve_volume_flow ) THEN
-         message_string = 'initializing_actions = "initialize_vortex"' // &
+         message_string = 'initializing_actions = "initialize_vortex"' //      &
                         ' ist not allowed with conserve_volume_flow = .T.'
        CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
@@ -938,28 +939,28 @@
 
 
-    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 ' //  &
+       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
 
-    IF ( cloud_physics  .AND.  .NOT. humidity )  THEN
-       WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ', &
+    IF ( cloud_physics  .AND.  .NOT.  humidity )  THEN
+       WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ',   &
               'not allowed with humidity = ', humidity
        CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
@@ -967,5 +968,5 @@
 
     IF ( precipitation  .AND.  .NOT.  cloud_physics )  THEN
-       WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ', &
+       WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ',   &
               'not allowed with cloud_physics = ', cloud_physics
        CALL message( 'check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
@@ -973,5 +974,5 @@
 
     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 )
@@ -979,5 +980,5 @@
 
     IF ( passive_scalar  .AND.  humidity )  THEN
-       message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' // &
+       message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' //    &
                         'is not allowed simultaneously'
        CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
@@ -1032,5 +1033,5 @@
 !--    calculated from the temperature/humidity gradients in the land surface
 !--    model
-       IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' )  THEN
+       IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
           message_string = 'lsm requires setting of'//                         &
                            'bc_pt_b = "dirichlet" and '//                      &
@@ -1039,5 +1040,5 @@
        ENDIF
 
-       IF ( .NOT. constant_flux_layer )  THEN
+       IF (  .NOT.  constant_flux_layer )  THEN
           message_string = 'lsm requires '//                                   &
                            'constant_flux_layer = .T.'
@@ -1046,11 +1047,19 @@
 
        IF ( topography /= 'flat' )  THEN
-          message_string = 'lsm cannot be used ' //  & 
+          message_string = 'lsm cannot be used ' //                            & 
                            'in combination with  topography /= "flat"'
           CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
        ENDIF
 
+       IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                       &
+              most_method == 'lookup' )  THEN
+           WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ',     &
+                                      'allowed in combination with ',          &
+                                      'most_method = ', most_method
+          CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
+       ENDIF
+
        IF ( veg_type == 0 )  THEN
-          IF ( SUM(root_fraction) /= 1.0_wp)  THEN
+          IF ( SUM( root_fraction ) /= 1.0_wp )  THEN
              message_string = 'veg_type = 0 (user_defined)'//                  &
                               'requires setting of root_fraction(0:3)'//       &
@@ -1059,5 +1068,5 @@
           ENDIF
  
-          IF ( min_canopy_resistance == 9999999.9_wp)  THEN
+          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
              message_string = 'veg_type = 0 (user defined)'//                  &
                               'requires setting of min_canopy_resistance'//    &
@@ -1066,5 +1075,5 @@
           ENDIF
 
-          IF ( leaf_area_index == 9999999.9_wp)  THEN
+          IF ( leaf_area_index == 9999999.9_wp )  THEN
              message_string = 'veg_type = 0 (user_defined)'//                  &
                               'requires setting of leaf_area_index'//          &
@@ -1073,5 +1082,5 @@
           ENDIF
 
-          IF ( vegetation_coverage == 9999999.9_wp)  THEN
+          IF ( vegetation_coverage == 9999999.9_wp )  THEN
              message_string = 'veg_type = 0 (user_defined)'//                  &
                               'requires setting of vegetation_coverage'//      &
@@ -1087,5 +1096,5 @@
           ENDIF
 
-          IF ( lambda_surface_stable == 9999999.9_wp)  THEN
+          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
              message_string = 'veg_type = 0 (user_defined)'//                  &
                               'requires setting of lambda_surface_stable'//    &
@@ -1094,5 +1103,5 @@
           ENDIF
 
-          IF ( lambda_surface_unstable == 9999999.9_wp)  THEN
+          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
              message_string = 'veg_type = 0 (user_defined)'//                  &
                               'requires setting of lambda_surface_unstable'//  &
@@ -1101,5 +1110,5 @@
           ENDIF
 
-          IF ( f_shortwave_incoming == 9999999.9_wp)  THEN
+          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
              message_string = 'veg_type = 0 (user_defined)'//                  &
                               'requires setting of f_shortwave_incoming'//     &
@@ -1108,5 +1117,5 @@
           ENDIF
 
-          IF ( z0_eb == 9999999.9_wp)  THEN
+          IF ( z0_eb == 9999999.9_wp )  THEN
              message_string = 'veg_type = 0 (user_defined)'//                  &
                               'requires setting of z0_eb'//                   &
@@ -1115,5 +1124,5 @@
           ENDIF
 
-          IF ( z0h_eb == 9999999.9_wp)  THEN
+          IF ( z0h_eb == 9999999.9_wp )  THEN
              message_string = 'veg_type = 0 (user_defined)'//                  &
                               'requires setting of z0h_eb'//                  &
@@ -1127,5 +1136,5 @@
        IF ( soil_type == 0 )  THEN
 
-          IF ( alpha_vangenuchten == 9999999.9_wp)  THEN
+          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
              message_string = 'soil_type = 0 (user_defined)'//                 &
                               'requires setting of alpha_vangenuchten'//       &
@@ -1134,5 +1143,5 @@
           ENDIF
 
-          IF ( l_vangenuchten == 9999999.9_wp)  THEN
+          IF ( l_vangenuchten == 9999999.9_wp )  THEN
              message_string = 'soil_type = 0 (user_defined)'//                 &
                               'requires setting of l_vangenuchten'//           &
@@ -1141,5 +1150,5 @@
           ENDIF
 
-          IF ( n_vangenuchten == 9999999.9_wp)  THEN
+          IF ( n_vangenuchten == 9999999.9_wp )  THEN
              message_string = 'soil_type = 0 (user_defined)'//                 &
                               'requires setting of n_vangenuchten'//           &
@@ -1148,5 +1157,5 @@
           ENDIF
 
-          IF ( hydraulic_conductivity == 9999999.9_wp)  THEN
+          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
              message_string = 'soil_type = 0 (user_defined)'//                 &
                               'requires setting of hydraulic_conductivity'//   &
@@ -1155,5 +1164,5 @@
           ENDIF
 
-          IF ( saturation_moisture == 9999999.9_wp)  THEN
+          IF ( saturation_moisture == 9999999.9_wp )  THEN
              message_string = 'soil_type = 0 (user_defined)'//                 &
                               'requires setting of saturation_moisture'//      &
@@ -1162,5 +1171,5 @@
           ENDIF
 
-          IF ( field_capacity == 9999999.9_wp)  THEN
+          IF ( field_capacity == 9999999.9_wp )  THEN
              message_string = 'soil_type = 0 (user_defined)'//                 &
                               'requires setting of field_capacity'//           &
@@ -1169,5 +1178,5 @@
           ENDIF
 
-          IF ( wilting_point == 9999999.9_wp)  THEN
+          IF ( wilting_point == 9999999.9_wp )  THEN
              message_string = 'soil_type = 0 (user_defined)'//                 &
                               'requires setting of wilting_point'//            &
@@ -1176,5 +1185,5 @@
           ENDIF
 
-          IF ( residual_moisture == 9999999.9_wp)  THEN
+          IF ( residual_moisture == 9999999.9_wp )  THEN
              message_string = 'soil_type = 0 (user_defined)'//                 &
                               'requires setting of residual_moisture'//        &
@@ -1185,5 +1194,5 @@
        ENDIF
 
-       IF ( .NOT. radiation )  THEN
+       IF (  .NOT.  radiation )  THEN
           message_string = 'lsm requires '//                                   &
                            'radiation = .T.'
@@ -1194,6 +1203,6 @@
 
     IF ( radiation )  THEN
-       IF ( radiation_scheme /= 'constant'  .AND.                              &
-            radiation_scheme /= 'clear-sky' .AND.                              &
+       IF ( radiation_scheme /= 'constant'   .AND.                             &
+            radiation_scheme /= 'clear-sky'  .AND.                             &
             radiation_scheme /= 'rrtmg' )  THEN
           message_string = 'unknown radiation_scheme = '//                     &
@@ -1215,5 +1224,5 @@
 
        ENDIF
-       IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND.                &
+       IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.             &
             radiation_scheme == 'clear-sky')  THEN
           message_string = 'radiation_scheme = "clear-sky" in combination' //  & 
@@ -1222,5 +1231,5 @@
           CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
        ENDIF
-       IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND.           &
+       IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.        &
           (    albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
           .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
@@ -1235,5 +1244,5 @@
        ENDIF
        IF ( topography /= 'flat' )  THEN
-          message_string = 'radiation scheme cannot be used ' //  & 
+          message_string = 'radiation scheme cannot be used ' //               & 
                            'in combination with  topography /= "flat"'
           CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
@@ -1244,5 +1253,5 @@
                  loop_optimization == 'vector' )                               &
          .AND.  cloud_physics  .AND.  icloud_scheme == 0 )  THEN
-       message_string = 'cloud_scheme = seifert_beheng requires ' // &
+       message_string = 'cloud_scheme = seifert_beheng requires ' //           &
                         'loop_optimization = "cache" or "vector"'
        CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
@@ -1270,11 +1279,11 @@
        gradient = 0.0_wp
 
-       IF ( .NOT. ocean )  THEN
+       IF (  .NOT.  ocean )  THEN
 
           ug_vertical_gradient_level_ind(1) = 0
           ug(0) = ug_surface
           DO  k = 1, nzt+1
-             IF ( i < 11 ) THEN
-                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
+             IF ( i < 11 )  THEN
+                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
@@ -1299,6 +1308,6 @@
           ug(nzt+1) = ug_surface
           DO  k = nzt, nzb, -1
-             IF ( i < 11 ) THEN
-                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND. &
+             IF ( i < 11 )  THEN
+                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
@@ -1334,11 +1343,11 @@
        gradient = 0.0_wp
 
-       IF ( .NOT. ocean )  THEN
+       IF (  .NOT.  ocean )  THEN
 
           vg_vertical_gradient_level_ind(1) = 0
           vg(0) = vg_surface
           DO  k = 1, nzt+1
-             IF ( i < 11 ) THEN 
-                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
+             IF ( i < 11 )  THEN 
+                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
@@ -1363,6 +1372,6 @@
           vg(nzt+1) = vg_surface
           DO  k = nzt, nzb, -1
-             IF ( i < 11 ) THEN
-                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND. &
+             IF ( i < 11 )  THEN
+                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
@@ -1415,5 +1424,5 @@
 
              IF ( kk < 100 )  THEN
-                DO WHILE ( uv_heights(kk+1) <= zu(k) )
+                DO  WHILE ( uv_heights(kk+1) <= zu(k) )
                    kk = kk + 1
                    IF ( kk == 100 )  EXIT
@@ -1421,5 +1430,5 @@
              ENDIF
 
-             IF ( kk < 100 .AND. uv_heights(kk+1) /= 9999999.9_wp )  THEN
+             IF ( kk < 100  .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) ) * &
@@ -1444,15 +1453,15 @@
 !
 !--    Compute initial temperature profile using the given temperature gradients
-       IF ( .NOT. neutral )  THEN
+       IF (  .NOT.  neutral )  THEN
 
           i = 1
           gradient = 0.0_wp
 
-          IF ( .NOT. ocean )  THEN
+          IF (  .NOT.  ocean )  THEN
 
              pt_vertical_gradient_level_ind(1) = 0
              DO  k = 1, nzt+1
-                IF ( i < 11 ) THEN
-                   IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
+                IF ( i < 11 )  THEN
+                   IF ( pt_vertical_gradient_level(i) < zu(k)  .AND.           &
                         pt_vertical_gradient_level(i) >= 0.0_wp )  THEN
                       gradient = pt_vertical_gradient(i) / 100.0_wp
@@ -1476,6 +1485,6 @@
              pt_vertical_gradient_level_ind(1) = nzt+1
              DO  k = nzt, 0, -1
-                IF ( i < 11 ) THEN
-                   IF ( pt_vertical_gradient_level(i) > zu(k)  .AND. &
+                IF ( i < 11 )  THEN
+                   IF ( pt_vertical_gradient_level(i) > zu(k)  .AND.           &
                         pt_vertical_gradient_level(i) <= 0.0_wp )  THEN
                       gradient = pt_vertical_gradient(i) / 100.0_wp
@@ -1533,6 +1542,6 @@
           q_vertical_gradient_level_ind(1) = 0
           DO  k = 1, nzt+1
-             IF ( i < 11 ) THEN
-                IF ( q_vertical_gradient_level(i) < zu(k)  .AND. &
+             IF ( i < 11 )  THEN
+                IF ( q_vertical_gradient_level(i) < zu(k)  .AND.               &
                      q_vertical_gradient_level(i) >= 0.0_wp )  THEN
                    gradient = q_vertical_gradient(i) / 100.0_wp
@@ -1579,6 +1588,6 @@
           sa_vertical_gradient_level_ind(1) = nzt+1
           DO  k = nzt, 0, -1
-             IF ( i < 11 ) THEN
-                IF ( sa_vertical_gradient_level(i) > zu(k)  .AND. &
+             IF ( i < 11 )  THEN
+                IF ( sa_vertical_gradient_level(i) > zu(k)  .AND.              &
                      sa_vertical_gradient_level(i) <= 0.0_wp )  THEN
                    gradient = sa_vertical_gradient(i) / 100.0_wp
@@ -1606,10 +1615,10 @@
 !
 !-- 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 ' // &
+    IF ( use_subsidence_tendencies  .AND.  .NOT.  large_scale_subsidence )  THEN
+       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 ' // &
+       message_string = 'The usage of use_subsidence_tendencies ' //           &
                             'requires large_scale_forcing = .T..'
        CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
@@ -1619,6 +1628,6 @@
 !-- 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
@@ -1627,15 +1636,15 @@
 !--    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.'
+       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 )
@@ -1659,5 +1668,5 @@
        vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
     ELSE
-       message_string = 'illegal value for reference_state: "' // &
+       message_string = 'illegal value for reference_state: "' //              &
                         TRIM( reference_state ) // '"'
        CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
@@ -1686,5 +1695,5 @@
     IF ( alpha_surface /= 0.0_wp )  THEN
        IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
-          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface, &
+          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface,   &
                                      ' ) must be < 90.0'
           CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
@@ -1716,5 +1725,5 @@
           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 )
@@ -1740,26 +1749,26 @@
 !-- 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.  &
+       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.  &
+       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' // &
+          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.  &
+       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' // &
+          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 ' // &
+          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 )
@@ -1782,27 +1791,27 @@
     IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
        IF ( psolver(1:9) /= 'multigrid' )  THEN
-          message_string = 'non-cyclic lateral boundaries do not allow ' // &
+          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' .AND.                               &
-            momentum_advec /= 'ws-scheme-mono' )                              &
+       IF ( momentum_advec /= 'pw-scheme'  .AND.                               &
+          ( momentum_advec /= 'ws-scheme'  .AND.                               &
+            momentum_advec /= 'ws-scheme-mono' )                               &
           )  THEN
 
-          message_string = 'non-cyclic lateral boundaries do not allow ' // &
+          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' .AND.                                  &
+       IF ( scalar_advec /= 'pw-scheme'  .AND.                                 &
+          ( scalar_advec /= 'ws-scheme'  .AND.                                 &
             scalar_advec /= 'ws-scheme-mono' )                                 &
           )  THEN
-          message_string = 'non-cyclic lateral boundaries do not allow ' // &
+          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 ' // &
+          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
                            'galilei_transformation = .T.'
           CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
@@ -1819,10 +1828,10 @@
           bc_e_b = 'neumann'
           ibc_e_b = 1
-          message_string = 'boundary condition bc_e_b changed to "' // &
+          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 = "' // &
+       message_string = 'unknown boundary condition: bc_e_b = "' //            &
                         TRIM( bc_e_b ) // '"'
        CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
@@ -1836,5 +1845,5 @@
        ibc_p_b = 1
     ELSE
-       message_string = 'unknown boundary condition: 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 )
@@ -1847,5 +1856,5 @@
        ibc_p_t = 1
     ELSE
-       message_string = 'unknown boundary condition: 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 )
@@ -1862,5 +1871,5 @@
           ibc_pt_b = 1
        ELSE
-          message_string = 'unknown boundary condition: 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 )
@@ -1877,5 +1886,5 @@
        ibc_pt_t = 3
     ELSE
-       message_string = 'unknown boundary condition: 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 )
@@ -1889,6 +1898,6 @@
           ELSEIF ( ibc_pt_b == 1 )  THEN
              constant_heatflux = .TRUE.
-             IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND.    &
-                  .NOT. land_surface )  THEN
+             IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.   &
+                  .NOT.  land_surface )  THEN
                 surface_heatflux = shf_surf(1)
              ELSE
@@ -1899,6 +1908,6 @@
     ELSE
         constant_heatflux = .TRUE.
-        IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND.         &
-             large_scale_forcing .AND. .NOT. land_surface )  THEN
+        IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.        &
+             large_scale_forcing  .AND.  .NOT.  land_surface )  THEN
            surface_heatflux = shf_surf(1)
         ENDIF
@@ -1909,5 +1918,11 @@
     IF ( neutral )  THEN
 
-       IF ( surface_heatflux /= 0.0_wp  .AND. surface_heatflux /= 9999999.9_wp ) &
+       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
           message_string = 'heatflux must not be set for pure neutral flow'
@@ -1915,18 +1930,12 @@
        ENDIF
 
-       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 )
-       ENDIF
-
-    ENDIF
-
-    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.  &
+    ENDIF
+
+    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 ' // &
+       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' //  &
                         'must be set'
        CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
@@ -1937,5 +1946,5 @@
 !-- temperature. In this case specification of a constant heat flux is
 !-- forbidden.
-    IF ( ibc_pt_b == 0  .AND.   constant_heatflux  .AND. &
+    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 ) //&
@@ -1944,6 +1953,6 @@
     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) = ', &
+       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 )
@@ -1954,5 +1963,5 @@
 !-- temperature. In this case specification of a constant heat flux is
 !-- forbidden.
-    IF ( ibc_pt_t == 0  .AND.   constant_top_heatflux  .AND. &
+    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 ) //&
@@ -1969,5 +1978,5 @@
           ibc_sa_t = 1
        ELSE
-          message_string = 'unknown boundary condition: bc_sa_t = "' // &
+          message_string = 'unknown boundary condition: bc_sa_t = "' //        &
                            TRIM( bc_sa_t ) // '"'
           CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
@@ -1975,7 +1984,7 @@
 
        IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
-       IF ( ibc_sa_t == 1  .AND.   top_salinityflux == 9999999.9_wp )  THEN
-          message_string = 'boundary condition: bc_sa_t = "' // &
-                           TRIM( bc_sa_t ) // '" requires to set ' // &
+       IF ( ibc_sa_t == 1  .AND.  top_salinityflux == 9999999.9_wp )  THEN
+          message_string = 'boundary condition: bc_sa_t = "' //                &
+                           TRIM( bc_sa_t ) // '" requires to set ' //          &
                            'top_salinityflux'
           CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
@@ -1986,8 +1995,8 @@
 !--    salinity. In this case specification of a constant salinity flux is
 !--    forbidden.
-       IF ( ibc_sa_t == 0  .AND.   constant_top_salinityflux  .AND. &
+       IF ( ibc_sa_t == 0  .AND.  constant_top_salinityflux  .AND.             &
             top_salinityflux /= 0.0_wp )  THEN
-          message_string = 'boundary condition: bc_sa_t = "' // &
-                           TRIM( bc_sa_t ) // '" is not allowed with ' // &
+          message_string = 'boundary condition: bc_sa_t = "' //                &
+                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
                            'constant_top_salinityflux = .TRUE.'
           CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
@@ -2010,5 +2019,5 @@
           ibc_q_b = 1
        ELSE
-          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
+          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
                            '_b ="' // TRIM( bc_q_b ) // '"'
           CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 )
@@ -2021,5 +2030,5 @@
           ibc_q_t = 3
        ELSE
-          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
+          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
                            '_t ="' // TRIM( bc_q_t ) // '"'
           CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 )
@@ -2040,6 +2049,6 @@
        ELSE
           constant_waterflux = .TRUE.
-          IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
-                 large_scale_forcing ) THEN
+          IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.      &
+               large_scale_forcing )  THEN
              surface_waterflux = qsws_surf(1)
           ENDIF
@@ -2058,5 +2067,5 @@
        IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0_wp )  THEN
           WRITE( message_string, * )  'a prescribed surface flux is not allo', &
-                 'wed with ', sq, '_surface_initial_change (/=0) = ', &
+                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
                  q_surface_initial_change
           CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
@@ -2071,5 +2080,5 @@
        ibc_uv_b = 1
        IF ( constant_flux_layer )  THEN
-          message_string = 'boundary condition: bc_uv_b = "' // &
+          message_string = 'boundary condition: bc_uv_b = "' //                &
                TRIM( bc_uv_b ) // '" is not allowed with constant_flux_layer'  &
                // ' = .TRUE.'
@@ -2077,5 +2086,5 @@
        ENDIF
     ELSE
-       message_string = 'unknown boundary condition: 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 )
@@ -2106,5 +2115,5 @@
           ibc_uv_t = 3
        ELSE
-          message_string = 'unknown boundary condition: 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 )
@@ -2117,7 +2126,7 @@
        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 = ', &
+       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 )
@@ -2126,5 +2135,5 @@
 
     IF ( rayleigh_damping_height == -1.0_wp )  THEN
-       IF ( .NOT. ocean )  THEN
+       IF (  .NOT.  ocean )  THEN
           rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
        ELSE
@@ -2132,15 +2141,15 @@
        ENDIF
     ELSE
-       IF ( .NOT. ocean )  THEN
-          IF ( rayleigh_damping_height < 0.0_wp  .OR. &
+       IF (  .NOT.  ocean )  THEN
+          IF ( rayleigh_damping_height < 0.0_wp  .OR.                          &
                rayleigh_damping_height > zu(nzt) )  THEN
-             WRITE( message_string, * )  'rayleigh_damping_height = ', &
+             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. &
+          IF ( rayleigh_damping_height > 0.0_wp  .OR.                          &
                rayleigh_damping_height < zu(nzb) )  THEN
-             WRITE( message_string, * )  'rayleigh_damping_height = ', &
+             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 )
@@ -2154,11 +2163,11 @@
 !-- be opened (cf. check_open)
     IF ( statistic_regions > 9  .OR.  statistic_regions < 0 )  THEN
-       WRITE ( message_string, * ) 'number of statistic_regions = ', &
+       WRITE ( message_string, * ) 'number of statistic_regions = ',           &
                    statistic_regions+1, ' but only 10 regions are allowed'
        CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( normalizing_region > statistic_regions  .OR. &
+    IF ( normalizing_region > statistic_regions  .OR.                          &
          normalizing_region < 0)  THEN
-       WRITE ( message_string, * ) 'normalizing_region = ', &
+       WRITE ( message_string, * ) 'normalizing_region = ',                    &
                 normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
                 ' (value of statistic_regions)'
@@ -2184,20 +2193,20 @@
 !
 !-- Set the default skip time intervals for data output, if necessary
-    IF ( skip_time_dopr    == 9999999.9_wp ) &
+    IF ( skip_time_dopr    == 9999999.9_wp )                                   &
                                        skip_time_dopr    = skip_time_data_output
-    IF ( skip_time_dosp    == 9999999.9_wp ) &
+    IF ( skip_time_dosp    == 9999999.9_wp )                                   &
                                        skip_time_dosp    = skip_time_data_output
-    IF ( skip_time_do2d_xy == 9999999.9_wp ) &
+    IF ( skip_time_do2d_xy == 9999999.9_wp )                                   &
                                        skip_time_do2d_xy = skip_time_data_output
-    IF ( skip_time_do2d_xz == 9999999.9_wp ) &
+    IF ( skip_time_do2d_xz == 9999999.9_wp )                                   &
                                        skip_time_do2d_xz = skip_time_data_output
-    IF ( skip_time_do2d_yz == 9999999.9_wp ) &
+    IF ( skip_time_do2d_yz == 9999999.9_wp )                                   &
                                        skip_time_do2d_yz = skip_time_data_output
-    IF ( skip_time_do3d    == 9999999.9_wp ) &
+    IF ( skip_time_do3d    == 9999999.9_wp )                                   &
                                        skip_time_do3d    = skip_time_data_output
-    IF ( skip_time_data_output_av == 9999999.9_wp ) &
+    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 ) &
+       IF ( skip_time_domask(mid) == 9999999.9_wp )                            &
                                 skip_time_domask(mid)    = skip_time_data_output
     ENDDO
@@ -2207,5 +2216,5 @@
 !-- spectra)
     IF ( averaging_interval > dt_data_output_av )  THEN
-       WRITE( message_string, * )  'averaging_interval = ', &
+       WRITE( message_string, * )  'averaging_interval = ',                    &
              averaging_interval, ' must be <= dt_data_output = ', dt_data_output
        CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
@@ -2217,5 +2226,5 @@
 
     IF ( averaging_interval_pr > dt_dopr )  THEN
-       WRITE( message_string, * )  'averaging_interval_pr = ', &
+       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 )
@@ -2227,5 +2236,5 @@
 
     IF ( averaging_interval_sp > dt_dosp )  THEN
-       WRITE( message_string, * )  'averaging_interval_sp = ', &
+       WRITE( message_string, * )  'averaging_interval_sp = ',                 &
              averaging_interval_sp, ' must be <= dt_dosp = ', dt_dosp
        CALL message( 'check_parameters', 'PA0087', 1, 2, 0, 6, 0 )
@@ -2252,6 +2261,6 @@
 !-- 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 = ', &
+       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 )
@@ -2259,5 +2268,5 @@
 
     IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
-       WRITE( message_string, * )  'dt_averaging_input_pr = ', &
+       WRITE( message_string, * )  'dt_averaging_input_pr = ',                 &
                 dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
                 averaging_interval_pr
@@ -2272,6 +2281,6 @@
        ELSE
           IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
-             WRITE( message_string, * )  'precipitation_amount_interval = ', &
-                 precipitation_amount_interval, ' must not be larger than ', &
+             WRITE( message_string, * )  'precipitation_amount_interval = ',   &
+                 precipitation_amount_interval, ' must not be larger than ',   &
                  'dt_do2d_xy = ', dt_do2d_xy
              CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
@@ -2532,5 +2541,5 @@
           CASE ( 's', '#s' )
              IF ( .NOT. passive_scalar )  THEN
-                message_string = 'data_output_pr = ' // &
+                message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for passive_scalar = .FALSE.'
@@ -2573,5 +2582,5 @@
           CASE ( 'lpt', '#lpt' )
              IF ( .NOT. cloud_physics ) THEN
-                message_string = 'data_output_pr = ' // &
+                message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for cloud_physics = .FALSE.'
@@ -2616,6 +2625,6 @@
 
           CASE ( 'w"q"' )
-             IF ( .NOT. humidity )  THEN
-                message_string = 'data_output_pr = ' // &
+             IF (  .NOT.  humidity )  THEN
+                message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for humidity = .FALSE.'
@@ -2628,6 +2637,6 @@
 
           CASE ( 'w*q*' )
-             IF ( .NOT. humidity )  THEN
-                message_string = 'data_output_pr = ' // &
+             IF (  .NOT.  humidity )  THEN
+                message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for humidity = .FALSE.'
@@ -2640,6 +2649,6 @@
 
           CASE ( 'wq' )
-             IF ( .NOT. humidity )  THEN
-                message_string = 'data_output_pr = ' // &
+             IF (  .NOT.  humidity )  THEN
+                message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for humidity = .FALSE.'
@@ -2652,6 +2661,6 @@
 
           CASE ( 'w"s"' )
-             IF ( .NOT. passive_scalar ) THEN
-                message_string = 'data_output_pr = ' // &
+             IF (  .NOT.  passive_scalar )  THEN
+                message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for passive_scalar = .FALSE.'
@@ -2664,6 +2673,6 @@
 
           CASE ( 'w*s*' )
-             IF ( .NOT. passive_scalar ) THEN
-                message_string = 'data_output_pr = ' // &
+             IF (  .NOT.  passive_scalar )  THEN
+                message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for passive_scalar = .FALSE.'
@@ -2676,6 +2685,6 @@
 
           CASE ( 'ws' )
-             IF ( .NOT. passive_scalar ) THEN
-                message_string = 'data_output_pr = ' // &
+             IF (  .NOT.  passive_scalar )  THEN
+                message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for passive_scalar = .FALSE.'
@@ -2688,15 +2697,14 @@
 
           CASE ( 'w"qv"' )
-             IF ( humidity  .AND.  .NOT. cloud_physics ) &
-             THEN
+             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
                 dopr_index(i) = 48
                 dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
-             ELSEIF( humidity .AND. cloud_physics ) THEN
+             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
                 dopr_index(i) = 51
                 dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
-                message_string = 'data_output_pr = ' // &
+                message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for cloud_physics = .FALSE. an&' // &
@@ -2706,5 +2714,5 @@
 
           CASE ( 'w*qv*' )
-             IF ( humidity  .AND.  .NOT. cloud_physics ) &
+             IF ( humidity  .AND.  .NOT. cloud_physics )                       &
              THEN
                 dopr_index(i) = 49
@@ -2716,5 +2724,5 @@
                 hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
-                message_string = 'data_output_pr = ' // &
+                message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
                                  'lemented for cloud_physics = .FALSE. an&' // &
@@ -2724,10 +2732,9 @@
 
           CASE ( 'wqv' )
-             IF ( humidity  .AND.  .NOT. cloud_physics ) &
-             THEN
+             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
                 dopr_index(i) = 50
                 dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
-             ELSEIF( humidity .AND. cloud_physics ) THEN
+             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
                 dopr_index(i) = 53
                 dopr_unit(i)  = 'kg/kg m/s'
@@ -2742,5 +2749,5 @@
 
           CASE ( 'ql' )
-             IF ( .NOT. cloud_physics  .AND.  .NOT. cloud_droplets )  THEN
+             IF (  .NOT.  cloud_physics  .AND.  .NOT.  cloud_droplets )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2801,5 +2808,5 @@
 
           CASE ( 'rho' )
-             IF ( .NOT. ocean ) THEN
+             IF (  .NOT.  ocean ) THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2819,5 +2826,5 @@
 
           CASE ( 'w"sa"' )
-             IF ( .NOT. ocean ) THEN
+             IF (  .NOT.  ocean ) THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2831,5 +2838,5 @@
 
           CASE ( 'w*sa*' )
-             IF ( .NOT. ocean ) THEN
+             IF (  .NOT. ocean  ) THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2843,5 +2850,5 @@
 
           CASE ( 'wsa' )
-             IF ( .NOT. ocean ) THEN
+             IF (  .NOT.  ocean ) THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2865,5 +2872,5 @@
 
           CASE ( 'q*2' )
-             IF ( .NOT. humidity )  THEN
+             IF (  .NOT.  humidity )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2877,5 +2884,5 @@
 
           CASE ( 'prho' )
-             IF ( .NOT. ocean ) THEN
+             IF (  .NOT.  ocean ) THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2894,5 +2901,5 @@
 
           CASE ( 'nr' )
-             IF ( .NOT. cloud_physics )  THEN
+             IF (  .NOT.  cloud_physics )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2904,5 +2911,5 @@
                                  'lemented for cloud_scheme /= seifert_beheng'
                 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
-             ELSEIF ( .NOT. precipitation )  THEN
+             ELSEIF (  .NOT.  precipitation )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2916,5 +2923,5 @@
 
           CASE ( 'qr' )
-             IF ( .NOT. cloud_physics )  THEN
+             IF (  .NOT.  cloud_physics )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2926,5 +2933,5 @@
                                  'lemented for cloud_scheme /= seifert_beheng'
                 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
-             ELSEIF ( .NOT. precipitation )  THEN
+             ELSEIF (  .NOT.  precipitation )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2938,5 +2945,5 @@
 
           CASE ( 'qc' )
-             IF ( .NOT. cloud_physics )  THEN
+             IF (  .NOT.  cloud_physics )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2955,5 +2962,5 @@
 
           CASE ( 'prr' )
-             IF ( .NOT. cloud_physics )  THEN
+             IF (  .NOT.  cloud_physics )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2965,5 +2972,5 @@
                                  'lemented for cloud_scheme /= seifert_beheng'
                 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
-             ELSEIF ( .NOT. precipitation )  THEN
+             ELSEIF (  .NOT.  precipitation )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -2988,5 +2995,5 @@
 
           CASE ( 'w_subs' )
-             IF ( .NOT. large_scale_subsidence )  THEN
+             IF (  .NOT.  large_scale_subsidence )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3000,5 +3007,5 @@
 
           CASE ( 'td_lsa_lpt' )
-             IF ( .NOT. large_scale_forcing )  THEN
+             IF (  .NOT.  large_scale_forcing )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3012,5 +3019,5 @@
 
           CASE ( 'td_lsa_q' )
-             IF ( .NOT. large_scale_forcing )  THEN
+             IF (  .NOT.  large_scale_forcing )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3024,5 +3031,5 @@
 
           CASE ( 'td_sub_lpt' )
-             IF ( .NOT. large_scale_forcing )  THEN
+             IF (  .NOT.  large_scale_forcing )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3036,5 +3043,5 @@
 
           CASE ( 'td_sub_q' )
-             IF ( .NOT. large_scale_forcing )  THEN
+             IF (  .NOT.  large_scale_forcing )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3048,5 +3055,5 @@
 
           CASE ( 'td_nud_lpt' )
-             IF ( .NOT. nudging )  THEN
+             IF (  .NOT.  nudging )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3060,5 +3067,5 @@
 
           CASE ( 'td_nud_q' )
-             IF ( .NOT. nudging )  THEN
+             IF (  .NOT.  nudging )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3072,5 +3079,5 @@
 
           CASE ( 'td_nud_u' )
-             IF ( .NOT. nudging )  THEN
+             IF (  .NOT.  nudging )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3084,5 +3091,5 @@
 
           CASE ( 'td_nud_v' )
-             IF ( .NOT. nudging )  THEN
+             IF (  .NOT.  nudging )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3096,5 +3103,5 @@
 
           CASE ( 't_soil', '#t_soil' )
-             IF ( .NOT. land_surface )  THEN
+             IF (  .NOT.  land_surface )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3113,5 +3120,5 @@
 
           CASE ( 'm_soil', '#m_soil' )
-             IF ( .NOT. land_surface )  THEN
+             IF (  .NOT.  land_surface )  THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not imp' // &
@@ -3130,5 +3137,6 @@
 
           CASE ( 'rad_net' )
-             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
+             THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not ava' // &
@@ -3143,5 +3151,6 @@
 
           CASE ( 'rad_lw_in' )
-             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
+             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
+             THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not ava' // &
@@ -3156,5 +3165,6 @@
 
           CASE ( 'rad_lw_out' )
-             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
+             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
+             THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not ava' // &
@@ -3169,5 +3179,6 @@
 
           CASE ( 'rad_sw_in' )
-             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
+             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
+             THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not ava' // &
@@ -3182,5 +3193,6 @@
 
           CASE ( 'rad_sw_out')
-             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' ) &
+             THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not ava' // &
@@ -3195,5 +3207,6 @@
 
           CASE ( 'rad_lw_cs_hr' )
-             IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not ava' // &
@@ -3208,5 +3221,6 @@
 
           CASE ( 'rad_lw_hr' )
-             IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not ava' // &
@@ -3221,5 +3235,6 @@
 
           CASE ( 'rad_sw_cs_hr' )
-             IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not ava' // &
@@ -3234,5 +3249,6 @@
 
           CASE ( 'rad_sw_hr' )
-             IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
                 message_string = 'data_output_pr = ' //                        &
                                  TRIM( data_output_pr(i) ) // ' is not ava' // &
@@ -3329,5 +3345,5 @@
 
           CASE ( 'lpt' )
-             IF ( .NOT. cloud_physics )  THEN
+             IF (  .NOT.  cloud_physics )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'res cloud_physics = .TRUE.'
@@ -3337,5 +3353,5 @@
 
           CASE ( 'm_soil' )
-             IF ( .NOT. land_surface )  THEN
+             IF (  .NOT.  land_surface )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'land_surface = .TRUE.'
@@ -3345,5 +3361,5 @@
 
           CASE ( 'nr' )
-             IF ( .NOT. cloud_physics )  THEN
+             IF (  .NOT.  cloud_physics )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'res cloud_physics = .TRUE.'
@@ -3357,5 +3373,5 @@
 
           CASE ( 'pc', 'pr' )
-             IF ( .NOT. particle_advection )  THEN
+             IF (  .NOT.  particle_advection )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requir' // &
                    'es a "particles_par"-NAMELIST in the parameter file (PARIN)'
@@ -3366,5 +3382,5 @@
 
           CASE ( 'prr' )
-             IF ( .NOT. cloud_physics )  THEN
+             IF (  .NOT.  cloud_physics )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'res cloud_physics = .TRUE.'
@@ -3374,5 +3390,5 @@
                          'res cloud_scheme = seifert_beheng'
                 CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
-             ELSEIF ( .NOT. precipitation )  THEN
+             ELSEIF (  .NOT.  precipitation )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res precipitation = .TRUE.'
@@ -3382,5 +3398,5 @@
 
           CASE ( 'q', 'vpt' )
-             IF ( .NOT. humidity )  THEN
+             IF (  .NOT.  humidity )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res humidity = .TRUE.'
@@ -3391,5 +3407,5 @@
 
           CASE ( 'qc' )
-             IF ( .NOT. cloud_physics )  THEN
+             IF (  .NOT.  cloud_physics )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'res cloud_physics = .TRUE.'
@@ -3403,5 +3419,5 @@
 
           CASE ( 'ql' )
-             IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
+             IF ( .NOT.  ( cloud_physics  .OR.  cloud_droplets ) )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.'
@@ -3411,5 +3427,5 @@
 
           CASE ( 'ql_c', 'ql_v', 'ql_vp' )
-             IF ( .NOT. cloud_droplets )  THEN
+             IF (  .NOT.  cloud_droplets )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res cloud_droplets = .TRUE.'
@@ -3421,5 +3437,5 @@
 
           CASE ( 'qr' )
-             IF ( .NOT. cloud_physics )  THEN
+             IF (  .NOT.  cloud_physics )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'res cloud_physics = .TRUE.'
@@ -3429,5 +3445,5 @@
                          'res cloud_scheme = seifert_beheng'
                 CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
-             ELSEIF ( .NOT. precipitation )  THEN
+             ELSEIF (  .NOT.  precipitation )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res precipitation = .TRUE.'
@@ -3437,5 +3453,5 @@
 
           CASE ( 'qv' )
-             IF ( .NOT. cloud_physics )  THEN
+             IF (  .NOT.  cloud_physics )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res cloud_physics = .TRUE.'
@@ -3446,5 +3462,5 @@
           CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_cs_hr', 'rad_lw_hr',       &
                  'rad_sw_in', 'rad_sw_out', 'rad_sw_cs_hr', 'rad_sw_hr' )
-             IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' )  THEN
+             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
                 message_string = '"output of "' // TRIM( var ) // '" requi' // &
                                  'res radiation = .TRUE. and ' //              &
@@ -3455,5 +3471,5 @@
 
           CASE ( 'rho' )
-             IF ( .NOT. ocean )  THEN
+             IF (  .NOT.  ocean )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res ocean = .TRUE.'
@@ -3463,5 +3479,5 @@
 
           CASE ( 's' )
-             IF ( .NOT. passive_scalar )  THEN
+             IF (  .NOT.  passive_scalar )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res passive_scalar = .TRUE.'
@@ -3471,5 +3487,5 @@
 
           CASE ( 'sa' )
-             IF ( .NOT. ocean )  THEN
+             IF (  .NOT.  ocean )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res ocean = .TRUE.'
@@ -3479,5 +3495,5 @@
 
           CASE ( 't_soil' )
-             IF ( .NOT. land_surface )  THEN
+             IF (  .NOT.  land_surface )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                          'land_surface = .TRUE.'
@@ -3491,5 +3507,6 @@
                  'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', 'rad_net*',  &
                  'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*',   &
-                 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*' )
+                 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*', &
+                 'z0q*' )
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
                 message_string = 'illegal value for data_output: "' //         &
@@ -3498,8 +3515,8 @@
                 CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" )  THEN
-                IF ( TRIM( var ) == 'rrtm_aldif*' .OR.                         &
-                     TRIM( var ) == 'rrtm_aldir*' .OR.                         &
-                     TRIM( var ) == 'rrtm_asdif*' .OR.                         &
+             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
+                IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
+                     TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
+                     TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
                      TRIM( var ) == 'rrtm_asdir*'      )                       &
                 THEN
@@ -3511,10 +3528,10 @@
              ENDIF
 
-             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT. land_surface )  THEN
+             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res land_surface = .TRUE.'
                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT. land_surface )  THEN
+             IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res land_surface = .TRUE.'
@@ -3526,10 +3543,10 @@
                 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT. land_surface )  THEN
+             IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res land_surface = .TRUE.'
                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT. land_surface )  THEN
+             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res land_surface = .TRUE.'
@@ -3541,5 +3558,5 @@
                 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT. land_surface )  THEN
+             IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res land_surface = .TRUE.'
@@ -3556,20 +3573,20 @@
                 CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT. precipitation )  THEN
+             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT.  precipitation )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res precipitation = .TRUE.'
                 CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT. humidity )  THEN
+             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT.  humidity )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res humidity = .TRUE.'
                 CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT. land_surface )  THEN
+             IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT.  land_surface )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
                                  'res land_surface = .TRUE.'
                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )  &
+             IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )   &
              THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -3577,5 +3594,5 @@
                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT. land_surface ) &
+             IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
              THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -3583,5 +3600,5 @@
                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )  &
+             IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
              THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -3589,5 +3606,5 @@
                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT. land_surface ) &
+             IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface ) &
              THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -3595,5 +3612,5 @@
                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
              ENDIF
-             IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT. land_surface ) &
+             IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface ) &
              THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -3757,5 +3774,5 @@
 !-- Check mask conditions
     DO mid = 1, max_masks
-       IF ( data_output_masks(mid,1) /= ' ' .OR.                               &
+       IF ( data_output_masks(mid,1) /= ' '  .OR.                              &
             data_output_masks_user(mid,1) /= ' ' ) THEN
           masks = masks + 1
@@ -3763,5 +3780,5 @@
     ENDDO
     
-    IF ( masks < 0 .OR. masks > max_masks )  THEN
+    IF ( masks < 0  .OR.  masks > max_masks )  THEN
        WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
             '<= ', max_masks, ' (=max_masks)'
@@ -3908,6 +3925,6 @@
 !
 !-- 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 = "' //    &
@@ -3919,5 +3936,5 @@
 !-- Determine upper and lower hight level indices for random perturbations
     IF ( disturbance_level_b == -9999999.9_wp )  THEN
-       IF ( ocean ) THEN
+       IF ( ocean )  THEN
           disturbance_level_b     = zu((nzt*2)/3)
           disturbance_level_ind_b = ( nzt * 2 ) / 3
@@ -4182,5 +4199,5 @@
 !
 !-- Check pressure gradient conditions
-    IF ( dp_external .AND. conserve_volume_flow )  THEN
+    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.'
@@ -4188,5 +4205,5 @@
     ENDIF
     IF ( dp_external )  THEN
-       IF ( dp_level_b < zu(nzb) .OR. dp_level_b > zu(nzt) )  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'
@@ -4199,5 +4216,5 @@
        ENDIF
     ENDIF
-    IF ( ANY( dpdxy /= 0.0_wp ) .AND. .NOT. dp_external )  THEN
+    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'
@@ -4230,6 +4247,6 @@
        ENDIF
     ENDIF
-    IF ( ( u_bulk /= 0.0_wp .OR. v_bulk /= 0.0_wp ) .AND.                      &
-         ( .NOT. conserve_volume_flow .OR.                                     &
+    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 ',          &
@@ -4270,5 +4287,5 @@
 !
 !-- Check nudging and large scale forcing from external file
-    IF ( nudging .AND. ( .NOT. large_scale_forcing ) )  THEN
+    IF ( nudging  .AND.  (  .NOT.  large_scale_forcing ) )  THEN
        message_string = 'Nudging requires large_scale_forcing = .T.. &'//      &
                         'Surface fluxes and geostrophic wind should be &'//    &
@@ -4277,5 +4294,5 @@
     ENDIF
 
-    IF ( large_scale_forcing .AND. ( bc_lr /= 'cyclic'  .OR.                   &
+    IF ( large_scale_forcing  .AND.  ( bc_lr /= 'cyclic'  .OR.                 &
                                     bc_ns /= 'cyclic' ) )  THEN
        message_string = 'Non-cyclic lateral boundaries do not allow for &' //  &
@@ -4284,5 +4301,5 @@
      ENDIF
 
-    IF ( large_scale_forcing .AND. ( .NOT. humidity ) )  THEN
+    IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
        message_string = 'The usage of large scale forcing from external &'//   & 
                         'file LSF_DATA requires humidity = .T..'
@@ -4290,5 +4307,5 @@
      ENDIF
 
-    IF ( large_scale_forcing .AND. topography /= 'flat' )  THEN
+    IF ( large_scale_forcing  .AND.  topography /= 'flat' )  THEN
        message_string = 'The usage of large scale forcing from external &'//   & 
                         'file LSF_DATA is not implemented for non-flat topography'
@@ -4296,5 +4313,5 @@
     ENDIF
 
-    IF ( large_scale_forcing .AND.  ocean  )  THEN
+    IF ( large_scale_forcing  .AND.  ocean  )  THEN
        message_string = 'The usage of large scale forcing from external &'//   & 
                         'file LSF_DATA is not implemented for ocean runs'
@@ -4320,7 +4337,7 @@
 !-- Check for valid setting of most_method
     IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
-         TRIM( most_method ) /= 'newton'  .AND.                                &
+         TRIM( most_method ) /= 'newton'    .AND.                              &
          TRIM( most_method ) /= 'lookup' )  THEN
-       message_string = 'most_method = "' // TRIM( most_method ) //      &
+       message_string = 'most_method = "' // TRIM( most_method ) //            &
                         '" is unknown'
        CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
@@ -4352,13 +4369,13 @@
        IF ( dt_do == 0.0_wp )  THEN
           IF ( dt_fixed )  THEN
-             WRITE( message_string, '(A,F9.4,A)' )  'Output at every '  //  &
-                    'timestep is desired (' // dt_do_name // ' = 0.0).&'//  &
-                    'Setting the output interval to the fixed timestep '//  &
+             WRITE( message_string, '(A,F9.4,A)' )  'Output at every '  //     &
+                    'timestep is desired (' // dt_do_name // ' = 0.0).&'//     &
+                    'Setting the output interval 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 ' //  &
+             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 )
