Changeset 1817 for palm/trunk/SOURCE/land_surface_model_mod.f90

Timestamp:

Apr 6, 2016 3:44:20 PM (8 years ago)

Author:

maronga

Message:

further modularization of land surface model

File:

• palm/trunk/SOURCE/land_surface_model_mod.f90 (moved) (moved from palm/trunk/SOURCE/land_surface_model.f90) (14 diffs)

palm/trunk/SOURCE/land_surface_model_mod.f90

@@ -1 @@
-!> @file land_surface_model.f90
+!> @file land_surface_model_mod.f90
 !--------------------------------------------------------------------------------!
 ! This file is part of PALM.
@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added interface for lsm_init_arrays. Added subroutines for check_parameters,
+! header, and parin. Renamed some subroutines.
 ! 
 ! Former revisions:
@@ -139 +140 @@
         ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,    &
                initializing_actions, intermediate_timestep_count_max,          &
-               max_masks, precipitation, pt_surface, rho_surface,              &
-               roughness_length, surface_pressure, timestep_scheme, tsc,       &
-               z0h_factor, time_since_reference_point
+               max_masks, precipitation, pt_surface,                           &
+               rho_surface, roughness_length, surface_pressure,                &
+               timestep_scheme, tsc, z0h_factor, time_since_reference_point
 
     USE indices,                                                               &
@@ -523 +524 @@
     PRIVATE
 
-
+    
+!
+!-- Public functions
+    PUBLIC lsm_check_data_output, lsm_check_data_output_pr,                    &
+           lsm_check_parameters, lsm_energy_balance, lsm_header, lsm_init,     &
+           lsm_init_arrays, lsm_parin, lsm_soil_model 
 !
 !-- Public parameters, constants and initial values
     PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient,       &
            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, 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, z0q_eb
+           f_shortwave_incoming, hydraulic_conductivity,                       &
+           lambda_surface_stable, lambda_surface_unstable,                     &
+           land_surface, leaf_area_index,                                      &
+           lsm_swap_timelevel, l_vangenuchten,                                 &
+           min_canopy_resistance, 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, z0q_eb
 
 !
@@ -554 +560 @@
     PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
 
-    INTERFACE init_lsm
-       MODULE PROCEDURE init_lsm
-    END INTERFACE init_lsm
-
+
+    INTERFACE lsm_check_data_output
+       MODULE PROCEDURE lsm_check_data_output
+    END INTERFACE lsm_check_data_output
+    
+    INTERFACE lsm_check_data_output_pr
+       MODULE PROCEDURE lsm_check_data_output_pr
+    END INTERFACE lsm_check_data_output_pr
+    
+    INTERFACE lsm_check_parameters
+       MODULE PROCEDURE lsm_check_parameters
+    END INTERFACE lsm_check_parameters
+    
     INTERFACE lsm_energy_balance
        MODULE PROCEDURE lsm_energy_balance
     END INTERFACE lsm_energy_balance
 
+    INTERFACE lsm_header
+       MODULE PROCEDURE lsm_header
+    END INTERFACE lsm_header
+    
+    INTERFACE lsm_init
+       MODULE PROCEDURE lsm_init
+    END INTERFACE lsm_init
+
+    INTERFACE lsm_init_arrays
+       MODULE PROCEDURE lsm_init_arrays
+    END INTERFACE lsm_init_arrays
+    
+    INTERFACE lsm_parin
+       MODULE PROCEDURE lsm_parin
+    END INTERFACE lsm_parin
+    
     INTERFACE lsm_soil_model
        MODULE PROCEDURE lsm_soil_model
