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 )
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 1787)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 1788)
@@ -19,6 +19,6 @@
 ! Current revisions:
 ! -----------------
+! Added output of z0q
 ! 
-!
 ! Former revisions:
 ! -----------------
@@ -129,5 +129,5 @@
     USE arrays_3d,                                                             &
         ONLY:  dzw, e, nr, ol, p, pt, q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws,  &
-               rho, sa, shf, tend, ts, u, us, v, vpt, w, z0, z0h, zu, zw
+               rho, sa, shf, tend, ts, u, us, v, vpt, w, z0, z0h, z0q, zu, zw
         
     USE averaging
@@ -1247,4 +1247,22 @@
                       DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) =  z0h_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'z0q*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  z0q(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  z0q_av(j,i)
                       ENDDO
                    ENDDO
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 1787)
+++ palm/trunk/SOURCE/header.f90	(revision 1788)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-!
+! Parameter dewfall removed
 !
 ! Former revisions:
@@ -248,5 +248,5 @@
    
     USE land_surface_model_mod,                                                &
-        ONLY:  conserve_water_content, dewfall, land_surface, nzb_soil,        &
+        ONLY:  conserve_water_content, land_surface, nzb_soil,                 &
                nzt_soil, root_fraction, soil_moisture, soil_temperature,       &
                soil_type, soil_type_name, veg_type, veg_type_name, zs
@@ -964,10 +964,4 @@
        ELSE
           WRITE( io, 441 )
-       ENDIF
-
-       IF ( dewfall )  THEN
-          WRITE( io, 442 )
-       ELSE
-          WRITE( io, 443 )
        ENDIF
 
@@ -2341,8 +2335,6 @@
             '       Root fraction: ',A,'  '/ &
             '       Gridpoint:     ',A)
-440 FORMAT (/' --> Dewfall is allowed (default)')
-441 FORMAT (' --> Dewfall is inhibited')
-442 FORMAT (' --> Soil bottom is closed (water content is conserved, default)')
-443 FORMAT (' --> Soil bottom is open (water content is not conserved)')
+440 FORMAT (' --> Soil bottom is closed (water content is conserved, default)')
+441 FORMAT (' --> Soil bottom is open (water content is not conserved)')
 444 FORMAT (//' Radiation model information:'/                                 &
               ' ----------------------------'/)
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 1787)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 1788)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! ------------------
-!
+! Added z0q. 
+! Syntax layout improved.
 !
 ! Former revisions:
@@ -320,72 +321,72 @@
 !
 !-- Allocate arrays
-    ALLOCATE( mean_surface_level_height(0:statistic_regions),               &
-              mean_surface_level_height_l(0:statistic_regions),             &
-              ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
-              ngp_3d(0:statistic_regions),                                  &
-              ngp_3d_inner(0:statistic_regions),                            &
-              ngp_3d_inner_l(0:statistic_regions),                          &
-              ngp_3d_inner_tmp(0:statistic_regions),                        &
-              sums_divnew_l(0:statistic_regions),                           &
+    ALLOCATE( mean_surface_level_height(0:statistic_regions),                  &
+              mean_surface_level_height_l(0:statistic_regions),                &
+              ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions),    &
+              ngp_3d(0:statistic_regions),                                     &
+              ngp_3d_inner(0:statistic_regions),                               &
+              ngp_3d_inner_l(0:statistic_regions),                             &
+              ngp_3d_inner_tmp(0:statistic_regions),                           &
+              sums_divnew_l(0:statistic_regions),                              &
               sums_divold_l(0:statistic_regions) )
     ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt) )
-    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                 &
-              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),               &
-              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),               &
-              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),             &
-              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),               &
-              sums(nzb:nzt+1,pr_palm+max_pr_user),                          &
-              sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),   &
-              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
-              sums_up_fraction_l(10,3,0:statistic_regions),                 &
-              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                &
+    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                    &
+              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),                  &
+              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),                  &
+              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),                &
+              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),                  &
+              sums(nzb:nzt+1,pr_palm+max_pr_user),                             &
+              sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),      &
+              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1),    &
+              sums_up_fraction_l(10,3,0:statistic_regions),                    &
+              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                   &
               ts_value(dots_max,0:statistic_regions) )
     ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
 
-    ALLOCATE( ol(nysg:nyng,nxlg:nxrg), shf(nysg:nyng,nxlg:nxrg),      &
-              ts(nysg:nyng,nxlg:nxrg), tswst(nysg:nyng,nxlg:nxrg),    &
-              us(nysg:nyng,nxlg:nxrg), usws(nysg:nyng,nxlg:nxrg),     &
-              uswst(nysg:nyng,nxlg:nxrg), vsws(nysg:nyng,nxlg:nxrg),  &
-              vswst(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg),    &
-              z0h(nysg:nyng,nxlg:nxrg) )
-
-    ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),         &
-              kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
-              km(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
-              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
+    ALLOCATE( ol(nysg:nyng,nxlg:nxrg), shf(nysg:nyng,nxlg:nxrg),               &
+              ts(nysg:nyng,nxlg:nxrg), tswst(nysg:nyng,nxlg:nxrg),             &
+              us(nysg:nyng,nxlg:nxrg), usws(nysg:nyng,nxlg:nxrg),              &
+              uswst(nysg:nyng,nxlg:nxrg), vsws(nysg:nyng,nxlg:nxrg),           &
+              vswst(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg),             &
+              z0h(nysg:nyng,nxlg:nxrg), z0q(nysg:nyng,nxlg:nxrg) )
+
+    ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr),                                    &
+              kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                               &
+              km(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                               &
+              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
               tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 
 #if defined( __nopointer )
-    ALLOCATE( e(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
-              e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
-              pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-              u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
-              u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
-              v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
-              w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-              tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
-              tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-              tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+    ALLOCATE( e(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
+              e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                               &
+              pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+              u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
+              u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
+              v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                                &
+              w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+              tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
+              tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+              tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
               tw_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #else
-    ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
+    ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
+              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                              &
               w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-    IF ( .NOT. neutral )  THEN
+    IF (  .NOT.  neutral )  THEN
        ALLOCATE( pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ENDIF
@@ -412,9 +413,9 @@
     ENDIF
 
-    IF ( humidity  .OR.  passive_scalar ) THEN
+    IF ( humidity  .OR.  passive_scalar )  THEN
 !
 !--    2D-humidity/scalar arrays
-       ALLOCATE ( qs(nysg:nyng,nxlg:nxrg),   &
-                  qsws(nysg:nyng,nxlg:nxrg), &
+       ALLOCATE ( qs(nysg:nyng,nxlg:nxrg),                                     &
+                  qsws(nysg:nyng,nxlg:nxrg),                                   &
                   qswst(nysg:nyng,nxlg:nxrg) )
 
@@ -422,10 +423,10 @@
 !--    3D-humidity/scalar arrays
 #if defined( __nopointer )
-       ALLOCATE( q(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-                 q_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+       ALLOCATE( q(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                             &
+                 q_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
                  tq_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #else
-       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
-                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
+                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
                  q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #endif
@@ -440,5 +441,5 @@
 #endif
 
-          IF ( cloud_physics ) THEN
+          IF ( cloud_physics )  THEN
 
 !
@@ -451,5 +452,5 @@
 !
 !--          Precipitation amount and rate (only needed if output is switched)
-             ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), &
+             ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg),              &
                        precipitation_rate(nysg:nyng,nxlg:nxrg) )
 
@@ -457,5 +458,5 @@
 !
 !--             1D-arrays 
-                ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), &
+                ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1),                 &
                            q_1d(nzb:nzt+1), qc_1d(nzb:nzt+1) ) 
 !
@@ -474,25 +475,25 @@
 !
 !--                2D-rain water content and rain drop concentration arrays
-                   ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),                 &
-                              qrsws(nysg:nyng,nxlg:nxrg),               &
-                              qrswst(nysg:nyng,nxlg:nxrg),              &
-                              nrs(nysg:nyng,nxlg:nxrg),                 &
-                              nrsws(nysg:nyng,nxlg:nxrg),               &
+                   ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),                        &
+                              qrsws(nysg:nyng,nxlg:nxrg),                      &
+                              qrswst(nysg:nyng,nxlg:nxrg),                     &
+                              nrs(nysg:nyng,nxlg:nxrg),                        &
+                              nrsws(nysg:nyng,nxlg:nxrg),                      &
                               nrswst(nysg:nyng,nxlg:nxrg) )
 !
 !--                3D-rain water content, rain drop concentration arrays
 #if defined( __nopointer )
-                   ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),         &
-                             nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
-                             qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),         &
-                             qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
-                             tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),      &
+                   ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
+                             nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                             qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
+                             qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                             tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),             &
                              tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #else
-                   ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
-                             nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
-                             nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
-                             qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
-                             qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                   ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                             nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                             nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                             qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                             qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
                              qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #endif
@@ -509,13 +510,13 @@
 !--          Liquid water content, change in liquid water content
 #if defined( __nopointer )
-             ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+             ALLOCATE ( ql(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                     &
                         ql_c(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #else
-             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                   &
                         ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #endif
 !
 !--          Real volume of particles (with weighting), volume of particles
-             ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+             ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                   &
                         ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
           ENDIF
@@ -526,17 +527,17 @@
 
     IF ( ocean )  THEN
-       ALLOCATE( saswsb(nysg:nyng,nxlg:nxrg), &
+       ALLOCATE( saswsb(nysg:nyng,nxlg:nxrg),                                  &
                  saswst(nysg:nyng,nxlg:nxrg) )
 #if defined( __nopointer )
-       ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-                 rho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
-                 sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
-                 sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
+       ALLOCATE( prho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
+                 rho(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                           &
+                 sa(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                            &
+                 sa_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
                  tsa_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #else
-       ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
-                 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-                 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
-                 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
+       ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                        &
+                 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                         &
+                 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
+                 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                          &
                  sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
        prho => prho_1
@@ -553,5 +554,5 @@
 !-- 3D-array for storing the dissipation, needed for calculating the sgs
 !-- particle velocities
-    IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence  .OR.  &
+    IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence  .OR.      &
          num_acc_per_node > 0 )  THEN
        ALLOCATE( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
@@ -559,5 +560,5 @@
 
     IF ( dt_dosp /= 9999999.9_wp )  THEN
-       ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
+       ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ),                             &
                  spectrum_y( 1:ny/2, 1:10, 1:10 ) )
        spectrum_x = 0.0_wp
@@ -595,29 +596,29 @@
 !-- are needed for radiation boundary conditions
     IF ( outflow_l )  THEN
-       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), &
-                 v_m_l(nzb:nzt+1,nysg:nyng,0:1), &
+       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2),                               &
+                 v_m_l(nzb:nzt+1,nysg:nyng,0:1),                               &
                  w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
     ENDIF
     IF ( outflow_r )  THEN
-       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
-                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
+       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
+                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx),                           &
                  w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
     ENDIF
     IF ( outflow_l  .OR.  outflow_r )  THEN
-       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), &
+       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng),           &
                  c_w(nzb:nzt+1,nysg:nyng) )
     ENDIF
     IF ( outflow_s )  THEN