@@ -571 +602 @@
 
  CONTAINS
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check data output for land surface model
+!------------------------------------------------------------------------------!
+    SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
+ 
+ 
+       USE control_parameters,                                                 &
+           ONLY: data_output, message_string
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=*) ::  unit     !< 
+       CHARACTER (LEN=*) ::  var      !<
+
+       INTEGER(iwp) :: i
+       INTEGER(iwp) :: ilen   
+       INTEGER(iwp) :: k
+
+       SELECT CASE ( TRIM( var ) )
+
+          CASE ( 'm_soil' )
+             IF (  .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
+             unit = 'm3/m3'
+             
+          CASE ( 't_soil' )
+             IF (  .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
+             unit = 'K'   
+             
+          CASE ( 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*',      &
+                 'm_liq_eb*', 'qsws_eb*',      &
+                 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', &
+                 'r_a*', 'r_s*', 'shf_eb*' )
+             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+                message_string = 'illegal value for data_output: "' //         &
+                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
+                                 'cross sections are allowed for this value'
+                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
+             ENDIF
+
+             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
+                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_veg*'  .AND.  .NOT. land_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+             ENDIF
+             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 ) == 'm_liq_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_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 )   &
+             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_soil_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_veg_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 ) == 'r_a*'  .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 ) == 'r_s*'  .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_liq*' )  unit = 'none'
+             IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
+             IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
+             IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
+             IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
+             IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
+             IF ( TRIM( var ) == 'r_s*')     unit = 's/m' 
+             IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
+             
+          CASE DEFAULT
+             unit = 'illegal'
+
+       END SELECT
+
+
+    END SUBROUTINE lsm_check_data_output
+
+ SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
+ 
+       USE control_parameters,                                                 &
+           ONLY: data_output_pr, message_string
+
+       USE indices
+
+       USE profil_parameter
+
+       USE statistics
+
+       IMPLICIT NONE
+   
+       CHARACTER (LEN=*) ::  unit      !< 
+       CHARACTER (LEN=*) ::  variable  !< 
+       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
+ 
+       INTEGER(iwp) ::  user_pr_index !< 
+       INTEGER(iwp) ::  var_count     !< 
+
+       SELECT CASE ( TRIM( variable ) )
+       
+          CASE ( 't_soil', '#t_soil' )
+             IF (  .NOT.  land_surface )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not implemented for land_surface = .FALSE.'
+                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 89
+                dopr_unit     = 'K'
+                hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
+                IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
+                   dopr_initial_index(var_count) = 90
+                   hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
+                   data_output_pr(var_count)     = data_output_pr(var_count)(2:)
+                ENDIF
+                unit = dopr_unit
+             ENDIF
+
+          CASE ( 'm_soil', '#m_soil' )
+             IF (  .NOT.  land_surface )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 ' not implemented for land_surface = .FALSE.'
+                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 91
+                dopr_unit     = 'm3/m3'
+                hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
+                IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
+                   dopr_initial_index(var_count) = 92
+                   hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
+                   data_output_pr(var_count)     = data_output_pr(var_count)(2:)
+                ENDIF
+                 unit = dopr_unit
+             ENDIF
+
+
+          CASE DEFAULT
+             unit = 'illegal'
+
+       END SELECT
+
+
+    END SUBROUTINE lsm_check_data_output_pr
+ 
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check parameters routine for land surface model
+!------------------------------------------------------------------------------!
+    SUBROUTINE lsm_check_parameters
+
+       USE control_parameters,                                                 &
+           ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string,         &
+                 most_method, topography
+                 
+       USE radiation_model_mod,                                                &
+           ONLY: radiation
+    
+    
+       IMPLICIT NONE
+
+ 
+!
+!--    Dirichlet boundary conditions are required as the surface fluxes are
+!--    calculated from the temperature/humidity gradients in the land surface
+!--    model
+       IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
+          message_string = 'lsm requires setting of'//                         &
+                           'bc_pt_b = "dirichlet" and '//                      &
+                           'bc_q_b  = "dirichlet"'
+          CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF (  .NOT.  constant_flux_layer )  THEN
+          message_string = 'lsm requires '//                                   &
+                           'constant_flux_layer = .T.'
+          CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( topography /= 'flat' )  THEN
+          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
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of root_fraction(0:3)'//       &
+                              '/= 9999999.9 and SUM(root_fraction) = 1'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+ 
+          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
+             message_string = 'veg_type = 0 (user defined)'//                  &
+                              'requires setting of min_canopy_resistance'//    &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( leaf_area_index == 9999999.9_wp )  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of leaf_area_index'//          &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( vegetation_coverage == 9999999.9_wp )  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of vegetation_coverage'//      &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of'//                          &
+                              'canopy_resistance_coefficient /= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of lambda_surface_stable'//    &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of lambda_surface_unstable'//  &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of f_shortwave_incoming'//     &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( z0_eb == 9999999.9_wp )  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of z0_eb'//                    &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( z0h_eb == 9999999.9_wp )  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of z0h_eb'//                   &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
+          ENDIF
+
+
+       ENDIF
+
+       IF ( soil_type == 0 )  THEN
+
+          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of alpha_vangenuchten'//       &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( l_vangenuchten == 9999999.9_wp )  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of l_vangenuchten'//           &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( n_vangenuchten == 9999999.9_wp )  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of n_vangenuchten'//           &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of hydraulic_conductivity'//   &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( saturation_moisture == 9999999.9_wp )  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of saturation_moisture'//      &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( field_capacity == 9999999.9_wp )  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of field_capacity'//           &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( wilting_point == 9999999.9_wp )  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of wilting_point'//            &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( residual_moisture == 9999999.9_wp )  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of residual_moisture'//        &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+          ENDIF
+
+       ENDIF
+
+       IF (  .NOT.  radiation )  THEN
+          message_string = 'lsm requires '//                                   &
+                           'radiation = .T.'
+          CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
+       ENDIF
+       
+       
+    END SUBROUTINE lsm_check_parameters
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Solver for the energy balance at the surface.
+!------------------------------------------------------------------------------!
+    SUBROUTINE lsm_energy_balance
+
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i         !< running index
+       INTEGER(iwp) ::  j         !< running index
+       INTEGER(iwp) ::  k, ks     !< running index
+
+       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
+                   m_min,       & !< minimum soil moisture
+                   e,           & !< water vapour pressure
+                   e_s,         & !< water vapour saturation pressure
+                   e_s_dt,      & !< derivate of e_s with respect to T
+                   tend,        & !< tendency
+                   dq_s_dt,     & !< derivate of q_s with respect to T
+                   coef_1,      & !< coef. for prognostic equation
+                   coef_2,      & !< coef. for prognostic equation
+                   f_qsws,      & !< factor for qsws_eb
+                   f_qsws_veg,  & !< factor for qsws_veg_eb
+                   f_qsws_soil, & !< factor for qsws_soil_eb
+                   f_qsws_liq,  & !< factor for qsws_liq_eb
+                   f_shf,       & !< factor for shf_eb
+                   lambda_surface, & !< Current value of lambda_surface
+                   m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
+                   pt1,         & !< potential temperature at first grid level
+                   qv1            !< specific humidity at first grid level
+
+!
+!--    Calculate the exner function for the current time step
+       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
+
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             k = nzb_s_inner(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
+
+                c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
+                lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
+
+             ENDIF
+
+!
+!--          First step: calculate aerodyamic resistance. As pt, us, ts
+!--          are not available for the prognostic time step, data from the last
+!--          time step is used here. Note that this formulation is the
+!--          equivalent to the ECMWF formulation using drag coefficients
+             IF ( cloud_physics )  THEN
+                pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
+                qv1 = q(k+1,j,i) - ql(k+1,j,i)
+             ELSE
+                pt1 = pt(k+1,j,i)
+                qv1 = q(k+1,j,i)
+             ENDIF
+
+             r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
+
+!
+!--          Make sure that the resistance does not drop to zero
+             IF ( ABS(r_a(j,i)) < 1.0E-10_wp )  r_a(j,i) = 1.0E-10_wp
+
+!
+!--          Second step: calculate canopy resistance r_canopy
+!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
+ 
+!--          f1: correction for incoming shortwave radiation (stomata close at 
+!--          night)
+             IF ( radiation_scheme /= 'constant' )  THEN
+                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /  &
+                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)          &
+                               + 1.0_wp)) )
+             ELSE
+                f1 = 1.0_wp
+             ENDIF
+
+
+!
+!--          f2: correction for soil moisture availability to plants (the 
+!--          integrated soil moisture must thus be considered here)
+!--          f2 = 0 for very dry soils
+             m_total = 0.0_wp
+             DO  ks = nzb_soil, nzt_soil
+                 m_total = m_total + root_fr(ks,j,i)                           &
+                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
+             ENDDO 
+
+             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
+                f2 = 1.0_wp
+             ELSE
+                f2 = 1.0E-20_wp
+             ENDIF
+
+!
+!--          Calculate water vapour pressure at saturation 
+             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
+                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
+
+!
+!--          f3: correction for vapour pressure deficit
+             IF ( g_d(j,i) /= 0.0_wp )  THEN
+!
+!--             Calculate vapour pressure
+                e  = qv1 * surface_pressure / 0.622_wp
+                f3 = EXP ( -g_d(j,i) * (e_s - e) )
+             ELSE
+                f3 = 1.0_wp
+             ENDIF
+
+!
+!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
+!--          this calculation is obsolete, as r_canopy is not used below.
+!--          To do: check for very dry soil -> r_canopy goes to infinity
+             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
+                                             + 1.0E-20_wp)
+
+!
+!--          Third step: calculate bare soil resistance r_soil. The Clapp & 
+!--          Hornberger parametrization does not consider c_veg.
+             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)
+             ELSE
+                m_min = m_wilt(j,i)
+             ENDIF
+
+             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
+             f2 = MAX(f2,1.0E-20_wp)
+             f2 = MIN(f2,1.0_wp)
+
+             r_soil(j,i) = r_soil_min(j,i) / f2
+
+!
+!--          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
+
+!
+!--          Calculate saturation specific humidity
+             q_s = 0.622_wp * e_s / surface_pressure
+
+!
+!--          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
+                r_canopy(j,i) = 0.0_wp
+                r_soil(j,i)   = 0.0_wp
+             ENDIF
+
+!
+!--          Calculate coefficients for the total evapotranspiration 
+!--          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)
+             ENDIF
+!
+!--          If soil moisture is below wilting point, plants do no longer
+!--          transpirate.
+!              IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
+!                 f_qsws_veg = 0.0_wp
+!              ENDIF
+
+             f_shf  = rho_cp / r_a(j,i)
+             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
+
+!
+!--          Calculate derivative of q_s for Taylor series expansion
+             e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
+                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
+                              / (t_surface(j,i) - 35.86_wp)**2 )
+
+             dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
+
+!
+!--          Add LW up so that it can be removed in prognostic equation
+             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 )
+!
+!--             Numerator of the prognostic equation
+                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
+                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
+                         + f_shf * pt1 + f_qsws * ( qv1 - q_s                  &
+                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
+                         * t_soil(nzb_soil,j,i)
+
+!
+!--             Denominator of the prognostic equation
+                coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt           &
+                         + lambda_surface + f_shf / exn
+#else
+
+!
+!--             Numerator of the prognostic equation
+                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
+                         * t_surface(j,i) ** 4                                 &
+                         + f_shf * pt1 + f_qsws * ( qv1                        &
+                         - q_s + dq_s_dt * t_surface(j,i) )                    &
+                         + lambda_surface * t_soil(nzb_soil,j,i)
+
+!
+!--             Denominator of the prognostic equation
+                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
+                         * dq_s_dt + lambda_surface + f_shf / exn
+ 
+#endif
+             ELSE
+
+#if defined ( __rrtmg )
+!
+!--             Numerator of the prognostic equation
+                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
+                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
+                         + f_shf * pt1  + lambda_surface                       &
+                         * t_soil(nzb_soil,j,i)
+
+!
+!--             Denominator of the prognostic equation
+                coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
+#else
+
+!
+!--             Numerator of the prognostic equation
+                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
+                          * t_surface(j,i) ** 4 + f_shf * pt1                  &
+                          + lambda_surface * t_soil(nzb_soil,j,i)
+
+!
+!--             Denominator of the prognostic equation
+                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
+                         + lambda_surface + f_shf / exn
+#endif
+             ENDIF
+
+             tend = 0.0_wp
+
+!
+!--          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_tmp *    &
+                                t_surface(j,i) ) / ( c_surface_tmp + coef_2    &
+                                * dt_3d * tsc(2) ) 
+
+!
+!--          Add RK3 term
+             IF ( c_surface_tmp /= 0.0_wp )  THEN
+
+                t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
+                                   * tt_surface_m(j,i)
+
+!
+!--             Calculate true tendency
+                tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)     &
+                       * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
+!
+!--             Calculate t_surface tendencies for the next Runge-Kutta step
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( intermediate_timestep_count == 1 )  THEN
+                      tt_surface_m(j,i) = tend
+                   ELSEIF ( intermediate_timestep_count <                      &
+                            intermediate_timestep_count_max )  THEN
+                      tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp        &
+                                          * tt_surface_m(j,i)
+                   ENDIF
+                ENDIF
+             ENDIF
+
+!
+!--          In case of fast changes in the skin temperature, it is possible to
+!--          update the radiative fluxes independently from the prescribed
+!--          radiation call frequency. This effectively prevents oscillations,
+!--          especially when setting skip_time_do_radiation /= 0. The threshold
+!--          value of 0.2 used here is just a first guess. This method should be
+!--          revised in the future as tests have shown that the threshold is 
+!--          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.     &
+                  unscheduled_radiation_calls )  THEN
+                force_radiation_call_l = .TRUE.
+             ENDIF
+
+             pt(k,j,i) = t_surface_p(j,i) / exn
+
+!
+!--          Calculate fluxes
+#if defined ( __rrtmg )
+             rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)      &
+                                * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
+                                - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
+
+             IF ( rrtm_idrv == 1 )  THEN
+                rad_net(j,i) = rad_net_l(j,i)
+                rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
+                                      + rad_lw_out_change_0(j,i)               &
+                                      * ( t_surface_p(j,i) - t_surface(j,i) )
+             ENDIF
+#else
+             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb             &
+                                * t_surface(j,i)**4 - 4.0_wp * sigma_sb        &
+                                * t_surface(j,i)**3 * t_surface_p(j,i)
+#endif
+
+             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
+                              - t_soil(nzb_soil,j,i))
+
+             shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
+
+             shf(j,i) = shf_eb(j,i) / rho_cp
+
+             IF ( humidity )  THEN
+                qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt            &
+                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
+
+                qsws(j,i) = qsws_eb(j,i) / rho_lv
+
+                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                &
+                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
+                                    * t_surface_p(j,i) )
+
+                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                &
+                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
+                                    * t_surface_p(j,i) )
+
+                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                &
+                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
+                                    * t_surface_p(j,i) )
+             ENDIF
+
+!
+!--          Calculate the true surface resistance
+             IF ( qsws_eb(j,i) == 0.0_wp )  THEN
+                r_s(j,i) = 1.0E10_wp
+             ELSE
+                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)
+             ENDIF
+
+!
+!--          Calculate change in liquid water reservoir due to dew fall or 
+!--          evaporation of liquid water
+             IF ( humidity )  THEN
+!
+!--             If precipitation is activated, add rain water to qsws_liq_eb
+!--             and qsws_soil_eb according the the vegetation coverage.
+!--             precipitation_rate is given in mm.
+                IF ( precipitation )  THEN
+
+!
+!--                Add precipitation to liquid water reservoir, if possible.
+!--                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)                      &
+                                       + c_veg(j,i) * prr(k,j,i) * hyrho(k)    &
+                                       * 0.001_wp * rho_l * l_v
+                   ELSE
+                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
+                                        + c_veg(j,i) * prr(k,j,i) * hyrho(k)   &
+                                        * 0.001_wp * rho_l * l_v
+                   ENDIF
+
+!--                Add precipitation to bare soil according to the bare soil
+!--                coverage.
+                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp             &
+                                       - c_veg(j,i)) * prr(k,j,i) * hyrho(k)   &
+                                       * 0.001_wp * rho_l * l_v
+                ENDIF
+
+!
+!--             If the air is saturated, check the reservoir water level
+                IF ( qsws_eb(j,i) < 0.0_wp )  THEN
+
+!
+!--                Check if reservoir is full (avoid values > m_liq_eb_max)
+!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this 
+!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
+!--                so that tend is zero and no further check is needed
+                   IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
+                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
+                                           + qsws_liq_eb(j,i)
+                      qsws_liq_eb(j,i)  = 0.0_wp
+                   ENDIF
+
+!
+!--                In case qsws_veg_eb becomes negative (unphysical behavior), 
+!--                let the water enter the liquid water reservoir as dew on the
+!--                plant
+                   IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
+                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
+                      qsws_veg_eb(j,i) = 0.0_wp
+                   ENDIF
+                ENDIF                    
+  
+                tend = - qsws_liq_eb(j,i) * drho_l_lv
+
+                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
+                                                   + tsc(3) * tm_liq_eb_m(j,i) )
+
+!
+!--             Check if reservoir is overfull -> reduce to maximum
+!--             (conservation of water is violated here)
+                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
+
+!
+!--             Check if reservoir is empty (avoid values < 0.0)
+!--             (conservation of water is violated here)
+                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
+
+
+!
+!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( intermediate_timestep_count == 1 )  THEN
+                      tm_liq_eb_m(j,i) = tend
+                   ELSEIF ( intermediate_timestep_count <                      &
+                            intermediate_timestep_count_max )  THEN
+                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
+                                      * tm_liq_eb_m(j,i)
+                   ENDIF
+                ENDIF
+
+             ENDIF
+
+          ENDDO
+       ENDDO
+
+!
+!--    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    &
+            == intermediate_timestep_count_max-1 )  THEN
+#if defined( __parallel )
+          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+          CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
+                              1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
+#else
+          force_radiation_call = force_radiation_call_l
+#endif
+          force_radiation_call_l = .FALSE.
+       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
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Header output for land surface model
+!------------------------------------------------------------------------------!
+    SUBROUTINE lsm_header ( io )
+
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=86) ::  t_soil_chr          !< String for soil temperature profile
+       CHARACTER (LEN=86) ::  roots_chr           !< String for root profile
+       CHARACTER (LEN=86) ::  vertical_index_chr  !< String for the vertical index
+       CHARACTER (LEN=86) ::  m_soil_chr          !< String for soil moisture
+       CHARACTER (LEN=86) ::  soil_depth_chr      !< String for soil depth
+       CHARACTER (LEN=10) ::  coor_chr            !< Temporary string
+    
+       INTEGER(iwp) ::  i                         !< Loop index over soil layers
+ 
+       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
+ 
+       t_soil_chr = ''
+       m_soil_chr    = ''
+       soil_depth_chr  = '' 
+       roots_chr        = '' 
+       vertical_index_chr   = ''
+
+       i = 1
+       DO i = nzb_soil, nzt_soil
+          WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i)
+          t_soil_chr = TRIM( t_soil_chr ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i)
+          m_soil_chr = TRIM( m_soil_chr ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F10.2,7X)')  - zs(i)
+          soil_depth_chr = TRIM( soil_depth_chr ) // ' '  // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F10.2,7X)')  root_fraction(i)
+          roots_chr = TRIM( roots_chr ) // ' '  // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(I10,7X)')  i
+          vertical_index_chr = TRIM( vertical_index_chr ) // ' '  //           &
+                               TRIM( coor_chr )
+       ENDDO
+
+!
+!--    Write land surface model header
+       WRITE( io,  1 )
+       IF ( conserve_water_content )  THEN
+          WRITE( io, 2 )
+       ELSE
+          WRITE( io, 3 )
+       ENDIF
+
+       WRITE( io, 4 ) TRIM( veg_type_name(veg_type) ),                         &
+                        TRIM (soil_type_name(soil_type) )
+       WRITE( io, 5 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
+                        TRIM( m_soil_chr ), TRIM( roots_chr ),                 &
+                        TRIM( vertical_index_chr )
+
+1   FORMAT (//' Land surface model information:'/                              &
+              ' ------------------------------'/)
+2   FORMAT (' --> Soil bottom is closed (water content is conserved, default)')
+3   FORMAT (' --> Soil bottom is open (water content is not conserved)')         
+4   FORMAT (' --> Land surface type  : ',A,/                                   &
+            ' --> Soil porosity type : ',A)
+5   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
+            '       Height:        ',A,'  m'/                                  &
+            '       Temperature:   ',A,'  K'/                                  &
+            '       Moisture:      ',A,'  m**3/m**3'/                          &
+            '       Root fraction: ',A,'  '/                                   &
+            '       Gridpoint:     ',A)
+
+
+
+    END SUBROUTINE lsm_header
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initialization of the land surface model
+!------------------------------------------------------------------------------!
+    SUBROUTINE lsm_init
+    
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i !< running index
+       INTEGER(iwp) ::  j !< running index
+       INTEGER(iwp) ::  k !< running index
+
+       REAL(wp) :: pt1   !< potential temperature at first grid level
+
+
+!
+!--    Calculate Exner function
+       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
+
+
+!
+!--    If no cloud physics is used, rho_surface has not been calculated before
+       IF (  .NOT.  cloud_physics )  THEN
+          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
+       ENDIF
+
+!
+!--    Calculate frequently used parameters
+       rho_cp    = cp * rho_surface
+       rd_d_rv   = r_d / r_v
+       rho_lv    = rho_surface * l_v
+       drho_l_lv = 1.0_wp / (rho_l * l_v)
+
+!
+!--    Set inital values for prognostic quantities
+       tt_surface_m = 0.0_wp
+       tt_soil_m    = 0.0_wp
+       tm_soil_m    = 0.0_wp
+       tm_liq_eb_m  = 0.0_wp
+       c_liq        = 0.0_wp
+
+       ghf_eb = 0.0_wp
+       shf_eb = rho_cp * shf
+
+       IF ( humidity )  THEN
+          qsws_eb = rho_l * l_v * qsws
+       ELSE
+          qsws_eb = 0.0_wp
+       ENDIF
+
+       qsws_liq_eb  = 0.0_wp
+       qsws_soil_eb = 0.0_wp
+       qsws_veg_eb  = 0.0_wp
+
+       r_a        = 50.0_wp
+       r_s        = 50.0_wp
+       r_canopy   = 0.0_wp
+       r_soil     = 0.0_wp
+
+!
+!--    Allocate 3D soil model arrays
+       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+
+       lambda_h = 0.0_wp
+!
+!--    If required, allocate humidity-related variables for the soil model
+       IF ( humidity )  THEN
+          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
+
+          lambda_w = 0.0_wp 
+       ENDIF
+
+!
+!--    Calculate grid spacings. Temperature and moisture are defined at
+!--    the edges of the soil layers (_stag), whereas gradients/fluxes are defined
+!--    at the centers
+       dz_soil(nzb_soil) = zs(nzb_soil)
+
+       DO  k = nzb_soil+1, nzt_soil
+          dz_soil(k) = zs(k) - zs(k-1)
+       ENDDO
+       dz_soil(nzt_soil+1) = dz_soil(nzt_soil)
+
+       DO  k = nzb_soil, nzt_soil-1
+          dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k))
+       ENDDO
+       dz_soil_stag(nzt_soil) = dz_soil(nzt_soil)
+
+       ddz_soil      = 1.0_wp / dz_soil
+       ddz_soil_stag = 1.0_wp / dz_soil_stag
+
+!
+!--    Initialize standard soil types. It is possible to overwrite each
+!--    parameter by setting the respecticy NAMELIST variable to a 
+!--    value /= 9999999.9. 
+       IF ( soil_type /= 0 )  THEN  
+ 
+          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
+             alpha_vangenuchten = soil_pars(0,soil_type)
+          ENDIF
+
+          IF ( l_vangenuchten == 9999999.9_wp )  THEN
+             l_vangenuchten = soil_pars(1,soil_type)
+          ENDIF
+
+          IF ( n_vangenuchten == 9999999.9_wp )  THEN
+             n_vangenuchten = soil_pars(2,soil_type)            
+          ENDIF
+
+          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
+             hydraulic_conductivity = soil_pars(3,soil_type)            
+          ENDIF
+
+          IF ( saturation_moisture == 9999999.9_wp )  THEN
+             saturation_moisture = m_soil_pars(0,soil_type)           
+          ENDIF
+
+          IF ( field_capacity == 9999999.9_wp )  THEN
+             field_capacity = m_soil_pars(1,soil_type)           
+          ENDIF
+
+          IF ( wilting_point == 9999999.9_wp )  THEN
+             wilting_point = m_soil_pars(2,soil_type)            
+          ENDIF
+
+          IF ( residual_moisture == 9999999.9_wp )  THEN
+             residual_moisture = m_soil_pars(3,soil_type)       
+          ENDIF
+
+       ENDIF    
+
+!
+!--    Map values to the respective 2D arrays
+       alpha_vg      = alpha_vangenuchten
+       l_vg          = l_vangenuchten
+       n_vg          = n_vangenuchten 
+       gamma_w_sat   = hydraulic_conductivity
+       m_sat         = saturation_moisture
+       m_fc          = field_capacity
+       m_wilt        = wilting_point
+       m_res         = residual_moisture
+       r_soil_min    = min_soil_resistance
+
+!
+!--    Initial run actions
+       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+
+          t_soil    = 0.0_wp
+          m_liq_eb  = 0.0_wp
+          m_soil    = 0.0_wp
+
+!
+!--       Map user settings of T and q for each soil layer
+!--       (make sure that the soil moisture does not drop below the permanent
+!--       wilting point) -> problems with devision by zero)
+          DO  k = nzb_soil, nzt_soil
+             t_soil(k,:,:)    = soil_temperature(k)
+             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
+             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
+          ENDDO
+          t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
+
+!
+!--       Calculate surface temperature
+          t_surface   = pt_surface * exn
+
+!
+!--       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
+                k = nzb_s_inner(j,i)
+
+                IF ( cloud_physics )  THEN
+                   pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
+                ELSE
+                   pt1 = pt(k+1,j,i)
+                ENDIF
+
+!
+!--             Assure that r_a cannot be zero at model start
+                IF ( pt1 == pt(k,j,i) )  pt1 = pt1 + 1.0E-10_wp
+
+                us(j,i)  = 0.1_wp
+                ts(j,i)  = (pt1 - pt(k,j,i)) / r_a(j,i)
+                shf(j,i) = - us(j,i) * ts(j,i)
+             ENDDO
+          ENDDO
+
+!
+!--    Actions for restart runs
+       ELSE
+
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                k = nzb_s_inner(j,i)                
+                t_surface(j,i) = pt(k,j,i) * exn
+             ENDDO
+          ENDDO
+
+       ENDIF
+
+       DO  k = nzb_soil, nzt_soil
+          root_fr(k,:,:) = root_fraction(k)
+       ENDDO
+
+       IF ( veg_type /= 0 )  THEN
+          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
+             min_canopy_resistance = veg_pars(0,veg_type)
+          ENDIF
+          IF ( leaf_area_index == 9999999.9_wp )  THEN
+             leaf_area_index = veg_pars(1,veg_type)         
+          ENDIF
+          IF ( vegetation_coverage == 9999999.9_wp )  THEN
+             vegetation_coverage = veg_pars(2,veg_type)     
+          ENDIF
+          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
+              canopy_resistance_coefficient= veg_pars(3,veg_type)     
+          ENDIF
+          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
+             lambda_surface_stable = surface_pars(0,veg_type)         
+          ENDIF
+          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
+             lambda_surface_unstable = surface_pars(1,veg_type)        
+          ENDIF
+          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
+             f_shortwave_incoming = surface_pars(2,veg_type)        
+          ENDIF
+          IF ( z0_eb == 9999999.9_wp )  THEN
+             roughness_length = roughness_par(0,veg_type) 
+             z0_eb            = roughness_par(0,veg_type) 
+          ENDIF
+          IF ( z0h_eb == 9999999.9_wp )  THEN
+             z0h_eb = roughness_par(1,veg_type)
+          ENDIF
+          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
+             DO  k = nzb_soil, nzt_soil
+                root_fr(k,:,:) = root_distribution(k,veg_type)
+                root_fraction(k) = root_distribution(k,veg_type)
+             ENDDO
+          ENDIF
+
+       ELSE
+
+          IF ( z0_eb == 9999999.9_wp )  THEN
+             z0_eb = roughness_length
+          ENDIF
+          IF ( z0h_eb == 9999999.9_wp )  THEN
+             z0h_eb = z0_eb * z0h_factor
+          ENDIF
+          IF ( z0q_eb == 9999999.9_wp )  THEN
+             z0q_eb = z0_eb * z0h_factor
+          ENDIF
+
+       ENDIF
+
+!
+!--    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
+       c_veg                = vegetation_coverage
+       g_d                  = canopy_resistance_coefficient
+       lambda_surface_s     = lambda_surface_stable
+       lambda_surface_u     = lambda_surface_unstable
+       f_sw_in              = f_shortwave_incoming
+       z0                   = z0_eb
+       z0h                  = z0h_eb
+       z0q                  = z0q_eb
+
+!
+!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
+       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
+       m_soil_p    = m_soil
+       m_liq_eb_p  = m_liq_eb
+       t_surface_p = t_surface
+
+
+
+!--    Store initial profiles of t_soil and m_soil (assuming they are 
+!--    horizontally homogeneous on this PE)
+       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
+                                                nysg,nxlg), 2,                 &
+                                                statistic_regions+1 )
+       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
+                                                nysg,nxlg), 2,                 &
+                                                statistic_regions+1 )
+
+    END SUBROUTINE lsm_init
 
 
@@ -578 +1920 @@
 !> Allocate land surface model arrays and define pointers
 !------------------------------------------------------------------------------!
-    SUBROUTINE init_lsm_arrays
+    SUBROUTINE lsm_init_arrays
     
 
@@ -658 +2000 @@
 
 
-    END SUBROUTINE init_lsm_arrays
+    END SUBROUTINE lsm_init_arrays
+
 
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Initialization of the land surface model
+!> Parin for &lsmpar for land surface model
 !------------------------------------------------------------------------------!
-    SUBROUTINE init_lsm
-    
+    SUBROUTINE lsm_parin
+
 
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i !< running index
-       INTEGER(iwp) ::  j !< running index
-       INTEGER(iwp) ::  k !< running index
-
-       REAL(wp) :: pt1   !< potential temperature at first grid level
-
-
-!
-!--    Calculate Exner function
-       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
-
-
-!
-!--    If no cloud physics is used, rho_surface has not been calculated before
-       IF (  .NOT.  cloud_physics )  THEN
-          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
-       ENDIF
-
-!
-!--    Calculate frequently used parameters
-       rho_cp    = cp * rho_surface
-       rd_d_rv   = r_d / r_v
-       rho_lv    = rho_surface * l_v
-       drho_l_lv = 1.0_wp / (rho_l * l_v)
-
-!
-!--    Set inital values for prognostic quantities
-       tt_surface_m = 0.0_wp
-       tt_soil_m    = 0.0_wp
-       tm_soil_m    = 0.0_wp
-       tm_liq_eb_m  = 0.0_wp
-       c_liq        = 0.0_wp
-
-       ghf_eb = 0.0_wp
-       shf_eb = rho_cp * shf
-
-       IF ( humidity )  THEN
-          qsws_eb = rho_l * l_v * qsws
-       ELSE
-          qsws_eb = 0.0_wp
-       ENDIF
-
-       qsws_liq_eb  = 0.0_wp
-       qsws_soil_eb = 0.0_wp
-       qsws_veg_eb  = 0.0_wp
-
-       r_a        = 50.0_wp
-       r_s        = 50.0_wp
-       r_canopy   = 0.0_wp
-       r_soil     = 0.0_wp
-
-!
-!--    Allocate 3D soil model arrays
-       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-
-       lambda_h = 0.0_wp
-!
-!--    If required, allocate humidity-related variables for the soil model
-       IF ( humidity )  THEN
-          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
-          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
-
-          lambda_w = 0.0_wp 
-       ENDIF
-
-!
-!--    Calculate grid spacings. Temperature and moisture are defined at
-!--    the edges of the soil layers (_stag), whereas gradients/fluxes are defined
-!--    at the centers
-       dz_soil(nzb_soil) = zs(nzb_soil)
-
-       DO  k = nzb_soil+1, nzt_soil
-          dz_soil(k) = zs(k) - zs(k-1)
+       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
+       
+       NAMELIST /lsm_par/         alpha_vangenuchten, c_surface,               &
+                                  canopy_resistance_coefficient,               &
+                                  conserve_water_content,                      &
+                                  f_shortwave_incoming, field_capacity,        & 
+                                  hydraulic_conductivity,                      &
+                                  lambda_surface_stable,                       &
+                                  lambda_surface_unstable, leaf_area_index,    &
+                                  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, z0q_eb
+       
+       line = ' '
+       
+!
+!--    Try to find land surface model package
+       REWIND ( 11 )
+       line = ' '
+       DO   WHILE ( INDEX( line, '&lsm_par' ) == 0 )
+          READ ( 11, '(A)', END=10 )  line
        ENDDO