-       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), &
-                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), &
+       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg),                               &
+                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg),                               &
                  w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
     ENDIF
     IF ( outflow_n )  THEN
-       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
-                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
+       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
+                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg),                           &
                  w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
     ENDIF
     IF ( outflow_s  .OR.  outflow_n )  THEN
-       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), &
+       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg),           &
                  c_w(nzb:nzt+1,nxlg:nxrg) )
     ENDIF
@@ -682,5 +683,5 @@
 !--        pres may need to be called  outside the RK-substeps. Antti will
 !--        check if this is really required.
-    ALLOCATE( weight_substep(1:intermediate_timestep_count_max), &
+    ALLOCATE( weight_substep(1:intermediate_timestep_count_max),               &
               weight_pres(0:intermediate_timestep_count_max) )
     weight_substep = 1.0_wp
@@ -691,5 +692,5 @@
 !
 !-- Initialize model variables
-    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
          TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
 !
@@ -819,5 +820,5 @@
           CALL location_message( 'finished', .TRUE. )
 
-       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
+       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
        THEN
 
@@ -825,5 +826,5 @@
 !
 !--       Overwrite initial profiles in case of nudging
-          IF ( nudging ) THEN
+          IF ( nudging )  THEN
              pt_init = ptnudge(:,1)
              u_init  = unudge(:,1)
@@ -833,5 +834,5 @@
              ENDIF
 
-             WRITE( message_string, * ) 'Initial profiles of u, v and ', &
+             WRITE( message_string, * ) 'Initial profiles of u, v and ',       &
                  'scalars from NUDGING_DATA are used.'
              CALL message( 'init_3d_model', 'PA0370', 0, 0, 0, 6, 0 )
@@ -936,5 +937,5 @@
           CALL location_message( 'finished', .TRUE. )
 
-       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
+       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 )                  &
        THEN
 
@@ -1018,5 +1019,6 @@
           DO j=0,ny
              DO i=0,nx
-                id_random_array(j,i) = random_dummy + 1E6 * ( ensemble_member_nr - 1000 )
+                id_random_array(j,i) = random_dummy + 1E6                      &
+                                       * ( ensemble_member_nr - 1000 )
                 random_dummy = random_dummy + 1
              END DO
@@ -1047,5 +1049,5 @@
 !
 !--             Initialize shf with data from external file LSF_DATA
-                IF ( large_scale_forcing .AND. lsf_surf ) THEN
+                IF ( large_scale_forcing .AND. lsf_surf )  THEN
                    CALL ls_forcing_surf ( simulated_time )
                 ENDIF
@@ -1068,5 +1070,5 @@
 !--       Determine the near-surface water flux
           IF ( humidity  .OR.  passive_scalar )  THEN
-             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
+             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
                   precipitation )  THEN
                 qrsws = 0.0_wp
@@ -1106,5 +1108,5 @@
              IF ( humidity  .OR.  passive_scalar )  THEN
                 qswst = 0.0_wp
-                IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
+                IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.           &
                      precipitation ) THEN
                    nrswst = 0.0_wp
@@ -1133,4 +1135,5 @@
           z0 = roughness_length
           z0h = z0h_factor * z0
+          z0q = z0h_factor * z0
 
           IF ( .NOT. constant_heatflux )  THEN 
@@ -1145,6 +1148,6 @@
 
           IF ( humidity  .OR.  passive_scalar )  THEN
-             IF ( .NOT. constant_waterflux )  qsws   = 0.0_wp
-             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
+             IF (  .NOT.  constant_waterflux )  qsws   = 0.0_wp
+             IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.              &
                   precipitation )  THEN
                 qrsws = 0.0_wp
@@ -1159,5 +1162,5 @@
 !--    the reference state will be set (overwritten) in init_ocean)
        IF ( use_single_reference_value )  THEN
-          IF ( .NOT. humidity )  THEN
+          IF (  .NOT.  humidity )  THEN
              ref_state(:) = pt_reference
           ELSE
@@ -1165,5 +1168,5 @@
           ENDIF
        ELSE
-          IF ( .NOT. humidity )  THEN
+          IF (  .NOT.  humidity )  THEN
              ref_state(:) = pt_init(:)
           ELSE
@@ -1216,5 +1219,5 @@
 !--    If required, change the surface humidity/scalar at the start of the 3D
 !--    run
-       IF ( ( humidity .OR. passive_scalar ) .AND. &
+       IF ( ( humidity .OR. passive_scalar ) .AND.                             &
             q_surface_initial_change /= 0.0_wp )  THEN
           q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
@@ -1229,5 +1232,5 @@
           tq_m = 0.0_wp
           q_p = q
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
+          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
                precipitation )  THEN
              tqr_m = 0.0_wp
@@ -1245,6 +1248,6 @@
        CALL location_message( 'finished', .TRUE. )
 
-    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.    &
-         TRIM( initializing_actions ) == 'cyclic_fill' )  &
+    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
+         TRIM( initializing_actions ) == 'cyclic_fill' )                       &
     THEN
 
@@ -1282,5 +1285,5 @@
 !
 !--    Initialization of the turbulence recycling method
-       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.  &
+       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
             turbulent_inflow )  THEN
 !
@@ -1346,6 +1349,6 @@
                 inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
              ELSE
-                WRITE( message_string, * ) 'inflow_damping_height must be ',&
-                     'explicitly specified because&the inversion height ', &
+                WRITE( message_string, * ) 'inflow_damping_height must be ',   &
+                     'explicitly specified because&the inversion height ',     &
                      'calculated by the prerun is zero.'
                 CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
@@ -1382,5 +1385,5 @@
 !
 !--    Inside buildings set velocities and TKE back to zero
-       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.  &
+       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
             topography /= 'flat' )  THEN
 !
@@ -1415,5 +1418,5 @@
        IF ( humidity  .OR.  passive_scalar )  THEN
           q_p = q
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
+          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
                precipitation )  THEN
              qr_p = qr
@@ -1430,5 +1433,5 @@
        IF ( humidity  .OR.  passive_scalar )  THEN
           tq_m = 0.0_wp
-          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
+          IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.                 &
                precipitation )  THEN
              tqr_m = 0.0_wp
@@ -1486,5 +1489,5 @@
              DO  j = nys, nyn
                 DO  k = nzb_2d(j,nx)+1, nzt
-                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
+                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
                                               u_init(k) * dzw(k)
                    volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
@@ -1522,5 +1525,5 @@
              DO  j = nys, nyn
                 DO  k = nzb_2d(j,nx)+1, nzt
-                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
+                   volume_flow_initial_l(1) = volume_flow_initial_l(1) +       &
                                               hom_sum(k,1,0) * dzw(k)
                    volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
@@ -1532,5 +1535,5 @@
              DO  i = nxl, nxr
                 DO  k = nzb_2d(ny,i)+1, nzt  
-                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
+                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
                                               hom_sum(k,2,0) * dzw(k)
                    volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
@@ -1568,5 +1571,5 @@
              DO  i = nxl, nxr
                 DO  k = nzb_2d(ny,i)+1, nzt  
-                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
+                   volume_flow_initial_l(2) = volume_flow_initial_l(2) +       &
                                               v(k,ny,i) * dzw(k)
                    volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
@@ -1605,6 +1608,6 @@
 !-- Impose random perturbation on the horizontal velocity field and then
 !-- remove the divergences from the velocity field at the initial stage
-    IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND.&
-         TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
+    IF ( create_disturbances  .AND.  disturbance_energy_limit /= 0.0_wp  .AND. &
+         TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
          TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
 
@@ -1718,10 +1721,10 @@
     rdf_sc = 0.0_wp
     IF ( rayleigh_damping_factor /= 0.0_wp )  THEN
-       IF ( .NOT. ocean )  THEN
+       IF (  .NOT.  ocean )  THEN
           DO  k = nzb+1, nzt
              IF ( zu(k) >= rayleigh_damping_height )  THEN
-                rdf(k) = rayleigh_damping_factor * &
+                rdf(k) = rayleigh_damping_factor *                             &
                       ( SIN( pi * 0.5_wp * ( zu(k) - rayleigh_damping_height ) &
-                                         / ( zu(nzt) - rayleigh_damping_height ) )&
+                             / ( zu(nzt) - rayleigh_damping_height ) )         &
                       )**2
              ENDIF
@@ -1730,7 +1733,7 @@
           DO  k = nzt, nzb+1, -1
              IF ( zu(k) <= rayleigh_damping_height )  THEN
-                rdf(k) = rayleigh_damping_factor * &
+                rdf(k) = rayleigh_damping_factor *                             &
                       ( SIN( pi * 0.5_wp * ( rayleigh_damping_height - zu(k) ) &
-                                         / ( rayleigh_damping_height - zu(nzb+1)))&
+                             / ( rayleigh_damping_height - zu(nzb+1) ) )       &
                       )**2
              ENDIF
@@ -1774,5 +1777,5 @@
              ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5_wp *              &
                             REAL( pt_damping_width - i * dx, KIND=wp ) / (     &
-                            REAL( pt_damping_width, KIND=wp )            ) ) )**2 
+                            REAL( pt_damping_width, KIND=wp ) ) ) )**2 
           ENDIF
        ENDDO
@@ -1831,6 +1834,6 @@
 
     IF ( dots_num > dots_max )  THEN
-       WRITE( message_string, * ) 'number of time series quantities exceeds', &
-                                  ' its maximum of dots_max = ', dots_max,    &
+       WRITE( message_string, * ) 'number of time series quantities exceeds',  &
+                                  ' its maximum of dots_max = ', dots_max,     &
                                   ' &Please increase dots_max in modules.f90.'
        CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )    
@@ -1869,6 +1872,6 @@
 !--             All xy-grid points
                 ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
-                mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) + &
-                                                  zw(nzb_s_inner(j,i))
+                mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) &
+                                                  + zw(nzb_s_inner(j,i))
 !
 !--             xy-grid points above topography
@@ -1881,6 +1884,6 @@
 !
 !--             All grid points of the total domain above topography
-                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + &
-                                     ( nz - nzb_s_inner(j,i) + 2 )
+                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr)                        &
+                                     + ( nz - nzb_s_inner(j,i) + 2 )
              ENDIF
           ENDDO
@@ -1891,19 +1894,19 @@
 #if defined( __parallel )
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,   &
+    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,    &
                         comm2d, ierr )
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr,  &
+    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr,   &
                         MPI_INTEGER, MPI_SUM, comm2d, ierr )
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),         &
+    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),          &
                         (nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL, &
+    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL,  &
                         MPI_SUM, comm2d, ierr )
     ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( mean_surface_level_height_l(0),                       &
-                        mean_surface_level_height(0), sr, MPI_REAL,           &
+    CALL MPI_ALLREDUCE( mean_surface_level_height_l(0),                        &
+                        mean_surface_level_height(0), sr, MPI_REAL,            &
                         MPI_SUM, comm2d, ierr )
     mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh )
@@ -1924,9 +1927,9 @@
 !-- the respective subdomain lie below the surface topography
     ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   ) 
-    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),            &
+    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),             &
                            ngp_3d_inner(:) )
     ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 
 
-    DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l,      &
+    DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l,       &
                 ngp_3d_inner_l, ngp_3d_inner_tmp )
 
Index: palm/trunk/SOURCE/land_surface_model.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model.f90	(revision 1787)
+++ palm/trunk/SOURCE/land_surface_model.f90	(revision 1788)
@@ -19,6 +19,15 @@
 ! Current revisions:
 ! -----------------
+! Bugfix: calculate lambda_surface based on temperature gradient between skin
+! layer and soil layer instead of Obukhov length
+! Changed: moved calculation of surface specific humidity to energy balance solver
+! New: water surfaces are available by using a fixed sea surface temperature.
+! The roughness lengths are calculated dynamically using the Charnock 
+! parameterization. This involves the new roughness length for moisture z0q.
+! New: modified solution of the energy balance solver and soil model for
+! paved surfaces (i.e. asphalt concrete).
+! Syntax layout improved.
+! Changed: parameter dewfall removed.
 ! 
-!
 ! Former revisions:
 ! -----------------
@@ -105,5 +114,5 @@
 !> @todo Consider partial absorption of the net shortwave radiation by the 
 !>       skin layer.
-!> @todo Allow for water surfaces
+!> @todo Improve surface water parameterization
 !> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, 
 !>       nzt_soil=3)).
@@ -118,5 +127,6 @@
  
     USE arrays_3d,                                                             &
-        ONLY:  hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h
+        ONLY:  hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h, &
+               z0q
 
     USE cloud_parameters,                                                      &
@@ -171,9 +181,15 @@
 !
 !-- LSM variables
-    INTEGER(iwp) :: veg_type  = 2, & !< vegetation type, 0: user-defined, 1-19: generic (see list)
-                    soil_type = 3    !< soil type, 0: user-defined, 1-6: generic (see list)
+    INTEGER(iwp) :: veg_type  = 2, & !< NAMELIST veg_type_2d
+                    soil_type = 3    !< NAMELIST soil_type_2d
+
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  soil_type_2d, &  !< soil type, 0: user-defined, 1-7: generic (see list)
+                                                  veg_type_2d      !< vegetation type, 0: user-defined, 1-19: generic (see list)
+
+    LOGICAL, DIMENSION(:,:), ALLOCATABLE :: water_surface,     & !< flag parameter for water surfaces (classes 14+15)
+                                            pave_surface,      & !< flag parameter for pavements (asphalt etc.) (class 20)
+                                            building_surface     !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
 
     LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
-               dewfall = .TRUE.,                 & !< allow/inhibit dewfall
                force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls
                land_surface = .FALSE.              !< flag parameter indicating wheather the lsm is used
@@ -200,4 +216,7 @@
                 m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
                 n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
+                pave_depth = 9999999.9_wp,              & !< depth of the pavement
+                pave_heat_capacity = 1.94E6_wp,         & !< volumetric heat capacity of pavement (e.g. roads)
+                pave_heat_conductivity = 1.00_wp,       & !< heat conductivity for pavements (e.g. roads)
                 q_s = 0.0_wp,                           & !< saturation specific humidity
                 residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
@@ -210,5 +229,6 @@
                 wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
                 z0_eb  = 9999999.9_wp,                  & !< NAMELIST z0 (lsm_par)
-                z0h_eb = 9999999.9_wp                     !< NAMELIST z0h (lsm_par)
+                z0h_eb = 9999999.9_wp,                  & !< NAMELIST z0h (lsm_par)
+                z0q_eb = 9999999.9_wp                     !< NAMELIST z0q (lsm_par)
 
     REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
@@ -295,7 +315,7 @@
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
-              lambda_h, &   !< heat conductivity of soil (?)                            
+              lambda_h, &   !< heat conductivity of soil (W/m/K)                            
               lambda_w, &   !< hydraulic diffusivity of soil (?)
-              gamma_w,  &   !< hydraulic conductivity of soil (?)
+              gamma_w,  &   !< hydraulic conductivity of soil (W/m/K)
               rho_c_total   !< volumetric heat capacity of the actual soil matrix (?) 
 
@@ -327,5 +347,5 @@
 !
 !-- Predefined Land surface classes (veg_type)
-    CHARACTER(26), DIMENSION(0:19), PARAMETER :: veg_type_name = (/ &
+    CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ &
                                    'user defined              ',    & ! 0 
                                    'crops, mixed farming      ',    & !  1
@@ -347,5 +367,6 @@
                                    'deciduous shrubs          ',    & ! 17
                                    'mixed forest/woodland     ',    & ! 18
-                                   'interrupted forest        '     & ! 19
+                                   'interrupted forest        ',    & ! 19
+                                   'pavements/roads           '     & ! 20
                                                                  /)
 
@@ -368,5 +389,5 @@
 !-- Land surface parameters I
 !--                          r_canopy_min,     lai,   c_veg,     g_d
-    REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: veg_pars = RESHAPE( (/ &
+    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ &
                                  180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
                                  110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
@@ -387,34 +408,36 @@
                                  225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
                                  250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
-                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp   & ! 19
-                                 /), (/ 4, 19 /) )
-
-!
-!-- Land surface parameters II          z0,         z0h
-    REAL(wp), DIMENSION(0:1,1:19), PARAMETER :: roughness_par = RESHAPE( (/ & 
-                                   0.25_wp,  0.25E-2_wp,                    & !  1
-                                   0.20_wp,  0.20E-2_wp,                    & !  2
-                                   2.00_wp,     2.00_wp,                    & !  3
-                                   2.00_wp,     2.00_wp,                    & !  4
-                                   2.00_wp,     2.00_wp,                    & !  5
-                                   2.00_wp,     2.00_wp,                    & !  6
-                                   0.47_wp,  0.47E-2_wp,                    & !  7
-                                  0.013_wp, 0.013E-2_wp,                    & !  8
-                                  0.034_wp, 0.034E-2_wp,                    & !  9
-                                    0.5_wp,  0.50E-2_wp,                    & ! 10
-                                   0.17_wp,  0.17E-2_wp,                    & ! 11
-                                 1.3E-3_wp,   1.3E-4_wp,                    & ! 12
-                                   0.83_wp,  0.83E-2_wp,                    & ! 13
-                                   0.00_wp,  0.00E-2_wp,                    & ! 14
-                                   0.00_wp,  0.00E-2_wp,                    & ! 15
-                                   0.10_wp,  0.10E-2_wp,                    & ! 16
-                                   0.25_wp,  0.25E-2_wp,                    & ! 17
-                                   2.00_wp,  2.00E-2_wp,                    & ! 18
-                                   1.10_wp,  1.10E-2_wp                     & ! 19
-                                 /), (/ 2, 19 /) )
+                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,  & ! 19
+                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp   & ! 20
+                                 /), (/ 4, 20 /) )
+
+!
+!-- Land surface parameters II          z0,         z0h,         z0q
+    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: roughness_par = RESHAPE( (/ & 
+                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & !  1
+                                   0.20_wp,  0.20E-2_wp,  0.20E-2_wp,       & !  2
+                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  3
+                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  4
+                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  5
+                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  6
+                                   0.47_wp,  0.47E-2_wp,  0.47E-2_wp,       & !  7
+                                  0.013_wp, 0.013E-2_wp, 0.013E-2_wp,       & !  8
+                                  0.034_wp, 0.034E-2_wp, 0.034E-2_wp,       & !  9
+                                    0.5_wp,  0.50E-2_wp,  0.50E-2_wp,       & ! 10
+                                   0.17_wp,  0.17E-2_wp,  0.17E-2_wp,       & ! 11
+                                 1.3E-3_wp,   1.3E-4_wp,   1.3E-4_wp,       & ! 12
+                                   0.83_wp,  0.83E-2_wp,  0.83E-2_wp,       & ! 13
+                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 14
+                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 15
+                                   0.10_wp,  0.10E-2_wp,  0.10E-2_wp,       & ! 16
+                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & ! 17
+                                   2.00_wp,  2.00E-2_wp,  2.00E-2_wp,       & ! 18
+                                   1.10_wp,  1.10E-2_wp,  1.10E-2_wp,       & ! 19
+                                 1.0E-4_wp,   1.0E-5_wp,   1.0E-5_wp        & ! 20
+                                 /), (/ 3, 20 /) )
 
 !
 !-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