-       dz_soil(nzt_soil+1) = dz_soil(nzt_soil)
-
-       DO  k = nzb_soil, nzt_soil-1
-          dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k))
-       ENDDO
-       dz_soil_stag(nzt_soil) = dz_soil(nzt_soil)
-
-       ddz_soil      = 1.0_wp / dz_soil
-       ddz_soil_stag = 1.0_wp / dz_soil_stag
-
-!
-!--    Initialize standard soil types. It is possible to overwrite each
-!--    parameter by setting the respecticy NAMELIST variable to a 
-!--    value /= 9999999.9. 
-       IF ( soil_type /= 0 )  THEN  
- 
-          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
-             alpha_vangenuchten = soil_pars(0,soil_type)
-          ENDIF
-
-          IF ( l_vangenuchten == 9999999.9_wp )  THEN
-             l_vangenuchten = soil_pars(1,soil_type)
-          ENDIF
-
-          IF ( n_vangenuchten == 9999999.9_wp )  THEN
-             n_vangenuchten = soil_pars(2,soil_type)            
-          ENDIF
-
-          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
-             hydraulic_conductivity = soil_pars(3,soil_type)            
-          ENDIF
-
-          IF ( saturation_moisture == 9999999.9_wp )  THEN
-             saturation_moisture = m_soil_pars(0,soil_type)           
-          ENDIF
-
-          IF ( field_capacity == 9999999.9_wp )  THEN
-             field_capacity = m_soil_pars(1,soil_type)           
-          ENDIF
-
-          IF ( wilting_point == 9999999.9_wp )  THEN
-             wilting_point = m_soil_pars(2,soil_type)            
-          ENDIF
-
-          IF ( residual_moisture == 9999999.9_wp )  THEN
-             residual_moisture = m_soil_pars(3,soil_type)       
-          ENDIF
-
-       ENDIF    
-
-!
-!--    Map values to the respective 2D arrays
-       alpha_vg      = alpha_vangenuchten
-       l_vg          = l_vangenuchten
-       n_vg          = n_vangenuchten 
-       gamma_w_sat   = hydraulic_conductivity
-       m_sat         = saturation_moisture
-       m_fc          = field_capacity
-       m_wilt        = wilting_point
-       m_res         = residual_moisture
-       r_soil_min    = min_soil_resistance
-
-!
-!--    Initial run actions
-       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
-
-          t_soil    = 0.0_wp
-          m_liq_eb  = 0.0_wp
-          m_soil    = 0.0_wp
-
-!
-!--       Map user settings of T and q for each soil layer
-!--       (make sure that the soil moisture does not drop below the permanent
-!--       wilting point) -> problems with devision by zero)
-          DO  k = nzb_soil, nzt_soil
-             t_soil(k,:,:)    = soil_temperature(k)
-             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
-             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
-          ENDDO
-          t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
-
-!
-!--       Calculate surface temperature
-          t_surface   = pt_surface * exn
-
-!
-!--       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
-                k = nzb_s_inner(j,i)
-
-                IF ( cloud_physics )  THEN
-                   pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
-                ELSE
-                   pt1 = pt(k+1,j,i)
-                ENDIF
-
-!
-!--             Assure that r_a cannot be zero at model start
-                IF ( pt1 == pt(k,j,i) )  pt1 = pt1 + 1.0E-10_wp
-
-                us(j,i)  = 0.1_wp
-                ts(j,i)  = (pt1 - pt(k,j,i)) / r_a(j,i)
-                shf(j,i) = - us(j,i) * ts(j,i)
-             ENDDO
-          ENDDO
-
-!
-!--    Actions for restart runs
-       ELSE
-
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                k = nzb_s_inner(j,i)                
-                t_surface(j,i) = pt(k,j,i) * exn
-             ENDDO
-          ENDDO
-
-       ENDIF
-
-       DO  k = nzb_soil, nzt_soil
-          root_fr(k,:,:) = root_fraction(k)
-       ENDDO
-
-       IF ( veg_type /= 0 )  THEN
-          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
-             min_canopy_resistance = veg_pars(0,veg_type)
-          ENDIF
-          IF ( leaf_area_index == 9999999.9_wp )  THEN
-             leaf_area_index = veg_pars(1,veg_type)         
-          ENDIF
-          IF ( vegetation_coverage == 9999999.9_wp )  THEN
-             vegetation_coverage = veg_pars(2,veg_type)     
-          ENDIF
-          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
-              canopy_resistance_coefficient= veg_pars(3,veg_type)     
-          ENDIF
-          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
-             lambda_surface_stable = surface_pars(0,veg_type)         
-          ENDIF
-          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
-             lambda_surface_unstable = surface_pars(1,veg_type)        
-          ENDIF
-          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
-             f_shortwave_incoming = surface_pars(2,veg_type)        
-          ENDIF
-          IF ( z0_eb == 9999999.9_wp )  THEN
-             roughness_length = roughness_par(0,veg_type) 
-             z0_eb            = roughness_par(0,veg_type) 
-          ENDIF
-          IF ( z0h_eb == 9999999.9_wp )  THEN
-             z0h_eb = roughness_par(1,veg_type)
-          ENDIF
-          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
-             DO  k = nzb_soil, nzt_soil
-                root_fr(k,:,:) = root_distribution(k,veg_type)
-                root_fraction(k) = root_distribution(k,veg_type)
-             ENDDO
-          ENDIF
-
-       ELSE
-
-          IF ( z0_eb == 9999999.9_wp )  THEN
-             z0_eb = roughness_length
-          ENDIF
-          IF ( z0h_eb == 9999999.9_wp )  THEN
-             z0h_eb = z0_eb * z0h_factor
-          ENDIF
-          IF ( z0q_eb == 9999999.9_wp )  THEN
-             z0q_eb = z0_eb * z0h_factor
-          ENDIF
-
-       ENDIF
-
-!
-!--    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
-       c_veg                = vegetation_coverage
-       g_d                  = canopy_resistance_coefficient
-       lambda_surface_s     = lambda_surface_stable
-       lambda_surface_u     = lambda_surface_unstable
-       f_sw_in              = f_shortwave_incoming
-       z0                   = z0_eb
-       z0h                  = z0h_eb
-       z0q                  = z0q_eb
-
-!
-!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
-       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
-       m_soil_p    = m_soil
-       m_liq_eb_p  = m_liq_eb
-       t_surface_p = t_surface
-
-
-
-!--    Store initial profiles of t_soil and m_soil (assuming they are 
-!--    horizontally homogeneous on this PE)
-       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
-                                                nysg,nxlg), 2,                 &
-                                                statistic_regions+1 )
-       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
-                                                nysg,nxlg), 2,                 &
-                                                statistic_regions+1 )
-
-    END SUBROUTINE init_lsm
-
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Solver for the energy balance at the surface.
-!------------------------------------------------------------------------------!
-    SUBROUTINE lsm_energy_balance
-
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) ::  i         !< running index
-       INTEGER(iwp) ::  j         !< running index
-       INTEGER(iwp) ::  k, ks     !< running index
-
-       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
-                   m_min,       & !< minimum soil moisture
-                   e,           & !< water vapour pressure
-                   e_s,         & !< water vapour saturation pressure
-                   e_s_dt,      & !< derivate of e_s with respect to T
-                   tend,        & !< tendency
-                   dq_s_dt,     & !< derivate of q_s with respect to T
-                   coef_1,      & !< coef. for prognostic equation
-                   coef_2,      & !< coef. for prognostic equation
-                   f_qsws,      & !< factor for qsws_eb
-                   f_qsws_veg,  & !< factor for qsws_veg_eb
-                   f_qsws_soil, & !< factor for qsws_soil_eb
-                   f_qsws_liq,  & !< factor for qsws_liq_eb
-                   f_shf,       & !< factor for shf_eb
-                   lambda_surface, & !< Current value of lambda_surface
-                   m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
-                   pt1,         & !< potential temperature at first grid level
-                   qv1            !< specific humidity at first grid level
-
-!
-!--    Calculate the exner function for the current time step
-       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
-
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             k = nzb_s_inner(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
-
-                c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
-                lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
-
-             ENDIF
-
-!
-!--          First step: calculate aerodyamic resistance. As pt, us, ts
-!--          are not available for the prognostic time step, data from the last
-!--          time step is used here. Note that this formulation is the
-!--          equivalent to the ECMWF formulation using drag coefficients
-             IF ( cloud_physics )  THEN
-                pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
-                qv1 = q(k+1,j,i) - ql(k+1,j,i)
-             ELSE
-                pt1 = pt(k+1,j,i)
-                qv1 = q(k+1,j,i)
-             ENDIF
-
-             r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
-
-!
-!--          Make sure that the resistance does not drop to zero
-             IF ( ABS(r_a(j,i)) < 1.0E-10_wp )  r_a(j,i) = 1.0E-10_wp
-
-!
-!--          Second step: calculate canopy resistance r_canopy
-!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
- 
-!--          f1: correction for incoming shortwave radiation (stomata close at 
-!--          night)
-             IF ( radiation_scheme /= 'constant' )  THEN
-                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /  &
-                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)          &
-                               + 1.0_wp)) )
-             ELSE
-                f1 = 1.0_wp
-             ENDIF
-
-
-!
-!--          f2: correction for soil moisture availability to plants (the 
-!--          integrated soil moisture must thus be considered here)
-!--          f2 = 0 for very dry soils
-             m_total = 0.0_wp
-             DO  ks = nzb_soil, nzt_soil
-                 m_total = m_total + root_fr(ks,j,i)                           &
-                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
-             ENDDO 
-
-             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
-                f2 = 1.0_wp
-             ELSE
-                f2 = 1.0E-20_wp
-             ENDIF
-
-!
-!--          Calculate water vapour pressure at saturation 
-             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
-                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
-
-!
-!--          f3: correction for vapour pressure deficit
-             IF ( g_d(j,i) /= 0.0_wp )  THEN
-!
-!--             Calculate vapour pressure
-                e  = qv1 * surface_pressure / 0.622_wp
-                f3 = EXP ( -g_d(j,i) * (e_s - e) )
-             ELSE
-                f3 = 1.0_wp
-             ENDIF
-
-!
-!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
-!--          this calculation is obsolete, as r_canopy is not used below.
-!--          To do: check for very dry soil -> r_canopy goes to infinity
-             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
-                                             + 1.0E-20_wp)
-
-!
-!--          Third step: calculate bare soil resistance r_soil. The Clapp & 
-!--          Hornberger parametrization does not consider c_veg.
-             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)
-             ELSE
-                m_min = m_wilt(j,i)
-             ENDIF
-
-             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
-             f2 = MAX(f2,1.0E-20_wp)
-             f2 = MIN(f2,1.0_wp)
-
-             r_soil(j,i) = r_soil_min(j,i) / f2
-
-!
-!--          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
-
-!
-!--          Calculate saturation specific humidity
-             q_s = 0.622_wp * e_s / surface_pressure
-
-!
-!--          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
-                r_canopy(j,i) = 0.0_wp
-                r_soil(j,i)   = 0.0_wp
-             ENDIF
-
-!
-!--          Calculate coefficients for the total evapotranspiration 
-!--          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)
-             ENDIF
-!
-!--          If soil moisture is below wilting point, plants do no longer
-!--          transpirate.
-!              IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
-!                 f_qsws_veg = 0.0_wp
-!              ENDIF
-
-             f_shf  = rho_cp / r_a(j,i)
-             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
-
-!
-!--          Calculate derivative of q_s for Taylor series expansion
-             e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
-                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
-                              / (t_surface(j,i) - 35.86_wp)**2 )
-
-             dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
-
-!
-!--          Add LW up so that it can be removed in prognostic equation
-             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 )
-!
-!--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
-                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
-                         + f_shf * pt1 + f_qsws * ( qv1 - q_s                  &
-                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
-                         * t_soil(nzb_soil,j,i)
-
-!
-!--             Denominator of the prognostic equation
-                coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt           &
-                         + lambda_surface + f_shf / exn
-#else
-
-!
-!--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
-                         * t_surface(j,i) ** 4                                 &
-                         + f_shf * pt1 + f_qsws * ( qv1                        &
-                         - q_s + dq_s_dt * t_surface(j,i) )                    &
-                         + lambda_surface * t_soil(nzb_soil,j,i)
-
-!
-!--             Denominator of the prognostic equation
-                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
-                         * dq_s_dt + lambda_surface + f_shf / exn
- 
-#endif
-             ELSE
-
-#if defined ( __rrtmg )
-!
-!--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
-                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
-                         + f_shf * pt1  + lambda_surface                       &
-                         * t_soil(nzb_soil,j,i)
-
-!
-!--             Denominator of the prognostic equation
-                coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
-#else
-
-!
-!--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
-                          * t_surface(j,i) ** 4 + f_shf * pt1                  &
-                          + lambda_surface * t_soil(nzb_soil,j,i)
-
-!
-!--             Denominator of the prognostic equation
-                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
-                         + lambda_surface + f_shf / exn
-#endif
-             ENDIF
-
-             tend = 0.0_wp
-
-!
-!--          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_tmp *    &
-                                t_surface(j,i) ) / ( c_surface_tmp + coef_2    &
-                                * dt_3d * tsc(2) ) 
-
-!
-!--          Add RK3 term
-             IF ( c_surface_tmp /= 0.0_wp )  THEN
-
-                t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
-                                   * tt_surface_m(j,i)
-
-!
-!--             Calculate true tendency
-                tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)     &
-                       * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
-!
-!--             Calculate t_surface tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      tt_surface_m(j,i) = tend
-                   ELSEIF ( intermediate_timestep_count <                      &
-                            intermediate_timestep_count_max )  THEN
-                      tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp        &
-                                          * tt_surface_m(j,i)
-                   ENDIF
-                ENDIF
-             ENDIF
-
-!
-!--          In case of fast changes in the skin temperature, it is possible to
-!--          update the radiative fluxes independently from the prescribed
-!--          radiation call frequency. This effectively prevents oscillations,
-!--          especially when setting skip_time_do_radiation /= 0. The threshold
-!--          value of 0.2 used here is just a first guess. This method should be
-!--          revised in the future as tests have shown that the threshold is 
-!--          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.     &
-                  unscheduled_radiation_calls )  THEN
-                force_radiation_call_l = .TRUE.
-             ENDIF
-
-             pt(k,j,i) = t_surface_p(j,i) / exn
-
-!
-!--          Calculate fluxes
-#if defined ( __rrtmg )
-             rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)      &
-                                * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
-                                - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
-
-             IF ( rrtm_idrv == 1 )  THEN
-                rad_net(j,i) = rad_net_l(j,i)
-                rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
-                                      + rad_lw_out_change_0(j,i)               &
-                                      * ( t_surface_p(j,i) - t_surface(j,i) )
-             ENDIF
-#else
-             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb             &
-                                * t_surface(j,i)**4 - 4.0_wp * sigma_sb        &
-                                * t_surface(j,i)**3 * t_surface_p(j,i)
-#endif
-
-             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
-                              - t_soil(nzb_soil,j,i))
-
-             shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
-
-             shf(j,i) = shf_eb(j,i) / rho_cp
-
-             IF ( humidity )  THEN
-                qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt            &
-                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
-
-                qsws(j,i) = qsws_eb(j,i) / rho_lv
-
-                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                &
-                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
-                                    * t_surface_p(j,i) )
-
-                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                &
-                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
-                                    * t_surface_p(j,i) )
-
-                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                &
-                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
-                                    * t_surface_p(j,i) )
-             ENDIF
-
-!
-!--          Calculate the true surface resistance
-             IF ( qsws_eb(j,i) == 0.0_wp )  THEN
-                r_s(j,i) = 1.0E10_wp
-             ELSE
-                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)
-             ENDIF
-
-!
-!--          Calculate change in liquid water reservoir due to dew fall or 
-!--          evaporation of liquid water
-             IF ( humidity )  THEN
-!
-!--             If precipitation is activated, add rain water to qsws_liq_eb
-!--             and qsws_soil_eb according the the vegetation coverage.
-!--             precipitation_rate is given in mm.
-                IF ( precipitation )  THEN
-
-!
-!--                Add precipitation to liquid water reservoir, if possible.
-!--                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)                      &
-                                       + c_veg(j,i) * prr(k,j,i) * hyrho(k)    &
-                                       * 0.001_wp * rho_l * l_v
-                   ELSE
-                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
-                                        + c_veg(j,i) * prr(k,j,i) * hyrho(k)   &
-                                        * 0.001_wp * rho_l * l_v
-                   ENDIF
-
-!--                Add precipitation to bare soil according to the bare soil
-!--                coverage.
-                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp             &
-                                       - c_veg(j,i)) * prr(k,j,i) * hyrho(k)   &
-                                       * 0.001_wp * rho_l * l_v
-                ENDIF
-
-!
-!--             If the air is saturated, check the reservoir water level
-                IF ( qsws_eb(j,i) < 0.0_wp )  THEN
-
-!
-!--                Check if reservoir is full (avoid values > m_liq_eb_max)
-!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this 
-!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
-!--                so that tend is zero and no further check is needed
-                   IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
-                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
-                                           + qsws_liq_eb(j,i)
-                      qsws_liq_eb(j,i)  = 0.0_wp
-                   ENDIF
-
-!
-!--                In case qsws_veg_eb becomes negative (unphysical behavior), 
-!--                let the water enter the liquid water reservoir as dew on the
-!--                plant
-                   IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
-                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
-                      qsws_veg_eb(j,i) = 0.0_wp
-                   ENDIF
-                ENDIF                    
-  
-                tend = - qsws_liq_eb(j,i) * drho_l_lv
-
-                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
-                                                   + tsc(3) * tm_liq_eb_m(j,i) )
-
-!
-!--             Check if reservoir is overfull -> reduce to maximum
-!--             (conservation of water is violated here)
-                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
-
-!
-!--             Check if reservoir is empty (avoid values < 0.0)
-!--             (conservation of water is violated here)
-                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
-
-
-!
-!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
-                IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                   IF ( intermediate_timestep_count == 1 )  THEN
-                      tm_liq_eb_m(j,i) = tend
-                   ELSEIF ( intermediate_timestep_count <                      &
-                            intermediate_timestep_count_max )  THEN
-                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
-                                      * tm_liq_eb_m(j,i)
-                   ENDIF
-                ENDIF
-
-             ENDIF
-
-          ENDDO
-       ENDDO
-
-!
-!--    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    &
-            == intermediate_timestep_count_max-1 )  THEN
-#if defined( __parallel )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
-                              1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
-#else
-          force_radiation_call = force_radiation_call_l
-#endif
-          force_radiation_call_l = .FALSE.
-       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
+       BACKSPACE ( 11 )
+
+!
+!--    Read user-defined namelist
+       READ ( 11, lsm_par )
+
+!
+!--    Set flag that indicates that the land surface model is switched on
+       land_surface = .TRUE.
+
+ 10    CONTINUE
+       
+
+    END SUBROUTINE lsm_parin
 
 
@@ -1555 +2108 @@
                                      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)))