-    REAL(wp), DIMENSION(0:2,1:19), PARAMETER :: surface_pars = RESHAPE( (/ &
+    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ &
                                       10.0_wp,       10.0_wp, 0.05_wp,     & !  1
                                       10.0_wp,       10.0_wp, 0.05_wp,     & !  2
@@ -430,15 +453,16 @@
                                       58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
                                       10.0_wp,       10.0_wp, 0.05_wp,     & ! 13 
-                                    1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 14
-                                    1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 15
+                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 14
+                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 15
                                       10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
                                       10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
                                       20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
-                                      20.0_wp,       15.0_wp, 0.03_wp      & ! 19
-                                      /), (/ 3, 19 /) )
+                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 19
+                                       0.0_wp,        0.0_wp, 0.00_wp      & ! 20
+                                      /), (/ 3, 20 /) )
 
 !
 !-- Root distribution (sum = 1)  level 1, level 2, level 3, level 4,
-    REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: root_distribution = RESHAPE( (/ &
+    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: root_distribution = RESHAPE( (/ &
                                  0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  1
                                  0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  2
@@ -459,6 +483,7 @@
                                  0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 17 
                                  0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 18
-                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 19
-                                 /), (/ 4, 19 /) )
+                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 19
+                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp             & ! 20
+                                 /), (/ 4, 20 /) )
 
 !
@@ -499,14 +524,15 @@
 !-- Public parameters, constants and initial values
     PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient,       &
-           conserve_water_content, dewfall, field_capacity,                    &
+           conserve_water_content, field_capacity,                             &
            f_shortwave_incoming, hydraulic_conductivity, init_lsm,             &
            init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable,    &
            land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model,  &
            lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance,          &
-           min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp,     &
+           min_soil_resistance, n_vangenuchten, pave_heat_capacity,            &
+           pave_depth, pave_heat_conductivity, residual_moisture, rho_cp,      &
            rho_lv, root_fraction, saturation_moisture, skip_time_do_lsm,       &
            soil_moisture, soil_temperature, soil_type, soil_type_name,         &
            vegetation_coverage, veg_type, veg_type_name, wilting_point, z0_eb, &
-           z0h_eb
+           z0h_eb, z0q_eb
 
 !
@@ -586,4 +612,5 @@
 !--    Allocate 2D vegetation model arrays
        ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( building_surface(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
@@ -601,4 +628,5 @@
        ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( pave_surface(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
@@ -613,4 +641,7 @@
        ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( soil_type_2d(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( veg_type_2d(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( water_surface(nysg:nyng,nxlg:nxrg) )
 
 #if ! defined( __nopointer )
@@ -650,5 +681,5 @@
 !
 !--    If no cloud physics is used, rho_surface has not been calculated before
-       IF ( .NOT. cloud_physics )  THEN
+       IF (  .NOT.  cloud_physics )  THEN
           rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
        ENDIF
@@ -762,4 +793,6 @@
        ENDIF    
 
+!
+!--    Map values to the respective 2D arrays
        alpha_vg      = alpha_vangenuchten
        l_vg          = l_vangenuchten
@@ -796,6 +829,6 @@
 
 !
-!--       Set artifical values for ts and us so that r_a has its initial value for 
-!--       the first time step
+!--       Set artifical values for ts and us so that r_a has its initial value 
+!--       for the first time step
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
@@ -864,5 +897,8 @@
              z0h_eb = roughness_par(1,veg_type)
           ENDIF
-          z0h_factor = z0h_eb / z0_eb
+          IF ( z0q_eb == 9999999.9_wp )  THEN
+             z0q_eb = roughness_par(2,veg_type)
+          ENDIF
+          z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp )
 
           IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
@@ -881,9 +917,29 @@
              z0h_eb = z0_eb * z0h_factor
           ENDIF
+          IF ( z0q_eb == 9999999.9_wp )  THEN
+             z0q_eb = z0_eb * z0h_factor
+          ENDIF
 
        ENDIF
 
 !
-!--    Initialize vegetation
+!--    For surfaces covered with pavement, set depth of the pavement (with dry
+!--    soil below). The depth must be greater than the first soil layer depth
+       IF ( veg_type == 20 )  THEN
+          IF ( pave_depth == 9999999.9_wp )  THEN
+             pave_depth = zs(nzb_soil)  
+          ELSE
+             pave_depth = MAX( zs(nzb_soil), pave_depth )
+          ENDIF
+       ENDIF
+
+!
+!--    Map vegetation and soil types to 2D array to allow for heterogeneous 
+!--    surfaces via user interface see below
+       veg_type_2d = veg_type
+       soil_type_2d = soil_type
+
+!
+!--    Map vegetation parameters to the respective 2D arrays
        r_canopy_min         = min_canopy_resistance
        lai                  = leaf_area_index
@@ -895,4 +951,5 @@
        z0                   = z0_eb
        z0h                  = z0h_eb
+       z0q                  = z0q_eb
 
 !
@@ -900,4 +957,22 @@
        CALL user_init_land_surface
 
+!
+!--    Set flag parameter if vegetation type was set to a water surface. Also
+!--    set temperature to a constant value in all "soil" layers.
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             IF ( veg_type_2d(j,i) == 14  .OR.  veg_type_2d(j,i) == 15 )  THEN
+                water_surface(j,i) = .TRUE.
+             ELSEIF ( veg_type_2d(j,i) == 20 )  THEN
+                pave_surface(j,i) = .TRUE.
+                m_soil(:,j,i) = 0.0_wp
+             ENDIF
+
+          ENDDO
+       ENDDO
+
+!
+!--    Calculate new roughness lengths (for water surfaces only)
+       CALL calc_z0_water_surface
 
        t_soil_p    = t_soil
@@ -935,5 +1010,6 @@
        INTEGER(iwp) ::  k, ks     !< running index
 
-       REAL(wp) :: f1,          & !< resistance correction term 1
+       REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
+                   f1,          & !< resistance correction term 1
                    f2,          & !< resistance correction term 2
                    f3,          & !< resistance correction term 3
@@ -965,9 +1041,19 @@
 
 !
-!--          Set lambda_surface according to stratification
-             IF ( ol(j,i) >= 0.0_wp )  THEN
-                lambda_surface = lambda_surface_s(j,i)
+!--          Set lambda_surface according to stratification between skin layer and soil
+             IF (  .NOT.  pave_surface(j,i) )  THEN
+
+                c_surface_tmp = c_surface
+
+                IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
+                   lambda_surface = lambda_surface_s(j,i)
+                ELSE
+                   lambda_surface = lambda_surface_u(j,i)
+                ENDIF
              ELSE
-                lambda_surface = lambda_surface_u(j,i)
+
+                c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
+                lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
+
              ENDIF
 
@@ -1016,5 +1102,5 @@
              ENDDO 
 
-             IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) )  THEN
+             IF ( m_total > m_wilt(j,i)  .AND.  m_total < m_fc(j,i) )  THEN
                 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
              ELSEIF ( m_total >= m_fc(j,i) )  THEN
@@ -1050,5 +1136,5 @@
 !--          Third step: calculate bare soil resistance r_soil. The Clapp & 
 !--          Hornberger parametrization does not consider c_veg.
-             IF ( soil_type /= 7 )  THEN
+             IF ( soil_type_2d(j,i) /= 7 )  THEN
                 m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
                         m_res(j,i)
@@ -1064,8 +1150,18 @@
 
 !
-!--          Calculate fraction of liquid water reservoir
-             m_liq_eb_max = m_max_depth * lai(j,i)
-             c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp))
-             
+!--          Calculate the maximum possible liquid water amount on plants and
+!--          bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
+!--          assumed, while paved surfaces might hold up 1 mm of water. The 
+!--          liquid water fraction for paved surfaces is calculated after 
+!--          Noilhan & Planton (1989), while the ECMWF formulation is used for
+!--          vegetated surfaces and bare soils.
+             IF ( pave_surface(j,i) )  THEN
+                m_liq_eb_max = m_max_depth * 5.0_wp
+                c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
+             ELSE
+                m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)           &
+                               + (1.0_wp - c_veg(j,i)) )
+                c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
+             ENDIF
 
 !
@@ -1076,5 +1172,5 @@
 !--          In case of dewfall, set evapotranspiration to zero
 !--          All super-saturated water is then removed from the air
-             IF ( humidity .AND. q_s <= qv1 )  THEN
+             IF ( humidity  .AND.  q_s <= qv1 )  THEN
                 r_canopy(j,i) = 0.0_wp
                 r_soil(j,i)   = 0.0_wp
@@ -1083,12 +1179,21 @@
 !
 !--          Calculate coefficients for the total evapotranspiration 
-             f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
-                           (r_a(j,i) + r_canopy(j,i))
-
-             f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
+!--          In case of water surface, set vegetation and soil fluxes to zero.
+!--          For pavements, only evaporation of liquid water is possible.
+             IF ( water_surface(j,i) )  THEN
+                f_qsws_veg  = 0.0_wp
+                f_qsws_soil = 0.0_wp
+                f_qsws_liq  = rho_lv / r_a(j,i)
+             ELSEIF ( pave_surface (j,i) )  THEN
+                f_qsws_veg  = 0.0_wp
+                f_qsws_soil = 0.0_wp
+                f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
+             ELSE
+                f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/     &
+                              (r_a(j,i) + r_canopy(j,i))
+                f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +     &
                                                           r_soil(j,i))
-             f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
-
-
+                f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
+             ENDIF
 !
 !--          If soil moisture is below wilting point, plants do no longer
@@ -1098,14 +1203,4 @@
 !              ENDIF
 
-!
-!--          If dewfall is deactivated, vegetation, soil and liquid water
-!--          reservoir are not allowed to take up water from the super-saturated
-!--          air.
-             IF ( humidity .AND. q_s <= qv1 .AND. .NOT. dewfall )  THEN
-                f_qsws_veg  = 0.0_wp
-                f_qsws_soil = 0.0_wp
-                f_qsws_liq  = 0.0_wp
-             ENDIF 
-
              f_shf  = rho_cp / r_a(j,i)
              f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
@@ -1123,4 +1218,6 @@
              rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
 
+!
+!--          Calculate new skin temperature
              IF ( humidity )  THEN
 #if defined ( __rrtmg )
@@ -1186,10 +1283,11 @@
 !--          Implicit solution when the surface layer has no heat capacity,
 !--          otherwise use RK3 scheme.
-             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface *        &
-                                t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d&
-                                * tsc(2) ) 
+             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *    &
+                                t_surface(j,i) ) / ( c_surface_tmp + coef_2    &
+                                * dt_3d * tsc(2) ) 
+
 !
 !--          Add RK3 term
-             IF ( c_surface /= 0.0_wp )  THEN
+             IF ( c_surface_tmp /= 0.0_wp )  THEN
 
                 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
@@ -1211,5 +1309,4 @@
                    ENDIF
                 ENDIF
-
              ENDIF
 
@@ -1223,5 +1320,5 @@
 !--          often reached, when no oscillations would occur (causes immense
 !--          computing time for the radiation code).
-             IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND.      &
+             IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.     &
                   unscheduled_radiation_calls )  THEN
                 force_radiation_call_l = .TRUE.
@@ -1249,5 +1346,4 @@
 #endif
 
-
              ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
                               - t_soil(nzb_soil,j,i))
@@ -1275,4 +1371,5 @@
                                     * t_surface_p(j,i) )
              ENDIF
+
 !
 !--          Calculate the true surface resistance
@@ -1280,5 +1377,5 @@
                 r_s(j,i) = 1.0E10_wp
              ELSE
-                r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                   &
+                r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                    &
                            * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
                            / qsws_eb(j,i) - r_a(j,i)
@@ -1297,5 +1394,7 @@
 !
 !--                Add precipitation to liquid water reservoir, if possible.
-!--                Otherwise, add the water to bare soil fraction.
+!--                Otherwise, add the water to soil. In case of
+!--                pavements, the exceeding water amount is implicitely removed 
+!--                as runoff as qsws_soil_eb is then not used in the soil model
                    IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
                       qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
@@ -1376,5 +1475,5 @@
 !--    Make a logical OR for all processes. Force radiation call if at
 !--    least one processor reached the threshold change in skin temperature
-       IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count      &
+       IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count    &
             == intermediate_timestep_count_max-1 )  THEN
 #if defined( __parallel )
@@ -1388,4 +1487,14 @@
        ENDIF
 
+!
+!--    Calculate surface specific humidity
+       IF ( humidity )  THEN
+          CALL calc_q_surface
+       ENDIF
+
+!
+!--    Calculate new roughness lengths (for water surfaces only)
+       CALL calc_z0_water_surface
+
 
     END SUBROUTINE lsm_energy_balance
@@ -1415,238 +1524,266 @@
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             DO  k = nzb_soil, nzt_soil
-!
-!--             Calculate volumetric heat capacity of the soil, taking into 
-!--             account water content 
-                rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i))       &
-                                     + rho_c_water * m_soil(k,j,i))
-
-!
-!--             Calculate soil heat conductivity at the center of the soil
-!--             layers
-                lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *          &
-                               lambda_h_water ** m_soil(k,j,i)
-
-                ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
-
-                lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +          &
-                                 lambda_h_dry
-
-             ENDDO
-
-!
-!--          Calculate soil heat conductivity (lambda_h) at the _stag level 
-!--          using linear interpolation
-             DO  k = nzb_soil, nzt_soil-1
-                lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) * 0.5_wp
-             ENDDO
-             lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
-
-!
-!--          Prognostic equation for soil temperature t_soil
-             tend(:) = 0.0_wp
-             tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *             &
-                       ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)     &
-                         - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)       &
-                         + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
-
-
-
-
-             DO  k = nzb_soil+1, nzt_soil
-                tend(k) = (1.0_wp/rho_c_total(k,j,i))                          &
-                          * (   lambda_h(k,j,i)                                &
-                              * ( t_soil(k+1,j,i) - t_soil(k,j,i) )            &
-                              * ddz_soil(k+1)                                  &
-                              - lambda_h(k-1,j,i)                              &
-                              * ( t_soil(k,j,i) - t_soil(k-1,j,i) )            &
-                              * ddz_soil(k)                                    &
-                            ) * ddz_soil_stag(k)
-             ENDDO
-
-             t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)   &
-                                             + dt_3d * ( tsc(2)                &
-                                             * tend(:) + tsc(3)                &
-                                             * tt_soil_m(:,j,i) )   
-
-!
-!--          Calculate t_soil tendencies for the next Runge-Kutta step
-             IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                IF ( intermediate_timestep_count == 1 )  THEN
-                   DO  k = nzb_soil, nzt_soil
-                      tt_soil_m(k,j,i) = tend(k)
-                   ENDDO
-                ELSEIF ( intermediate_timestep_count <                         &
-                         intermediate_timestep_count_max )  THEN
-                   DO  k = nzb_soil, nzt_soil
-                      tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
-                                         * tt_soil_m(k,j,i)
-                   ENDDO
-                ENDIF
-             ENDIF
-
-
-             DO  k = nzb_soil, nzt_soil
-
-!
-!--             Calculate soil diffusivity at the center of the soil layers
-                lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat          &
-                                  / m_sat(j,i) ) * ( MAX(m_soil(k,j,i),        &
-                                  m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp)
-
-!
-!--             Parametrization of Van Genuchten
-                IF ( soil_type /= 7 )  THEN
-!
-!--                Calculate the hydraulic conductivity after Van Genuchten 
-!--                (1980)
-                   h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -       &
-                              MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i)   &
-                              / (n_vg(j,i)-1.0_wp)) - 1.0_wp                   &
-                          )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i)
-
-
-                   gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +            &
-                                   (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp   &
-                                   -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg    &
-                                   )**(n_vg(j,i)-1.0_wp))**2 )                 &
-                                   / ( (1.0_wp + (alpha_vg(j,i)*h_vg           &
-                                   )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) &
-                                   *(l_vg(j,i) + 2.0_wp)) )
-
-!
-!--             Parametrization of Clapp & Hornberger
-                ELSE
-                   gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i)           &
-                                   / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
-                ENDIF
-
-             ENDDO
-
-
-             IF ( humidity )  THEN
-!
-!--             Calculate soil diffusivity (lambda_w) at the _stag level 
-!--             using linear interpolation. To do: replace this with
-!--             ECMWF-IFS Eq. 8.81
+
+             IF ( pave_surface(j,i) )  THEN
+                rho_c_total(nzb_soil,j,i) = pave_heat_capacity
+                lambda_temp(nzb_soil)     = pave_heat_conductivity
+             ENDIF
+
+             IF (  .NOT.  water_surface(j,i) )  THEN
+                DO  k = nzb_soil, nzt_soil
+
+
+                   IF ( pave_surface(j,i)  .AND.  zs(k) <= pave_depth )  THEN
+                    
+                      rho_c_total(k,j,i) = pave_heat_capacity
+                      lambda_temp(k)     = pave_heat_conductivity   
+
+                   ELSE            
+!
+!--                   Calculate volumetric heat capacity of the soil, taking 
+!--                   into account water content 
+                      rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) &
+                                           + rho_c_water * m_soil(k,j,i))
+
+!
+!--                   Calculate soil heat conductivity at the center of the soil
+!--                   layers
+                      lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *    &
+                                     lambda_h_water ** m_soil(k,j,i)
+
+                      ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
+
+                      lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +    &
+                                       lambda_h_dry
+                   ENDIF
+
+                ENDDO
+
+!
+!--             Calculate soil heat conductivity (lambda_h) at the _stag level 
+!--             using linear interpolation. For pavement surface, the
+!--             true pavement depth is considered
                 DO  k = nzb_soil, nzt_soil-1
-                     
-                   lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )     &
-                                     * 0.5_wp
-                   gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )       &
-                                     * 0.5_wp
-
+                   IF ( pave_surface(j,i)  .AND.  zs(k)   < pave_depth         &
+                                           .AND.  zs(k+1) > pave_depth )  THEN
+                      lambda_h(k,j,i) = ( pave_depth - zs(k) ) / dz_soil(k+1)  &
+                                        * lambda_temp(k)                       &
+                                        + ( 1.0_wp - ( pave_depth - zs(k) )    &
+                                        / dz_soil(k+1) ) * lambda_temp(k+1)
+                   ELSE
+                      lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
+                                        * 0.5_wp
+                   ENDIF
                 ENDDO
-
-!
-!
-!--             In case of a closed bottom (= water content is conserved), set
-!--             hydraulic conductivity to zero to that no water will be lost
-!--             in the bottom layer.
-                IF ( conserve_water_content )  THEN
-                   gamma_w(nzt_soil,j,i) = 0.0_wp
-                ELSE
-                   gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
-                ENDIF     
-
-!--             The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
-!--             ensures the mass conservation for water. The transpiration of
-!--             plants equals the cumulative withdrawals by the roots in the 
-!--             soil. The scheme takes into account the availability of water
-!--             in the soil layers as well as the root fraction in the
-!--             respective layer. Layer with moisture below wilting point will
-!--             not contribute, which reflects the preference of plants to
-!--             take water from moister layers.
-
-!
-!--             Calculate the root extraction (ECMWF 7.69, the sum of root_extr
-!--             = 1). The energy balance solver guarantees a positive
-!--             transpiration, so that there is no need for an additional check.
-                DO  k = nzb_soil, nzt_soil
-                    IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
-                       m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
-                    ENDIF
-                ENDDO  
-
-                IF ( m_total > 0.0_wp )  THEN
-                   DO  k = nzb_soil, nzt_soil
-                      IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
-                         root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
-                      ELSE
-                         root_extr(k) = 0.0_wp
-                      ENDIF
-                   ENDDO
-                ENDIF
-
-!
-!--             Prognostic equation for soil water content m_soil
+                lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
+
+
+
+
+!
+!--             Prognostic equation for soil temperature t_soil
                 tend(:) = 0.0_wp