+                      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) +    &
@@ -1687 +2241 @@
 !
 !
-!--                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.
+!--                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
@@ -1696 +2250 @@
                    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.
+!--                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
@@ -1718 +2273 @@
                       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
+                            root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i)      &
+                                                            / m_total
                          ELSE
                             root_extr(k) = 0.0_wp
@@ -1792 +2348 @@
     END SUBROUTINE lsm_soil_model
 
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Swapping of timelevels
+!------------------------------------------------------------------------------!
+    SUBROUTINE lsm_swap_timelevel ( mod_count )
+
+       IMPLICIT NONE
+
+       INTEGER, INTENT(IN) :: mod_count
+
+#if defined( __nopointer )
+
+       t_surface    = t_surface_p
+       t_soil       = t_soil_p
+       IF ( humidity )  THEN
+          m_soil    = m_soil_p
+          m_liq_eb  = m_liq_eb_p
+       ENDIF
+
+#else
+    
+       SELECT CASE ( mod_count )
+
+          CASE ( 0 )
+
+             t_surface  => t_surface_1; t_surface_p  => t_surface_2
+             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
+             IF ( humidity )  THEN
+                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
+                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
+             ENDIF
+
+
+          CASE ( 1 )
+
+             t_surface  => t_surface_2; t_surface_p  => t_surface_1
+             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
+             IF ( humidity )  THEN
+                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
+                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
+             ENDIF
+
+       END SELECT
+#endif
+
+    END SUBROUTINE lsm_swap_timelevel
+
 
 !------------------------------------------------------------------------------!
@@ -1909 +2514 @@
 
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Swapping of timelevels
-!------------------------------------------------------------------------------!
-    SUBROUTINE lsm_swap_timelevel ( mod_count )
-
-       IMPLICIT NONE
-
-       INTEGER, INTENT(IN) :: mod_count
-
-#if defined( __nopointer )
-
-       t_surface    = t_surface_p
-       t_soil       = t_soil_p
-       IF ( humidity )  THEN
-          m_soil    = m_soil_p
-          m_liq_eb  = m_liq_eb_p
-       ENDIF
-
-#else
-    
-       SELECT CASE ( mod_count )
-
-          CASE ( 0 )
-
-             t_surface  => t_surface_1; t_surface_p  => t_surface_2
-             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
-             IF ( humidity )  THEN
-                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
-                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
-             ENDIF
-
-
-          CASE ( 1 )
-
-             t_surface  => t_surface_2; t_surface_p  => t_surface_1
-             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
-             IF ( humidity )  THEN
-                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
-                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
-             ENDIF
-
-       END SELECT
-#endif
-
-    END SUBROUTINE lsm_swap_timelevel
-
-
  END MODULE land_surface_model_mod