-                tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (                  &
-                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
-                            * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
-                            root_extr(nzb_soil) * qsws_veg_eb(j,i)             &
-                            + qsws_soil_eb(j,i) ) * drho_l_lv )                &
-                            * ddz_soil_stag(nzb_soil)
-
-                DO  k = nzb_soil+1, nzt_soil-1
-                   tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)             &
-                             - m_soil(k,j,i) ) * ddz_soil(k+1) - gamma_w(k,j,i)&
-                             - lambda_w(k-1,j,i) * (m_soil(k,j,i) -            &
-                             m_soil(k-1,j,i)) * ddz_soil(k)                    &
-                             + gamma_w(k-1,j,i) - (root_extr(k)                &
-                             * qsws_veg_eb(j,i) * drho_l_lv)                   &
-                             ) * ddz_soil_stag(k)
+
+                tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *          &
+                          ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)  &
+                            - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)    &
+                            + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
+
+                DO  k = nzb_soil+1, nzt_soil
+                   tend(k) = (1.0_wp/rho_c_total(k,j,i))                       &
+                             * (   lambda_h(k,j,i)                             &
+                                 * ( t_soil(k+1,j,i) - t_soil(k,j,i) )         &
+                                 * ddz_soil(k+1)                               &
+                                 - lambda_h(k-1,j,i)                           &
+                                 * ( t_soil(k,j,i) - t_soil(k-1,j,i) )         &
+                                 * ddz_soil(k)                                 &
+                               ) * ddz_soil_stag(k)
 
                 ENDDO
-                tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                     &
-                                        - lambda_w(nzt_soil-1,j,i)             &
-                                        * (m_soil(nzt_soil,j,i)                &
-                                        - m_soil(nzt_soil-1,j,i))              &
-                                        * ddz_soil(nzt_soil)                   &
-                                        + gamma_w(nzt_soil-1,j,i) - (          &
-                                          root_extr(nzt_soil)                  &
-                                        * qsws_veg_eb(j,i) * drho_l_lv  )      &
-                                      ) * ddz_soil_stag(nzt_soil)             
-
-                m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
-                                                + dt_3d * ( tsc(2) * tend(:)   &
-                                                + tsc(3) * tm_soil_m(:,j,i) )   
-
-!
-!--             Account for dry soils (find a better solution here!)
-                DO  k = nzb_soil, nzt_soil
-                   IF ( m_soil_p(k,j,i) < 0.0_wp )  m_soil_p(k,j,i) = 0.0_wp
-                ENDDO
-
-
-
-
-!
-!--             Calculate m_soil tendencies for the next Runge-Kutta step
+
+                t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)&
+                                                  + dt_3d * ( tsc(2)           &
+                                                  * tend(nzb_soil:nzt_soil)    & 
+                                                  + tsc(3)                     &
+                                                  * tt_soil_m(:,j,i) )   
+
+!
+!--             Calculate t_soil tendencies for the next Runge-Kutta step
                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
                    IF ( intermediate_timestep_count == 1 )  THEN
                       DO  k = nzb_soil, nzt_soil
-                         tm_soil_m(k,j,i) = tend(k)
+                         tt_soil_m(k,j,i) = tend(k)
                       ENDDO
                    ELSEIF ( intermediate_timestep_count <                      &
                             intermediate_timestep_count_max )  THEN
                       DO  k = nzb_soil, nzt_soil
-                         tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
-                                     * tm_soil_m(k,j,i)
+                         tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
+                                         * tt_soil_m(k,j,i)
                       ENDDO
                    ENDIF
                 ENDIF
 
+
+                DO  k = nzb_soil, nzt_soil
+
+!
+!--                Calculate soil diffusivity at the center of the soil layers
+                   lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat       &
+                                     / m_sat(j,i) ) * ( MAX( m_soil(k,j,i),    &
+                                     m_wilt(j,i) ) / m_sat(j,i) )**(           &
+                                     b_ch + 2.0_wp )
+
+!
+!--                Parametrization of Van Genuchten
+                   IF ( soil_type /= 7 )  THEN
+!
+!--                   Calculate the hydraulic conductivity after Van Genuchten 
+!--                   (1980)
+                      h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -    &
+                                 MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**(      &
+                                 n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp  &
+                             )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i)
+
+
+                      gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +         &
+                                      ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**(  &
+                                      1.0_wp - 1.0_wp / n_vg(j,i) ) - (        &
+                                      alpha_vg(j,i) * h_vg )**( n_vg(j,i)      &
+                                      - 1.0_wp) )**2 )                         &
+                                      / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg    &
+                                      )**n_vg(j,i) )**( ( 1.0_wp  - 1.0_wp     &
+                                      / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) )
+
+!
+!--                Parametrization of Clapp & Hornberger
+                   ELSE
+                      gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i)       &
+                                      / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
+                   ENDIF
+
+                ENDDO
+
+!
+!--             Prognostic equation for soil moisture content. Only performed,
+!--             when humidity is enabled in the atmosphere and the surface type
+!--             is not pavement (implies dry soil below).
+                IF ( humidity  .AND.  .NOT.  pave_surface(j,i) )  THEN
+!
+!--                Calculate soil diffusivity (lambda_w) at the _stag level 
+!--                using linear interpolation. To do: replace this with
+!--                ECMWF-IFS Eq. 8.81
+                   DO  k = nzb_soil, nzt_soil-1
+                     
+                      lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
+                                        * 0.5_wp
+                      gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )    &
+                                        * 0.5_wp
+
+                   ENDDO
+
+!
+!
+!--                In case of a closed bottom (= water content is conserved), set
+!--                hydraulic conductivity to zero to that no water will be lost
+!--                in the bottom layer.
+                   IF ( conserve_water_content )  THEN
+                      gamma_w(nzt_soil,j,i) = 0.0_wp
+                   ELSE
+                      gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
+                   ENDIF     
+
+!--                The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
+!--                ensures the mass conservation for water. The transpiration of
+!--                plants equals the cumulative withdrawals by the roots in the 
+!--                soil. The scheme takes into account the availability of water
+!--                in the soil layers as well as the root fraction in the
+!--                respective layer. Layer with moisture below wilting point will
+!--                not contribute, which reflects the preference of plants to
+!--                take water from moister layers.
+
+!
+!--                Calculate the root extraction (ECMWF 7.69, the sum of root_extr
+!--                = 1). The energy balance solver guarantees a positive
+!--                transpiration, so that there is no need for an additional check.
+                   DO  k = nzb_soil, nzt_soil
+                       IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
+                          m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
+                       ENDIF
+                   ENDDO  
+
+                   IF ( m_total > 0.0_wp )  THEN
+                      DO  k = nzb_soil, nzt_soil
+                         IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
+                            root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
+                         ELSE
+                            root_extr(k) = 0.0_wp
+                         ENDIF
+                      ENDDO
+                   ENDIF
+
+!
+!--                Prognostic equation for soil water content m_soil.
+                   tend(:) = 0.0_wp
+
+                   tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (               &
+                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
+                            * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
+                               root_extr(nzb_soil) * qsws_veg_eb(j,i)          &
+                               + qsws_soil_eb(j,i) ) * drho_l_lv )             &
+                               * ddz_soil_stag(nzb_soil)
+
+                   DO  k = nzb_soil+1, nzt_soil-1
+                      tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)          &
+                                - m_soil(k,j,i) ) * ddz_soil(k+1)              &
+                                - gamma_w(k,j,i)                               &
+                                - lambda_w(k-1,j,i) * (m_soil(k,j,i) -         &
+                                m_soil(k-1,j,i)) * ddz_soil(k)                 &
+                                + gamma_w(k-1,j,i) - (root_extr(k)             &
+                                * qsws_veg_eb(j,i) * drho_l_lv)                &
+                                ) * ddz_soil_stag(k)
+
+                   ENDDO
+                   tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                  &
+                                           - lambda_w(nzt_soil-1,j,i)          &
+                                           * (m_soil(nzt_soil,j,i)             &
+                                           - m_soil(nzt_soil-1,j,i))           &
+                                           * ddz_soil(nzt_soil)                &
+                                           + gamma_w(nzt_soil-1,j,i) - (       &
+                                             root_extr(nzt_soil)               &
+                                           * qsws_veg_eb(j,i) * drho_l_lv  )   &
+                                     ) * ddz_soil_stag(nzt_soil)             
+
+                   m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
+                                                   + dt_3d * ( tsc(2) * tend(:)   &
+                                                   + tsc(3) * tm_soil_m(:,j,i) )   
+   
+!
+!--                Account for dry soils (find a better solution here!)
+                   DO  k = nzb_soil, nzt_soil
+                      IF ( m_soil_p(k,j,i) < 0.0_wp )  m_soil_p(k,j,i) = 0.0_wp
+                   ENDDO
+
+!
+!--                Calculate m_soil tendencies for the next Runge-Kutta step
+                   IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                      IF ( intermediate_timestep_count == 1 )  THEN
+                         DO  k = nzb_soil, nzt_soil
+                            tm_soil_m(k,j,i) = tend(k)
+                         ENDDO
+                      ELSEIF ( intermediate_timestep_count <                   &
+                               intermediate_timestep_count_max )  THEN
+                         DO  k = nzb_soil, nzt_soil
+                            tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp&
+                                     * tm_soil_m(k,j,i)
+                         ENDDO
+                      ENDIF
+                   ENDIF
+
+                ENDIF
+
              ENDIF
 
           ENDDO
        ENDDO
-
-!
-!--    Calculate surface specific humidity
-       IF ( humidity )  THEN
-          CALL calc_q_surface
-       ENDIF
-
 
     END SUBROUTINE lsm_soil_model
@@ -1656,5 +1793,71 @@
 ! Description:
 ! ------------
-!> Calculation of specific humidity of the surface layer (surface)
+!> Calculation of roughness length for open water (lakes, ocean). The
+!> parameterization follows Charnock (1955). Two different implementations
+!> are available: as in ECMWF-IFS (Beljaars 1994) or as in FLake (Subin et al.
+!> 2012)
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_z0_water_surface
+
+       USE control_parameters,                                                 &
+           ONLY: g, kappa, molecular_viscosity
+
+       IMPLICIT NONE
+
+       INTEGER :: i  !< running index
+       INTEGER :: j  !< running index
+
+       REAL(wp), PARAMETER :: alpha_ch  = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF
+!       REAL(wp), PARAMETER :: pr_number = 0.71_wp !< molecular Prandtl number in the Charnock parameterization (differs from prandtl_number)
+!       REAL(wp), PARAMETER :: sc_number = 0.66_wp !< molecular Schmidt number in the Charnock parameterization
+!       REAL(wp) :: re_0 !< near-surface roughness Reynolds number
+
+
+       DO  i = nxlg, nxrg    
+          DO  j = nysg, nyng
+             IF ( water_surface(j,i) )  THEN
+
+!
+!--             Disabled: FLake parameterization. Ideally, the Charnock 
+!--             coefficient should depend on the water depth and the fetch
+!--             length
+!                re_0 = z0(j,i) * us(j,i) / molecular_viscosity
+!        
+!                z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
+!                              alpha_ch * us(j,i) / g )
+!
+!                z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) )
+!                z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) )
+
+!
+!--              Set minimum roughness length for u* > 0.2
+!                IF ( us(j,i) > 0.2_wp )  THEN
+!                   z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) )
+!                   z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) )
+!                ENDIF
+
+!
+!--             ECMWF IFS model parameterization after Beljaars (1994). At low
+!--             wind speed, the sea surface becomes aerodynamically smooth and
+!--             the roughness scales with the viscosity. At high wind speed, the
+!--             Charnock relation is used.
+                z0(j,i) =   ( 0.11_wp * molecular_viscosity / us(j,i) )        &
+                          + ( alpha_ch * us(j,i)**2 / g )
+
+                z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i)
+                z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i)
+
+             ENDIF
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE calc_z0_water_surface
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculation of specific humidity of the skin layer (surface). It is assumend
+!> that the skin is always saturated.
 !------------------------------------------------------------------------------!
     SUBROUTINE calc_q_surface
@@ -1665,4 +1868,5 @@
        INTEGER :: j              !< running index
        INTEGER :: k              !< running index
+
        REAL(wp) :: resistance    !< aerodynamic and soil resistance term
 
@@ -1701,4 +1905,5 @@
     END SUBROUTINE calc_q_surface
 
+
 !------------------------------------------------------------------------------!
 ! Description:
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 1787)
+++ palm/trunk/SOURCE/modules.f90	(revision 1788)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-!
+! Added roughness length for moisture (z0q)
 !
 ! Former revisions:
@@ -384,5 +384,8 @@
           saswsb, saswst, shf, tnudge, td_lsa_lpt, td_lsa_q, td_sub_lpt,       &
           td_sub_q, total_2d_a, total_2d_o, ts, tswst, ug_vert, unudge, us,    &
-          usws, uswst, vnudge, vg_vert, vsws, vswst, wnudge, wsubs_vert, z0, z0h
+          usws, uswst, vnudge, vg_vert, vsws, vswst, wnudge, wsubs_vert,       &
+          z0,                                                                  & !< roughness length for momentum
+          z0h,                                                                 & !< roughness length for heat
+          z0q                                                                    !< roughness length for moisture
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
@@ -432,13 +435,14 @@
     USE kinds
 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lwp_av,                & !> Avg. liquid water path
-                                              precipitation_rate_av, & !> Avg. of precipitation rate
-                                              ol_av,                 & !> Avg. of Obukhov length
-                                              qsws_av,               & !> Avg. of surface moisture flux
-                                              shf_av,                & !> Avg. of surface heat flux
-                                              ts_av,                 & !> Avg. of characteristic temperature scale
-                                              us_av,                 & !> Avg. of friction velocity
-                                              z0_av,                 & !> Avg. of roughness length for momentum
-                                              z0h_av                   !> Avg. of roughness length for heat and moisture
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lwp_av,                & !< Avg. liquid water path
+                                              precipitation_rate_av, & !< Avg. of precipitation rate
+                                              ol_av,                 & !< Avg. of Obukhov length
+                                              qsws_av,               & !< Avg. of surface moisture flux
+                                              shf_av,                & !< Avg. of surface heat flux
+                                              ts_av,                 & !< Avg. of characteristic temperature scale
+                                              us_av,                 & !< Avg. of friction velocity
+                                              z0_av,                 & !< Avg. of roughness length for momentum
+                                              z0h_av,                & !< Avg. of roughness length for heat
+                                              z0q_av                   !< Avg. of roughness length for moisture
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
Index: palm/trunk/SOURCE/package_parin.f90
===================================================================
--- palm/trunk/SOURCE/package_parin.f90	(revision 1787)
+++ palm/trunk/SOURCE/package_parin.f90	(revision 1788)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Parameter dewfall removed.
 !
 ! Former revisions:
@@ -129,12 +129,13 @@
     USE land_surface_model_mod,                                                &
         ONLY: alpha_vangenuchten, c_surface, canopy_resistance_coefficient,    &
-              conserve_water_content, dewfall, field_capacity,                 &
+              conserve_water_content, field_capacity,                          &
               f_shortwave_incoming, hydraulic_conductivity,                    &
               lambda_surface_stable, lambda_surface_unstable, leaf_area_index, &
               land_surface, l_vangenuchten, min_canopy_resistance,             &
-              min_soil_resistance, n_vangenuchten, residual_moisture,          &
+              min_soil_resistance, n_vangenuchten, pave_heat_capacity,         &
+              pave_depth, pave_heat_conductivity, residual_moisture,           &
               root_fraction, saturation_moisture, wilting_point,               &
               skip_time_do_lsm, soil_moisture, soil_temperature, soil_type,    &
-              vegetation_coverage, veg_type, zs, z0_eb, z0h_eb
+              vegetation_coverage, veg_type, zs, z0_eb, z0h_eb, z0q_eb
 
     USE kinds
@@ -205,5 +206,5 @@
     NAMELIST /lsm_par/            alpha_vangenuchten, c_surface,               &
                                   canopy_resistance_coefficient,               &
-                                  conserve_water_content, dewfall,             &
+                                  conserve_water_content,                      &
                                   f_shortwave_incoming, field_capacity,        & 
                                   hydraulic_conductivity,                      &
@@ -212,9 +213,11 @@
                                   l_vangenuchten, min_canopy_resistance,       &
                                   min_soil_resistance, n_vangenuchten,         &
+                                  pave_depth, pave_heat_capacity,              &
+                                  pave_heat_conductivity,                      &
                                   residual_moisture, root_fraction,            &
                                   saturation_moisture, skip_time_do_lsm,       &
                                   soil_moisture, soil_temperature, soil_type,  &
                                   vegetation_coverage, veg_type, wilting_point,& 
-                                  zs, z0_eb, z0h_eb
+                                  zs, z0_eb, z0h_eb, z0q_eb
 
 
Index: palm/trunk/SOURCE/radiation_model.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model.f90	(revision 1787)
+++ palm/trunk/SOURCE/radiation_model.f90	(revision 1788)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added new albedo class for pavements / roads.
 !
 ! Former revisions:
@@ -112,4 +112,8 @@
 
 #if defined ( __rrtmg )
+
+!     USE netcdf_interface,                                                      &
+!         ONLY:  nc_stat, netcdf_handle_error
+
     USE parrrsw,                                                               &
         ONLY:  naerec, nbndsw
@@ -139,5 +143,5 @@
 !
 !-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
-    CHARACTER(37), DIMENSION(0:16), PARAMETER :: albedo_type_name = (/      &
+    CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/      &
                                    'user defined                         ', & !  0 
                                    'ocean                                ', & !  1
@@ -156,5 +160,6 @@
                                    'land ice                             ', & ! 14 
                                    'sea ice                              ', & ! 15 
-                                   'snow                                 '  & ! 16
+                                   'snow                                 ', & ! 16
+                                   'pavement/roads                       '  & ! 17
                                                          /)
 
@@ -163,9 +168,4 @@
                     day_init     = 172,  & !< day of the year at model start (21/06)
                     dots_rad     = 0       !< starting index for timeseries output 
-
-
-
-
-
 
     LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
@@ -214,5 +214,5 @@
 !-- Land surface albedos for solar zenith angle of 60Â° after Briegleb (1992)     
 !-- (shortwave, longwave, broadband):   sw,      lw,      bb,
-    REAL(wp), DIMENSION(0:2,1:16), PARAMETER :: albedo_pars = RESHAPE( (/& 
+    REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/& 
                                    0.06_wp, 0.06_wp, 0.06_wp,            & !  1
                                    0.09_wp, 0.28_wp, 0.19_wp,            & !  2
@@ -230,6 +230,7 @@
                                    0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
                                    0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
-                                   0.95_wp, 0.70_wp, 0.82_wp             & ! 16
-                                 /), (/ 3, 16 /) )
+                                   0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
+                                   0.08_wp, 0.08_wp, 0.08_wp             & ! 17
+                                 /), (/ 3, 17 /) )
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
@@ -275,4 +276,6 @@
                     rrtm_iaer = 0,     & !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
                     rrtm_idrv = 1        !< longwave upward flux calculation option (0,1)
+
+    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
 
     LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
@@ -1030,4 +1033,10 @@
                  rrtm_aldir(0,:,:) = aldif
                  rrtm_asdir(0,:,:) = asdif
+
+!
+!--        Asphalt
+           ELSEIF ( albedo_type == 17 )  THEN
+                 rrtm_aldir(0,:,:) = aldif
+                 rrtm_asdir(0,:,:) = asdif
 !
 !--        Land surfaces
Index: palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/read_3d_binary.f90	(revision 1787)
+++ palm/trunk/SOURCE/read_3d_binary.f90	(revision 1788)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added z0q and z0q_av
 ! 
 ! Former revisions:
@@ -97,5 +97,5 @@
                rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us,    &
                usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst,   &
-               w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h
+               w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
 
     USE averaging
@@ -354,5 +354,5 @@
 !--    First compare the version numbers
        READ ( 13 )  version_on_file
-       binary_version = '4.2'
+       binary_version = '4.3'
        IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
           WRITE( message_string, * ) 'version mismatch concerning data ',      &
@@ -1358,4 +1358,17 @@
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
                    z0h_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                       tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'z0q' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   z0q(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                     tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'z0q_av' )
+                   IF ( .NOT. ALLOCATED( z0q_av ) )  THEN
+                      ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   z0q_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
                        tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1787)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 1788)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added z0q and z0q_av
 ! 
 ! Former revisions:
@@ -92,5 +92,5 @@
     USE arrays_3d,                                                             &
         ONLY:  dzw, e, nr, ol, p, pt, q, qc, ql, ql_c, ql_v, qr, qsws, rho, sa,&
-               shf, ts, u, us, v, vpt, w, z0, z0h
+               shf, ts, u, us, v, vpt, w, z0, z0h, z0q
 
     USE averaging,                                                             &
@@ -98,5 +98,6 @@
                precipitation_rate_av, pt_av, q_av, qc_av, ql_av, ql_c_av,      &
                ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_av, s_av, sa_av,  &
-               shf_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av
+               shf_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av,  &
+               z0q_av
 
     USE cloud_parameters,                                                      &
@@ -494,4 +495,10 @@
                 ENDIF
                 z0h_av = 0.0_wp
+
+             CASE ( 'z0q*' )
+                IF ( .NOT. ALLOCATED( z0q_av ) )  THEN
+                   ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                z0q_av = 0.0_wp
 
              CASE DEFAULT
@@ -1006,4 +1013,11 @@
              ENDDO
 
+          CASE ( 'z0q*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   z0q_av(j,i) = z0q_av(j,i) + z0q(j,i)
+                ENDDO
+             ENDDO
+
           CASE DEFAULT
 !
Index: palm/trunk/SOURCE/surface_layer_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes.f90	(revision 1787)
+++ palm/trunk/SOURCE/surface_layer_fluxes.f90	(revision 1788)
@@ -19,5 +19,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! Added parameter z0q which replaces z0h in the similarity functions for
+! humidity.
+! Syntax layout improved.
 ! 
 ! Former revisions:
@@ -117,5 +119,5 @@
     USE arrays_3d,                                                             &
         ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
-               shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
+               shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h, z0q
 
     USE cloud_parameters,                                                      &
@@ -234,5 +236,5 @@
 !--    Use either Newton iteration or a lookup table for the bulk Richardson
 !--    number to calculate the Obukhov length 
-       ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' )  THEN
+       ELSEIF ( most_method == 'newton'  .OR.  most_method == 'lookup' )  THEN
 
           CALL calc_uv_total
@@ -302,5 +304,5 @@
 !--       Check for roughness heterogeneity. In that case terminate run and
 !--       inform user
-          IF ( MINVAL( z0h ) /= MAXVAL( z0h ) .OR.                             &
+          IF ( MINVAL( z0h ) /= MAXVAL( z0h )  .OR.                            &
                MINVAL( z0  ) /= MAXVAL( z0  ) )  THEN
              terminate_run_l = .TRUE.
@@ -347,5 +349,5 @@
 !
 !--       Calculate stretching factor for the free convection range
-          DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
+          DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
              regr_old = regr
              regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
@@ -533,5 +535,5 @@
 !--             Ensure that the bulk Richardson number and the Obukhov 
 !--             lengtH have the same sign
-                IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR.                        &
+                IF ( rib(j,i) * ol(j,i) < 0.0_wp  .OR.                         &
                      ABS( ol(j,i) ) == ol_max )  THEN
                    IF ( rib(j,i) > 0.0_wp ) ol(j,i) =  0.01_wp
@@ -638,5 +640,5 @@
 
           !$OMP PARALLEL DO PRIVATE( k, z_mo )
-          !# WARNING: does not work on GPU so far because of DO WHILE construct
+          !# WARNING: does not work on GPU so far because of DO  WHILE construct
           !!!!!!$acc kernels loop
           DO  i = nxl, nxr
@@ -656,10 +658,10 @@
                 l = l_bnd
                 IF ( rib_tab(l) - rib(j,i) > 0.0_wp )  THEN
-                   DO WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp .AND. l > 0 )
+                   DO  WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp  .AND.  l > 0 )
                       l = l-1
                    ENDDO
                 ELSE
-                   DO WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                &
-                              .AND. l < num_steps-1 )
+                   DO  WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                  &
+                              .AND.  l < num_steps-1 )
                       l = l+1
                    ENDDO
@@ -669,6 +671,6 @@
 !
 !--             Linear interpolation to find the correct value of z/L
-                ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )       &
-                            / (  rib_tab(l) - rib_tab(l-1) )                &
+                ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )          &
+                            / (  rib_tab(l) - rib_tab(l-1) )                   &
                             * ( rib(j,i) - rib_tab(l-1) ) )
 
@@ -797,5 +799,5 @@
 !
 !--       For a given surface temperature:
-          IF ( large_scale_forcing .AND. lsf_surf )  THEN
+          IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
              !$OMP PARALLEL DO
              !$acc kernels loop private( j, k )
@@ -837,5 +839,5 @@
           IF ( constant_waterflux )  THEN
 !
-!--          For a given water flux in the Prandtl layer:
+!--          For a given water flux in the surface layer
              !$OMP PARALLEL DO
              !$acc kernels loop private( j )
@@ -847,8 +849,8 @@
 
           ELSE
-             coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND.      &
+             coupled_run = ( coupling_mode == 'atmosphere_to_ocean'  .AND.     &
                              run_coupled )
 
-             IF ( large_scale_forcing .AND. lsf_surf )  THEN
+             IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
                 !$OMP PARALLEL DO
                 !$acc kernels loop private( j, k )
@@ -862,5 +864,5 @@
 
              !$OMP PARALLEL DO PRIVATE( e_s, k, z_mo )
-             !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0h ) private( e_s, j, k, z_mo )
+             !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0q ) private( e_s, j, k, z_mo )
              DO  i = nxlg, nxrg
                 !$acc loop independent
@@ -882,13 +884,13 @@
                    IF ( cloud_physics )  THEN
                       qs(j,i) = kappa * ( qv1(j,i) - q(k,j,i) )                &
-                                        / ( LOG( z_mo / z0h(j,i) )             &
+                                        / ( LOG( z_mo / z0q(j,i) )             &
                                             - psi_h( z_mo / ol(j,i) )          &
-                                            + psi_h( z0h(j,i) / ol(j,i) ) )
+                                            + psi_h( z0q(j,i) / ol(j,i) ) )
 
                    ELSE
                       qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) )              &
-                                        / ( LOG( z_mo / z0h(j,i) )             &
+                                        / ( LOG( z_mo / z0q(j,i) )             &
                                             - psi_h( z_mo / ol(j,i) )          &
-                                            + psi_h( z0h(j,i) / ol(j,i) ) )
+                                            + psi_h( z0q(j,i) / ol(j,i) ) )
                    ENDIF
 
@@ -901,8 +903,9 @@
 !
 !--    If required compute qr* and nr*
-       IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )   &
+       THEN
 
           !$OMP PARALLEL DO PRIVATE( k, z_mo )
-          !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0h ) private( j, k, z_mo )
+          !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0q ) private( j, k, z_mo )
           DO  i = nxlg, nxrg
              !$acc loop independent
@@ -913,12 +916,12 @@
 
                 qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) )              &
-                                 / ( LOG( z_mo / z0h(j,i) )                 &
+                                 / ( LOG( z_mo / z0q(j,i) )                 &
                                      - psi_h( z_mo / ol(j,i) )              &
-                                     + psi_h( z0h(j,i) / ol(j,i) ) )
+                                     + psi_h( z0q(j,i) / ol(j,i) ) )
 
                 nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) )              &
-                                 / ( LOG( z_mo / z0h(j,i) )                 &
+                                 / ( LOG( z_mo / z0q(j,i) )                 &
                                      - psi_h( z_mo / ol(j,i) )              &
-                                     + psi_h( z0h(j,i) / ol(j,i) ) )
+                                     + psi_h( z0q(j,i) / ol(j,i) ) )
              ENDDO
           ENDDO
@@ -1003,6 +1006,6 @@
 !
 !--    Compute the vertical kinematic heat flux
-       IF ( .NOT. constant_heatflux .AND. ( simulated_time <=                &
-            skip_time_do_lsm .OR. .NOT. land_surface ) )  THEN
+       IF (  .NOT.  constant_heatflux .AND.  ( simulated_time <=               &
+            skip_time_do_lsm  .OR.  .NOT.  land_surface ) )  THEN
           !$OMP PARALLEL DO
           !$acc kernels loop independent present( shf, ts, us )
@@ -1018,7 +1021,7 @@
 !
 !-- Compute the vertical water/scalar flux
-       IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar )    &
-            .AND. ( simulated_time <= skip_time_do_lsm .OR. .NOT.            &
-            land_surface ) )  THEN
+       IF (  .NOT.  constant_waterflux  .AND.  ( humidity  .OR.                &
+             passive_scalar )  .AND.  ( simulated_time <= skip_time_do_lsm     &
+            .OR.  .NOT.  land_surface ) )  THEN
           !$OMP PARALLEL DO
           !$acc kernels loop independent present( qs, qsws, us )
@@ -1089,5 +1092,5 @@
 
        IF ( zeta < 0.0_wp )  THEN
-          x = SQRT( SQRT(1.0_wp  - 16.0_wp * zeta ) )
+          x = SQRT( SQRT( 1.0_wp  - 16.0_wp * zeta ) )
           psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2    &
                   * ( 1.0_wp + x**2 ) * 0.125_wp )
@@ -1126,5 +1129,5 @@
 
        IF ( zeta < 0.0_wp )  THEN
-          x = SQRT(1.0_wp  - 16.0_wp * zeta )
+          x = SQRT( 1.0_wp  - 16.0_wp * zeta )
           psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
        ELSE
Index: palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/write_3d_binary.f90	(revision 1787)
+++ palm/trunk/SOURCE/write_3d_binary.f90	(revision 1788)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added z0q and z0q_av
 ! 
 ! Former revisions:
@@ -90,5 +90,5 @@
                rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us,    &
                usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst,   &
-               w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h
+               w, w_m_l, w_m_n, w_m_r, w_m_s, z0, z0h, z0q
         
     USE averaging
@@ -142,5 +142,5 @@
 !
 !-- Write arrays.
-    binary_version = '4.2'
+    binary_version = '4.3'
 
     WRITE ( 14 )  binary_version
@@ -470,4 +470,8 @@
        WRITE ( 14 )  'z0h_av              ';  WRITE ( 14 )  z0h_av
     ENDIF
+    WRITE ( 14 )  'z0q                 ';  WRITE ( 14 )  z0q
+    IF ( ALLOCATED( z0q_av ) )  THEN
+       WRITE ( 14 )  'z0q_av              ';  WRITE ( 14 )  z0q_av
+    ENDIF
 
 !
