source: palm/trunk/SOURCE/land_surface_model_mod.f90 @ 4110

!> @file land_surface_model_mod.f90
!------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2019 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
! 
! 
! Former revisions:
! -----------------
! $Id: land_surface_model_mod.f90 4110 2019-07-22 17:05:21Z suehring $
! Relax checks for non-consistent initialization in case static or dynamic 
! input is provided. For example, soil_temperature or deep_soil_temperature
! is not mandatory any more if dynamic input is available. Also, improper 
! settings of x_type in namelist are only checked if no static file is 
! available. 
! 
! 4109 2019-07-22 17:00:34Z suehring
! Further revision of last commit in order to avoid any side effects when 
! albedo type is not set in namelist and default albedo type changes. 
! 
! 4024 2019-06-12 14:06:46Z suehring
! Bugfix in albedo initialization, caused crashes in rrtmg calls
! 
! 3987 2019-05-22 09:52:13Z kanani
! Introduce alternative switch for debug output during timestepping
! 
! 3964 2019-05-09 09:48:32Z suehring
! In a nested child domain, distinguish between soil moisture and temperature 
! initialization from parent via dynamic input file. Further, initialize soil
! moisture/temperature from dynamic input file only when initialization via
! 'inifor' is desired.
! 
! 3943 2019-05-02 09:50:41Z maronga
! Removed extra blank character
! 
! 3941 2019-04-30 09:48:33Z suehring
! Check that at least one surface type is set at surface element.
! 
! 3933 2019-04-25 12:33:20Z kanani
! Remove unused subroutine and allocation of pt_2m, this is done in surface_mod
! now (surfaces%pt_2m)
! 
! 
! Changes related to global restructuring of location messages and introduction 
! of additional debug messages
! 
! 3881 2019-04-10 09:31:22Z suehring
! Bugfix in level 3 initialization of pavement albedo type and pavement 
! emissivity
! 
! 3868 2019-04-08 11:52:36Z suehring
! More strict limitation of roughness length when it is in the order of the 
! vertical grid spacing 
! 
! 3856 2019-04-03 11:06:59Z suehring
! Bugfix in lsm_init in case no surface-fractions are provided
! 
! 3847 2019-04-01 14:51:44Z suehring
! Adjust message-call for checks that are especially carried out locally. 
! 
! 3832 2019-03-28 13:16:58Z raasch
! instrumented with openmp directives
! 
! 3786 2019-03-06 16:58:03Z raasch
! further unused variables removed
! 
! 3767 2019-02-27 08:18:02Z raasch
! unused variable for file index removed from rrd-subroutines parameter list
! 
! 3715 2019-02-04 17:34:55Z suehring
! Revise check for saturation moisture
! 
! 3710 2019-01-30 18:11:19Z suehring
! Check if soil-, water-, pavement- and vegetation types are set within a valid
! range.
! 
! 3692 2019-01-23 14:45:49Z suehring
! Revise check for soil moisture higher than its saturation value
! 
! 3685 2019-01-21 01:02:11Z knoop
! Some interface calls moved to module_interface + cleanup
! 
! 3677 2019-01-17 09:07:06Z moh.hefny
! Removed most_method
! 
! 3655 2019-01-07 16:51:22Z knoop
! nopointer option removed
! 
! 3620 2018-12-11 12:29:43Z moh.hefny
! update the 3d rad_lw_out array
!
! 3597 2018-12-04 08:40:18Z maronga 
! Added pt_2m / theta_2m. Removed unncessary _eb strings.
! 
! 3486 2018-11-05 06:20:18Z maronga
! Bugfix for liquid water treatment on pavements
! 
! 3361 2018-10-16 20:39:37Z knoop
! Bugfix in initialization of soil properties from dynamic input file
! 
! 3347 2018-10-15 14:21:08Z suehring
! Assign real value instead of integer
! 
! 3341 2018-10-15 10:31:27Z suehring
! Modularization of all bulk cloud physics code components
! 
! 3271 2018-09-24 08:20:34Z suehring
! Several bugfixes:
! - Initialization of pt_surface array with soil temperature in the uppermost
!   soil layer, else heat fluxes at the very first time step become quite large
! - Default initialization of vertical surface elements in special case terrain
!   height changes are larger than adjacent building heights. 
! 
! 3256 2018-09-17 12:20:07Z suehring
! Enable initialization of z0q for vegetation, pavement and water surfaces via
! namelist input. 
! 
! 3248 2018-09-14 09:42:06Z sward
! Minor formating changes
! 
! 3246 2018-09-13 15:14:50Z sward
! Added error handling for input namelist via parin_fail_message
! 
! 3241 2018-09-12 15:02:00Z raasch
! unused variables removed
! 
! 3233 2018-09-07 13:21:24Z schwenkel
! Adapted for the use of cloud_droplets 
! 
! 3222 2018-08-30 13:35:35Z suehring
! - Introduction of surface array for type and its names
! - Bugfix in intialization of pavement surfaces
! 
! 3215 2018-08-29 09:58:59Z suehring
! Enable optional initialization of soil properties directly from dynamic 
! input file.
! 
! 3209 2018-08-27 16:58:37Z suehring
! Added maximum aerodynamic resistance of 300.
! 
! 3161 2018-07-23 09:30:10Z maronga
! Increased roughness of asphalt surfaces to account for turbulence production
! by traffic and other obstacles
! 
! 3152 2018-07-19 13:26:52Z suehring
! Further adjustments for q_surface
! 
! 3147 2018-07-18 22:38:11Z maronga
! Adjustments for new surface structure
! 
! 3146 2018-07-18 22:36:19Z maronga
! Modified calculation of the surface resistance r_s
! 
! 3142 2018-07-17 15:27:45Z suehring
! Minor bugfix for last commit.
! 
! 3138 2018-07-17 08:21:20Z maronga
! Bugfix: limit roughness lengths in case of sea surface with constant_roughness
! = .F.
! 
! 3136 2018-07-16 14:48:21Z suehring
! Limit also roughness length for heat and moisture where necessary;
! Limit surface resistance to positive values
! 
! 3133 2018-07-16 11:46:50Z maronga
! Bugfix for last commit.
! 
! Some adjustments for pavement parameters
! Limit magnus formula to avoid negative q_s (leads to model crash)
! 
! 3091 2018-06-28 16:20:35Z suehring
! Add check for local roughness length not exceeding surface-layer height and
! limit roughness length where necessary.
! 
! 3051 2018-05-30 17:43:55Z suehring
! Bugfix in surface-element loops for pavement surfaces 
! 
! 3045 2018-05-28 07:55:41Z Giersch
! Error messages revised
! 
! 3045 2018-05-28 07:55:41Z Giersch
! Error messages revised and added
! 
! 3026 2018-05-22 10:30:53Z schwenkel
! Changed the name specific humidity to mixing ratio, since we are computing
! mixing ratios.
! 
! 3014 2018-05-09 08:42:38Z maronga
! Bugfix: set some initial values
! Bugfix: domain bounds of local_pf corrected
! 
! 3004 2018-04-27 12:33:25Z Giersch
! Further allocation checks implemented (averaged data will be assigned to fill
! values if no allocation happened so far) 
! 
! 2968 2018-04-13 11:52:24Z suehring
! Bugfix in initialization in case of elevated model surface
! 
! 2963 2018-04-12 14:47:44Z suehring
! - In initialization of surface types, consider the case that surface_fractions
!   is not given in static input file.
! - Introduce index for vegetation/wall, pavement/green-wall and water/window 
!   surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
! 
! 2938 2018-03-27 15:52:42Z suehring
! Initialization of soil moisture and temperature via Inifor-provided data also
! in nested child domains, even if no dynamic input file is available for 
! child domain. 1D soil profiles are received from parent model.  
! 
! 2936 2018-03-27 14:49:27Z suehring
! renamed lsm_par to land_surface_parameters. Bugfix in message calls
! 
! 2921 2018-03-22 15:05:23Z Giersch
! The activation of spinup has been moved to parin
! 
! 2894 2018-03-15 09:17:58Z Giersch
! Calculations of the index range of the subdomain on file which overlaps with
! the current subdomain are already done in read_restart_data_mod,
! lsm_read/write_restart_data was renamed to lsm_r/wrd_local, USE kinds has 
! been removed in several routines, variable named found has been
! introduced for checking if restart data was found, reading of restart strings
! has been moved completely to read_restart_data_mod, lsm_rrd_local is already
! inside the overlap loop programmed in read_restart_data_mod, the marker ***
! end lsm *** is not necessary anymore, strings and their respective lengths 
! are written out and read now in case of restart runs to get rid of prescribed
! character lengths, SAVE attribute added where necessary, deallocation and 
! allocation of some arrays have been changed to take care of different restart
! files that can be opened (index i)
! 
! 2881 2018-03-13 16:24:40Z suehring
! Bugfix: wrong loop structure for soil moisture calculation
! 
! 2805 2018-02-14 17:00:09Z suehring
! Bugfix in initialization of roughness over water surfaces
! 
! 2798 2018-02-09 17:16:39Z suehring
! Minor bugfix for initialization of pt_surface
! 
! 2797 2018-02-08 13:24:35Z suehring
! Move output of ghf to general 2D output to output ghf also at urban-type 
! surfaces. 
! Move restart data of ghf_av to read/write_3d_binary, as this is not a 
! exclusively LSM variable anymore.   
! 
! 2765 2018-01-22 11:34:58Z maronga
! Major bugfix in calculation of f_shf for vertical surfaces
! 
! 2735 2018-01-11 12:01:27Z suehring
! output of r_a moved from land-surface to consider also urban-type surfaces
! 
! 2729 2018-01-09 11:22:28Z maronga
! Separated deep soil temperature from soil_temperature array
! 
! 2724 2018-01-05 12:12:38Z maronga
! Added security check for insufficient soil_temperature values
! 
! 2723 2018-01-05 09:27:03Z maronga
! Bugfix for spinups (end_time was increased twice in case of LSM + USM runs)
! 
! 2718 2018-01-02 08:49:38Z maronga
! Corrected "Former revisions" section
! 
! 2707 2017-12-18 18:34:46Z suehring
! Changes from last commit documented
! 
! 2706 2017-12-18 18:33:49Z suehring
! Bugfix, read surface temperature in case of restart runs.
!
! 2705 2017-12-18 11:26:23Z maronga
! Bugfix in binary output (wrong sequence)
! 
! 2696 2017-12-14 17:12:51Z kanani
! Change in file header (GPL part)
! Bugfix: missing USE statement for calc_mean_profile
! do not write surface temperatures onto pt array as this might cause 
! problems with nesting (MS)
! Revised calculation of pt1 and qv1 (now done in surface_layer_fluxes). Bugfix
! in calculation of surface density (cannot be done via an surface non-air 
! temperature) (BM)
! Bugfix: g_d was NaN for non-vegetaed surface types (BM)
! Bugfix initialization of c_veg and lai
! Revise data output to enable _FillValues
! Bugfix in calcultion of r_a and rad_net_l (MS)
! Bugfix: rad_net is not updated in case of radiation_interaction and must thu
! be calculated again from the radiative fluxes
! Temporary fix for cases where no soil model is used on some PEs (BM)
! Revised input and initialization of soil and surface paramters
! pavement_depth is variable for each surface element
! radiation quantities belong to surface type now
! surface fractions initialized
! Rename lsm_last_actions into lsm_wrd_subdomain (MS)
! 
! 2608 2017-11-13 14:04:26Z schwenkel
! Calculation of magnus equation in external module (diagnostic_quantities_mod).
! Adjust calculation of vapor pressure and saturation mixing ratio that it is
! consistent with formulations in other parts of PALM.
! 
! 2575 2017-10-24 09:57:58Z maronga
! Pavement parameterization revised
! 
! 2573 2017-10-20 15:57:49Z scharf
! bugfixes in last_actions
! 
! 2548 2017-10-16 13:18:20Z suehring
! extended by cloud_droplets option
! 
! 2532 2017-10-11 16:00:46Z scharf
! bugfixes in data_output_3d
! 
! 2516 2017-10-04 11:03:04Z suehring
! Remove tabs
! 
! 2514 2017-10-04 09:52:37Z suehring
! upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny
! no output of ghost layer data
! 
! 2504 2017-09-27 10:36:13Z maronga
! Support roots and water under pavement. Added several pavement types.
! 
! 2476 2017-09-18 07:54:32Z maronga
! Bugfix for last commit
! 
! 2475 2017-09-18 07:42:36Z maronga
! Bugfix: setting of vegetation_pars for bare soil corrected.
! 
! 2354 2017-08-17 10:49:36Z schwenkel
! minor bugfixes
! 
! 2340 2017-08-07 17:11:13Z maronga
! Revised root_distribution tabel and implemented a pseudo-generic root fraction
! calculation
! 
! 2333 2017-08-04 09:08:26Z maronga
! minor bugfixes
! 
! 2332 2017-08-03 21:15:22Z maronga
! bugfix in pavement_pars
! 
! 2328 2017-08-03 12:34:22Z maronga
! Revised skin layer concept.
! Bugfix for runs with pavement surface and humidity
! Revised some standard values in vegetation_pars
! Added emissivity and default albedo_type as variable to tables
! Changed default surface type to vegetation
! Revised input of soil layer configuration
! 
! 2307 2017-07-07 11:32:10Z suehring
! Bugfix, variable names corrected
! 
! 2299 2017-06-29 10:14:38Z maronga
! Removed pt_p from USE statement. Adjusted call to lsm_soil_model to allow
! spinups without soil moisture prediction
! 
! 2298 2017-06-29 09:28:18Z raasch
! type of write_binary changed from CHARACTER to LOGICAL
! 
! 2296 2017-06-28 07:53:56Z maronga
! Bugfix in calculation of bare soil heat capacity.
! Bugfix in calculation of shf
! Added support for spinups
! 
! 2282 2017-06-13 11:38:46Z schwenkel
! Bugfix for check of saturation moisture
! 
! 2273 2017-06-09 12:46:06Z sward
! Error number changed
! 
! 2270 2017-06-09 12:18:47Z maronga
! Revised parameterization of heat conductivity between skin layer and soil.
! Temperature and moisture are now defined at the center of the layers.
! Renamed veg_type to vegetation_type and pave_type to pavement_type_name
! Renamed and reduced the number of look-up tables (vegetation_pars, soil_pars)
! Revised land surface model initialization
! Removed output of shf_eb and qsws_eb and removed _eb throughout code
! Removed Clapp & Hornberger parameterization
! 
! 2249 2017-06-06 13:58:01Z sward
! 
! 2248 2017-06-06 13:52:54Z sward $
! Error no changed
!
! 2246 2017-06-06 13:09:34Z sward
! Error no changed
!
! Changed soil configuration to 8 layers. The number of soil layers is now
! freely adjustable via the NAMELIST.
! 
! 2237 2017-05-31 10:34:53Z suehring
! Bugfix in write restart data
! 
! 2233 2017-05-30 18:08:54Z suehring
!
! 2232 2017-05-30 17:47:52Z suehring
! Adjustments to new topography and surface concept
!   - now, also vertical walls are possible
!   - for vertical walls, parametrization of r_a (aerodynamic resisistance) is 
!     implemented. 
!
! Add check for soil moisture, it must not exceed its saturation value.
! 
! 2149 2017-02-09 16:57:03Z scharf
! Land surface parameters II corrected for vegetation_type 18 and 19
! 
! 2031 2016-10-21 15:11:58Z knoop
! renamed variable rho to rho_ocean
! 
! 2000 2016-08-20 18:09:15Z knoop
! Forced header and separation lines into 80 columns
! 
! 1978 2016-07-29 12:08:31Z maronga
! Bugfix: initial values of pave_surface and water_surface were not set.
! 
! 1976 2016-07-27 13:28:04Z maronga
! Parts of the code have been reformatted. Use of radiation model output is
! generalized and simplified. Added more output quantities due to modularization
! 
! 1972 2016-07-26 07:52:02Z maronga
! Further modularization: output of cross sections and 3D data is now done in this
! module. Moreover, restart data is written and read directly within this module.
!
! 
! 1966 2016-07-18 11:54:18Z maronga
! Bugfix: calculation of m_total in soil model was not set to zero at model start
! 
! 1949 2016-06-17 07:19:16Z maronga
! Bugfix: calculation of qsws_soil_eb with precipitation = .TRUE. gave
! qsws_soil_eb = 0 due to a typo
! 
! 1856 2016-04-13 12:56:17Z maronga
! Bugfix: for water surfaces, the initial water surface temperature is set equal
! to the intital skin temperature. Moreover, the minimum value of r_a is now
! 1.0 to avoid too large fluxes at the first model time step
! 
! 1849 2016-04-08 11:33:18Z hoffmann
! prr moved to arrays_3d
!
! 1826 2016-04-07 12:01:39Z maronga
! Cleanup after modularization
! 
! 1817 2016-04-06 15:44:20Z maronga
! Added interface for lsm_init_arrays. Added subroutines for check_parameters,
! header, and parin. Renamed some subroutines.
! 
! 1788 2016-03-10 11:01:04Z maronga
! Bugfix: calculate lambda_surface based on temperature gradient between skin
! layer and soil layer instead of Obukhov length
! Changed: moved calculation of surface specific humidity to energy balance solver
! New: water surfaces are available by using a fixed sea surface temperature.
! The roughness lengths are calculated dynamically using the Charnock 
! parameterization. This involves the new roughness length for moisture z0q.
! New: modified solution of the energy balance solver and soil model for
! paved surfaces (i.e. asphalt concrete).
! Syntax layout improved.
! Changed: parameter dewfall removed.
! 
! 1783 2016-03-06 18:36:17Z raasch
! netcdf variables moved to netcdf module
!
! 1757 2016-02-22 15:49:32Z maronga
! Bugfix: set tm_soil_m to zero after allocation. Added parameter 
! unscheduled_radiation_calls to control calls of the radiation model based on
! the skin temperature change during one time step (preliminary version). Set 
! qsws_soil_eb to zero at model start (previously set to qsws_eb). Removed MAX
! function as it cannot be vectorized.
! 
! 1709 2015-11-04 14:47:01Z maronga
! Renamed pt_1 and qv_1 to pt1 and qv1.
! Bugfix: set initial values for t_surface_p in case of restart runs
! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky'
! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky'
! Added todo action
!
! 1697 2015-10-28 17:14:10Z raasch
! bugfix: misplaced cpp-directive
!
! 1695 2015-10-27 10:03:11Z maronga
! Bugfix: REAL constants provided with KIND-attribute in call of 
! Replaced rif with ol
! 
! 1691 2015-10-26 16:17:44Z maronga
! Added skip_time_do_lsm to allow for spin-ups without LSM. Various bugfixes:
! Soil temperatures are now defined at the edges of the layers, calculation of
! shb_eb corrected, prognostic equation for skin temperature corrected. Surface
! fluxes are now directly transfered to atmosphere
! 
! 1682 2015-10-07 23:56:08Z knoop
! Code annotations made doxygen readable 
! 
! 1590 2015-05-08 13:56:27Z maronga
! Bugfix: definition of character strings requires same length for all elements
! 
! 1585 2015-04-30 07:05:52Z maronga
! Modifications for RRTMG. Changed tables to PARAMETER type. 
! 
! 1571 2015-03-12 16:12:49Z maronga
! Removed upper-case variable names. Corrected distribution of precipitation to 
! the liquid water reservoir and the bare soil fractions.
! 
! 1555 2015-03-04 17:44:27Z maronga
! Added output of r_a and r_s
! 
! 1553 2015-03-03 17:33:54Z maronga
! Improved better treatment of roughness lengths. Added default soil temperature
! profile
! 
! 1551 2015-03-03 14:18:16Z maronga
! Flux calculation is now done in prandtl_fluxes. Added support for data output.
! Vertical indices have been replaced. Restart runs are now possible. Some
! variables have beem renamed. Bugfix in the prognostic equation for the surface
! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
! the hydraulic conductivity. Bugfix for root fraction and extraction
! calculation
! 
! intrinsic function MAX and MIN
!
! 1500 2014-12-03 17:42:41Z maronga
! Corrected calculation of aerodynamic resistance (r_a).
! Precipitation is now added to liquid water reservoir using LE_liq.
! Added support for dry runs.
! 
! 1496 2014-12-02 17:25:50Z maronga
! Initial revision
! 
!
! Description:
! ------------
!> Land surface model, consisting of a solver for the energy balance at the
!> surface and a multi layer soil scheme. The scheme is similar to the TESSEL
!> scheme implemented in the ECMWF IFS model, with modifications according to
!> H-TESSEL. The implementation is based on the formulation implemented in the 
!> DALES and UCLA-LES models.
!>
!> @todo Extensive verification energy-balance solver for vertical surfaces, 
!>       e.g. parametrization of r_a
!> @todo Revise single land-surface processes for vertical surfaces, e.g. 
!>       treatment of humidity, etc. 
!> @todo Consider partial absorption of the net shortwave radiation by the 
!>       skin layer.
!> @todo Improve surface water parameterization
!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, 
!>       nzt_soil=3)).
!> @todo Implement surface runoff model (required when performing long-term LES 
!>       with considerable precipitation.
!> @todo Revise calculation of f2 when wilting point is non-constant in the
!>       soil
!> @todo Allow for zero soil moisture (currently, it is set to wilting point)
!> @note No time step criterion is required as long as the soil layers do not
!>       become too thin.
!> @todo Attention, pavement_subpars_1/2 are hardcoded to 8 levels, in case 
!>       more levels are used this may cause an potential bug
!> @todo Routine calc_q_surface required?
!> @todo Allow for precipitation water to enter pavements that are semi-pervious
!------------------------------------------------------------------------------!
 MODULE land_surface_model_mod
 
    USE arrays_3d,                                                             &
        ONLY:  hyp, pt, prr, q, q_p, ql, vpt, u, v, w, hyrho, exner, d_exner

    USE basic_constants_and_equations_mod,                                     &
        ONLY:  c_p, g, lv_d_cp, l_v, kappa, magnus, rho_l, r_d, r_v, rd_d_rv

    USE calc_mean_profile_mod,                                                 &
        ONLY:  calc_mean_profile

    USE control_parameters,                                                    &
        ONLY:  cloud_droplets, coupling_start_time,                            &
               debug_output, debug_output_timestep, debug_string,              &
               dt_3d,                                                          &
               end_time, humidity, intermediate_timestep_count,                &
               initializing_actions, intermediate_timestep_count_max,          &
               land_surface, max_masks, pt_surface,                            &
               rho_surface, spinup, spinup_pt_mean, spinup_time,               &
               surface_pressure, timestep_scheme, tsc,                         &
               time_since_reference_point

    USE indices,                                                               &
        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb

    USE bulk_cloud_model_mod,                                                  &
        ONLY: bulk_cloud_model, precipitation

    USE netcdf_data_input_mod,                                                 &
        ONLY :  building_type_f, init_3d,                                      &
                input_pids_dynamic,                                            &
                input_pids_static,                                             &
                netcdf_data_input_interpolate, netcdf_data_input_init_lsm,     &
                pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,  &
                root_area_density_lsm_f, soil_pars_f, soil_type_f,             &
                surface_fraction_f, vegetation_pars_f, vegetation_type_f,      &
                water_pars_f, water_type_f

    USE kinds

    USE pegrid

    USE radiation_model_mod,                                                   &
        ONLY:  albedo, albedo_type, emissivity, force_radiation_call,          &
               radiation, radiation_scheme, unscheduled_radiation_calls
        
    USE statistics,                                                            &
        ONLY:  hom, statistic_regions

    USE surface_mod,                                                           &
        ONLY :  ind_pav_green, ind_veg_wall, ind_wat_win,                      &
                surf_lsm_h, surf_lsm_v, surf_type, surface_restore_elements

    IMPLICIT NONE

    TYPE surf_type_lsm
       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  var_1d !< 1D prognostic variable
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2d !< 2D prognostic variable
    END TYPE surf_type_lsm

!
!-- LSM model constants

    REAL(wp), PARAMETER  ::                    &
              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil (W/m/K)  
              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix (W/m/K)
              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water (W/m/K)
              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil (J/m3/K)
              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water (J/m3/K)
              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)


    REAL(wp), DIMENSION(0:7), PARAMETER  :: dz_soil_default =                  & ! default soil layer configuration
                                            (/ 0.01_wp, 0.02_wp, 0.04_wp,      &
                                               0.06_wp, 0.14_wp, 0.26_wp,      &
                                               0.54_wp, 1.86_wp/)

    REAL(wp), DIMENSION(0:3), PARAMETER  :: dz_soil_ref =                      & ! reference four layer soil configuration used for estimating the root fractions
                                            (/ 0.07_wp, 0.21_wp, 0.72_wp,      &
                                               1.89_wp /)

    REAL(wp), DIMENSION(0:3), PARAMETER  :: zs_ref =                           & ! reference four layer soil configuration used for estimating the root fractions
                                            (/ 0.07_wp, 0.28_wp, 1.0_wp,       &
                                               2.89_wp /)


!
!-- LSM variables
    CHARACTER(10) :: surface_type = 'netcdf'      !< general classification. Allowed are:
                                                  !< 'vegetation', 'pavement', ('building'), 
                                                  !< 'water', and 'netcdf'



    INTEGER(iwp) :: nzb_soil = 0,             & !< bottom of the soil model (Earth's surface)
                    nzt_soil = 7,             & !< top of the soil model
                    nzt_pavement = 0,         & !< top of the pavement within the soil
                    nzs = 8,                  & !< number of soil layers
                    pavement_depth_level = 0, & !< default NAMELIST nzt_pavement
                    pavement_type = 1,        & !< default NAMELIST pavement_type                 
                    soil_type = 3,            & !< default NAMELIST soil_type
                    vegetation_type = 2,      & !< default NAMELIST vegetation_type
                    water_type = 1              !< default NAMELISt water_type
                    
   
        
    LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
               constant_roughness = .FALSE.,     & !< use fixed/dynamic roughness lengths for water surfaces
               force_radiation_call_l = .FALSE., & !< flag to force calling of radiation routine
               aero_resist_kray = .TRUE.           !< flag to control parametrization of aerodynamic resistance at vertical surface elements

!   value 9999999.9_wp -> generic available or user-defined value must be set
!   otherwise -> no generic variable and user setting is optional
    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
                canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
                c_surface = 9999999.9_wp,               & !< Surface (skin) heat capacity (J/m2/K)
                deep_soil_temperature =  9999999.9_wp,  & !< Deep soil temperature (bottom boundary condition)
                drho_l_lv,                              & !< (rho_l * l_v)**-1
                field_capacity = 9999999.9_wp,          & !< NAMELIST m_fc
                f_shortwave_incoming = 9999999.9_wp,    & !< NAMELIST f_sw_in
                hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
                ke = 0.0_wp,                            & !< Kersten number
                lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil (W/m/K)
                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s (W/m2/K)
                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u (W/m2/K)
                leaf_area_index = 9999999.9_wp,         & !< NAMELIST lai
                l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
                min_canopy_resistance = 9999999.9_wp,   & !< NAMELIST r_canopy_min
                min_soil_resistance = 50.0_wp,          & !< NAMELIST r_soil_min
                m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
                pavement_heat_capacity = 9999999.9_wp,  & !< volumetric heat capacity of pavement (e.g. roads) (J/m3/K)
                pavement_heat_conduct  = 9999999.9_wp,  & !< heat conductivity for pavements (e.g. roads) (W/m/K)
                q_s = 0.0_wp,                           & !< saturation water vapor mixing ratio
                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
                rho_cp,                                 & !< rho_surface * cp
                rho_lv,                                 & !< rho_ocean * l_v
                saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
                skip_time_do_lsm = 0.0_wp,              & !< LSM is not called before this time
                vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
                water_temperature = 9999999.9_wp,       & !< water temperature
                wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
                z0_vegetation  = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
                z0h_vegetation = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
                z0q_vegetation = 9999999.9_wp,          & !< NAMELIST z0q (lsm_par)
                z0_pavement    = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
                z0h_pavement   = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
                z0q_pavement   = 9999999.9_wp,          & !< NAMELIST z0q (lsm_par)
                z0_water       = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
                z0h_water      = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
                z0q_water      = 9999999.9_wp             !< NAMELIST z0q (lsm_par)  
                
                
    REAL(wp), DIMENSION(:), ALLOCATABLE  :: ddz_soil_center, & !< 1/dz_soil_center
                                            ddz_soil,        & !< 1/dz_soil
                                            dz_soil_center,  & !< soil grid spacing (center-center)
                                            zs,              & !< depth of the temperature/moisute levels
                                            root_extr          !< root extraction


                                            
    REAL(wp), DIMENSION(0:20)  ::  root_fraction = 9999999.9_wp,     & !< (NAMELIST) distribution of root surface area to the individual soil layers
                                   soil_moisture = 0.0_wp,           & !< NAMELIST soil moisture content (m3/m3)
                                   soil_temperature = 9999999.9_wp,  & !< NAMELIST soil temperature (K) +1
                                   dz_soil  = 9999999.9_wp,          & !< (NAMELIST) soil layer depths (spacing)
                                   zs_layer = 9999999.9_wp             !< soil layer depths (edge)
                                 
    TYPE(surf_type_lsm), POINTER ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
                                     t_soil_h_p,  & !< Prog. soil temperature (K), horizontal surface elements
                                     m_soil_h,    & !< Soil moisture (m3/m3), horizontal surface elements
                                     m_soil_h_p     !< Prog. soil moisture (m3/m3), horizontal surface elements 

    TYPE(surf_type_lsm), TARGET  ::  t_soil_h_1,  & !<
                                     t_soil_h_2,  & !<
                                     m_soil_h_1,  & !<
                                     m_soil_h_2     !<

    TYPE(surf_type_lsm), DIMENSION(:), POINTER :: &
                                     t_soil_v,    & !< Soil temperature (K), vertical surface elements
                                     t_soil_v_p,  & !< Prog. soil temperature (K), vertical surface elements
                                     m_soil_v,    & !< Soil moisture (m3/m3), vertical surface elements
                                     m_soil_v_p     !< Prog. soil moisture (m3/m3), vertical surface elements   

    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::&
                                     t_soil_v_1,  & !<
                                     t_soil_v_2,  & !<
                                     m_soil_v_1,  & !<
                                     m_soil_v_2     !<

    TYPE(surf_type_lsm), POINTER  ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements 
                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements 
                                      m_liq_h,        & !< liquid water reservoir (m), horizontal surface elements 
                                      m_liq_h_p         !< progn. liquid water reservoir (m), horizontal surface elements 

    TYPE(surf_type_lsm), TARGET   ::  t_surface_h_1,  & !<
                                      t_surface_h_2,  & !<
                                      m_liq_h_1,      & !<
                                      m_liq_h_2         !<

    TYPE(surf_type_lsm), DIMENSION(:), POINTER  ::    &
                                      t_surface_v,    & !< surface temperature (K), vertical surface elements 
                                      t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements 
                                      m_liq_v,        & !< liquid water reservoir (m), vertical surface elements 
                                      m_liq_v_p         !< progn. liquid water reservoir (m), vertical surface elements 

    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
                                      t_surface_v_1,  & !<
                                      t_surface_v_2,  & !<
                                      m_liq_v_1,      & !<
                                      m_liq_v_2         !<

    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av

    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  t_soil_av, & !< Average of t_soil
                                                        m_soil_av    !< Average of m_soil

    TYPE(surf_type_lsm), TARGET ::  tm_liq_h_m      !< liquid water reservoir tendency (m), horizontal surface elements 
    TYPE(surf_type_lsm), TARGET ::  tt_surface_h_m  !< surface temperature tendency (K), horizontal surface elements 
    TYPE(surf_type_lsm), TARGET ::  tt_soil_h_m     !< t_soil storage array, horizontal surface elements 
    TYPE(surf_type_lsm), TARGET ::  tm_soil_h_m     !< m_soil storage array, horizontal surface elements 

    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_liq_v_m      !< liquid water reservoir tendency (m), vertical surface elements 
    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_surface_v_m  !< surface temperature tendency (K), vertical surface elements 
    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_soil_v_m     !< t_soil storage array, vertical surface elements 
    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_soil_v_m     !< m_soil storage array, vertical surface elements 

!
!-- Energy balance variables                
    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
              c_liq_av,         & !< average of c_liq
              c_soil_av,        & !< average of c_soil
              c_veg_av,         & !< average of c_veg
              lai_av,           & !< average of lai       
              qsws_liq_av,      & !< average of qsws_liq
              qsws_soil_av,     & !< average of qsws_soil
              qsws_veg_av,      & !< average of qsws_veg
              r_s_av              !< average of r_s
  
!
!-- Predefined Land surface classes (vegetation_type)
    CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ &
                                   'user defined              ',           & !  0 
                                   'bare soil                 ',           & !  1                           
                                   'crops, mixed farming      ',           & !  2
                                   'short grass               ',           & !  3
                                   'evergreen needleleaf trees',           & !  4
                                   'deciduous needleleaf trees',           & !  5
                                   'evergreen broadleaf trees ',           & !  6
                                   'deciduous broadleaf trees ',           & !  7
                                   'tall grass                ',           & !  8
                                   'desert                    ',           & !  9
                                   'tundra                    ',           & ! 10
                                   'irrigated crops           ',           & ! 11
                                   'semidesert                ',           & ! 12
                                   'ice caps and glaciers     ',           & ! 13
                                   'bogs and marshes          ',           & ! 14
                                   'evergreen shrubs          ',           & ! 15
                                   'deciduous shrubs          ',           & ! 16
                                   'mixed forest/woodland     ',           & ! 17
                                   'interrupted forest        '            & ! 18
                                                                 /)

!
!-- Soil model classes (soil_type)
    CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ &
                                   'user defined',                  & ! 0 
                                   'coarse      ',                  & ! 1
                                   'medium      ',                  & ! 2
                                   'medium-fine ',                  & ! 3
                                   'fine        ',                  & ! 4
                                   'very fine   ',                  & ! 5
                                   'organic     '                   & ! 6
                                                                 /)

!
!-- Pavement classes
    CHARACTER(29), DIMENSION(0:15), PARAMETER :: pavement_type_name = (/ &
                                   'user defined                 ', & ! 0 
                                   'asphalt/concrete mix         ', & ! 1
                                   'asphalt (asphalt concrete)   ', & ! 2
                                   'concrete (Portland concrete) ', & ! 3
                                   'sett                         ', & ! 4
                                   'paving stones                ', & ! 5
                                   'cobblestone                  ', & ! 6
                                   'metal                        ', & ! 7
                                   'wood                         ', & ! 8
                                   'gravel                       ', & ! 9
                                   'fine gravel                  ', & ! 10
                                   'pebblestone                  ', & ! 11
                                   'woodchips                    ', & ! 12
                                   'tartan (sports)              ', & ! 13
                                   'artifical turf (sports)      ', & ! 14
                                   'clay (sports)                '  & ! 15
                                                                 /)                                                              
                                                                 
!
!-- Water classes
    CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ &
                                   'user defined',                   & ! 0 
                                   'lake        ',                   & ! 1
                                   'river       ',                   & ! 2
                                   'ocean       ',                   & ! 3
                                   'pond        ',                   & ! 4
                                   'fountain    '                    & ! 5
                                                                  /)                                                                                  
                    
!
!-- Land surface parameters according to the respective classes (vegetation_type)
    INTEGER(iwp) ::  ind_v_rc_min = 0    !< index for r_canopy_min in vegetation_pars
    INTEGER(iwp) ::  ind_v_rc_lai = 1    !< index for LAI in vegetation_pars
    INTEGER(iwp) ::  ind_v_c_veg   = 2   !< index for c_veg in vegetation_pars
    INTEGER(iwp) ::  ind_v_gd  = 3       !< index for g_d in vegetation_pars
    INTEGER(iwp) ::  ind_v_z0 = 4        !< index for z0 in vegetation_pars
    INTEGER(iwp) ::  ind_v_z0qh = 5      !< index for z0h / z0q in vegetation_pars
    INTEGER(iwp) ::  ind_v_lambda_s = 6  !< index for lambda_s_s in vegetation_pars
    INTEGER(iwp) ::  ind_v_lambda_u = 7  !< index for lambda_s_u in vegetation_pars
    INTEGER(iwp) ::  ind_v_f_sw_in = 8   !< index for f_sw_in in vegetation_pars
    INTEGER(iwp) ::  ind_v_c_surf = 9    !< index for c_surface in vegetation_pars
    INTEGER(iwp) ::  ind_v_at = 10       !< index for albedo_type in vegetation_pars
    INTEGER(iwp) ::  ind_v_emis = 11     !< index for emissivity in vegetation_pars

    INTEGER(iwp) ::  ind_w_temp     = 0    !< index for temperature in water_pars
    INTEGER(iwp) ::  ind_w_z0       = 1    !< index for z0 in water_pars
    INTEGER(iwp) ::  ind_w_z0h      = 2    !< index for z0h in water_pars
    INTEGER(iwp) ::  ind_w_lambda_s = 3    !< index for lambda_s_s in water_pars
    INTEGER(iwp) ::  ind_w_lambda_u = 4    !< index for lambda_s_u in water_pars
    INTEGER(iwp) ::  ind_w_at       = 5    !< index for albedo type in water_pars
    INTEGER(iwp) ::  ind_w_emis     = 6    !< index for emissivity in water_pars

    INTEGER(iwp) ::  ind_p_z0       = 0    !< index for z0 in pavement_pars
    INTEGER(iwp) ::  ind_p_z0h      = 1    !< index for z0h in pavement_pars
    INTEGER(iwp) ::  ind_p_at       = 2    !< index for albedo type in pavement_pars
    INTEGER(iwp) ::  ind_p_emis     = 3    !< index for emissivity in pavement_pars
    INTEGER(iwp) ::  ind_p_lambda_h = 0    !< index for lambda_h in pavement_subsurface_pars
    INTEGER(iwp) ::  ind_p_rho_c    = 1    !< index for rho_c in pavement_pars
!
!-- Land surface parameters
!-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface, albedo_type, emissivity
    REAL(wp), DIMENSION(0:11,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  0.005_wp,   0.5E-4_wp,     0.0_wp,    0.0_wp, 0.00_wp, 0.00_wp, 17.0_wp, 0.94_wp, & !  1
        180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp,   0.10_wp,    0.001_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  2.0_wp, 0.95_wp, & !  2
        110.0_wp, 2.00_wp, 1.00_wp, 0.00_wp,   0.03_wp,   0.3E-4_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  2.0_wp, 0.95_wp, & !  3
        500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  5.0_wp, 0.97_wp, & !  4
        500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  6.0_wp, 0.97_wp, & !  5
        175.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  8.0_wp, 0.97_wp, & !  6
        240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  9.0_wp, 0.97_wp, & !  7
        100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,   0.47_wp,  0.47E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  8.0_wp, 0.97_wp, & !  8
        250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  0.013_wp, 0.013E-2_wp,    15.0_wp,   15.0_wp, 0.00_wp, 0.00_wp,  3.0_wp, 0.94_wp, & !  9
         80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  0.034_wp, 0.034E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp, 11.0_wp, 0.97_wp, & ! 10
        180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp,    0.5_wp,  0.50E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp, 13.0_wp, 0.97_wp, & ! 11
        150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,   0.17_wp,  0.17E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  2.0_wp, 0.97_wp, & ! 12
          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 1.3E-3_wp,   1.3E-4_wp,    58.0_wp,   58.0_wp, 0.00_wp, 0.00_wp, 11.0_wp, 0.97_wp, & ! 13
        240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,   0.83_wp,  0.83E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  4.0_wp, 0.97_wp, & ! 14
        225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,   0.10_wp,  0.10E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  4.0_wp, 0.97_wp, & ! 15
        225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  4.0_wp, 0.97_wp, & ! 16
        250.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  7.0_wp, 0.97_wp, & ! 17
        175.0_wp, 2.50_wp, 1.00_wp, 0.03_wp,   1.10_wp,     1.10_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  8.0_wp, 0.97_wp  & ! 18
                                                               /), (/ 12, 18 /) )

                                    
!
!-- Root distribution for default soil layer configuration (sum = 1)
!--                                level 1 - level 4 according to zs_ref
    REAL(wp), DIMENSION(0:3,1:18), PARAMETER :: root_distribution = RESHAPE( (/ &
                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  1
                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  2
                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  3
                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  4
                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  5
                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  6
                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  7
                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  8
                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  9
                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & ! 10
                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 11
                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 12
                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 13
                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 14
                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 15 
                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16 
                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 17
                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 18
                                 /), (/ 4, 18 /) )

!
!-- Soil parameters according to the following porosity classes (soil_type)

!
!-- Soil parameters  alpha_vg,      l_vg,    n_vg, gamma_w_sat,    m_sat,     m_fc,   m_wilt,    m_res 
    REAL(wp), DIMENSION(0:7,1:6), PARAMETER :: soil_pars = RESHAPE( (/     &
                      3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp,& ! 1
                      3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp,& ! 2
                      0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp,& ! 3
                      3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp,& ! 4
                      2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp,& ! 5
                      1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp & ! 6
                                                                     /), (/ 8, 6 /) )


!
!-- TO BE FILLED
!-- Pavement parameters      z0,       z0h, albedo_type, emissivity  
    REAL(wp), DIMENSION(0:3,1:15), PARAMETER :: pavement_pars = RESHAPE( (/ &
                      5.0E-2_wp, 5.0E-4_wp,     18.0_wp,    0.97_wp,  & !  1
                      5.0E-2_wp, 5.0E-4_wp,     19.0_wp,    0.94_wp,  & !  2
                      1.0E-2_wp, 1.0E-4_wp,     20.0_wp,    0.98_wp,  & !  3                                  
                      1.0E-2_wp, 1.0E-4_wp,     21.0_wp,    0.93_wp,  & !  4
                      1.0E-2_wp, 1.0E-4_wp,     22.0_wp,    0.97_wp,  & !  5
                      1.0E-2_wp, 1.0E-4_wp,     23.0_wp,    0.97_wp,  & !  6
                      1.0E-2_wp, 1.0E-4_wp,     24.0_wp,    0.97_wp,  & !  7
                      1.0E-2_wp, 1.0E-4_wp,     25.0_wp,    0.94_wp,  & !  8
                      1.0E-2_wp, 1.0E-4_wp,     26.0_wp,    0.98_wp,  & !  9                                  
                      1.0E-2_wp, 1.0E-4_wp,     27.0_wp,    0.93_wp,  & ! 10
                      1.0E-2_wp, 1.0E-4_wp,     28.0_wp,    0.97_wp,  & ! 11
                      1.0E-2_wp, 1.0E-4_wp,     29.0_wp,    0.97_wp,  & ! 12
                      1.0E-2_wp, 1.0E-4_wp,     30.0_wp,    0.97_wp,  & ! 13
                      1.0E-2_wp, 1.0E-4_wp,     31.0_wp,    0.94_wp,  & ! 14
                      1.0E-2_wp, 1.0E-4_wp,     32.0_wp,    0.98_wp   & ! 15
                      /), (/ 4, 15 /) )                             
!
!-- Pavement subsurface parameters part 1: thermal conductivity (W/m/K)
!--   0.0-0.01, 0.01-0.03, 0.03-0.07, 0.07-0.15, 0.15-0.30, 0.30-0.50,    0.50-1.25,    1.25-3.00
    REAL(wp), DIMENSION(0:7,1:15), PARAMETER :: pavement_subsurface_pars_1 = RESHAPE( (/ &
       0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp, 9999999.9_wp, 9999999.9_wp, & !  1
       0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp, 9999999.9_wp, 9999999.9_wp, & !  2
       0.89_wp,   0.89_wp,   0.89_wp,   0.89_wp,   0.89_wp,   0.89_wp, 9999999.9_wp, 9999999.9_wp, & !  3
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & !  4
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & !  5
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & !  6
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & !  7
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & !  8
       0.70_wp,   0.70_wp,   0.70_wp,   0.70_wp,   0.70_wp,   0.70_wp, 9999999.9_wp, 9999999.9_wp, & !  9
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 10
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 11
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 12
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 13
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & ! 14
       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp  & ! 15
       /), (/ 8, 15 /) )

!
!-- Pavement subsurface parameters part 2: volumetric heat capacity (J/m3/K)
!--     0.0-0.01, 0.01-0.03, 0.03-0.07, 0.07-0.15, 0.15-0.30, 0.30-0.50,    0.50-1.25,    1.25-3.00
    REAL(wp), DIMENSION(0:7,1:15), PARAMETER :: pavement_subsurface_pars_2 = RESHAPE( (/ &
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & !  1
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & !  2
       1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 9999999.9_wp, 9999999.9_wp, & !  3
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & !  4
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & !  5
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & !  6
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & !  7
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & !  8
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & !  9
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 10
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 11
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 12
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 13
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & ! 14
       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp  & ! 15
                           /), (/ 8, 15 /) )
 
!
!-- TO BE FILLED
!-- Water parameters                    temperature,     z0,      z0h, albedo_type, emissivity,
    REAL(wp), DIMENSION(0:6,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 1
       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 2
       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 3
       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 4
       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp  & ! 5
                                                                     /), (/ 7, 5 /) )                                                                   
                                                                                                                                      
    SAVE


    PRIVATE

    
!
!-- Public functions
    PUBLIC lsm_boundary_condition, lsm_check_data_output,                      &
           lsm_check_data_output_pr,                                           &
           lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,& 
           lsm_data_output_2d, lsm_data_output_3d, lsm_energy_balance,         &
           lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model,   &
           lsm_swap_timelevel, lsm_rrd_local, lsm_wrd_local
! !vegetat
!-- Public parameters, constants and initial values
    PUBLIC aero_resist_kray, skip_time_do_lsm

!
!-- Public grid variables
    PUBLIC nzb_soil, nzs, nzt_soil, zs

!
!-- Public prognostic variables
    PUBLIC m_soil_h, t_soil_h

    INTERFACE lsm_boundary_condition
       MODULE PROCEDURE lsm_boundary_condition
    END INTERFACE lsm_boundary_condition

    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_3d_data_averaging
       MODULE PROCEDURE lsm_3d_data_averaging
    END INTERFACE lsm_3d_data_averaging

    INTERFACE lsm_data_output_2d
       MODULE PROCEDURE lsm_data_output_2d
    END INTERFACE lsm_data_output_2d

    INTERFACE lsm_data_output_3d
       MODULE PROCEDURE lsm_data_output_3d
    END INTERFACE lsm_data_output_3d

    INTERFACE lsm_define_netcdf_grid
       MODULE PROCEDURE lsm_define_netcdf_grid
    END INTERFACE lsm_define_netcdf_grid

    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
    END INTERFACE lsm_soil_model

    INTERFACE lsm_swap_timelevel
       MODULE PROCEDURE lsm_swap_timelevel
    END INTERFACE lsm_swap_timelevel

    INTERFACE lsm_rrd_local
       MODULE PROCEDURE lsm_rrd_local
    END INTERFACE lsm_rrd_local

    INTERFACE lsm_wrd_local
       MODULE PROCEDURE lsm_wrd_local
    END INTERFACE lsm_wrd_local

 CONTAINS


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Set internal Neumann boundary condition at outer soil grid points 
!> for temperature and humidity. 
!------------------------------------------------------------------------------!
 SUBROUTINE lsm_boundary_condition
 
    IMPLICIT NONE

    INTEGER(iwp) :: i      !< grid index x-direction
    INTEGER(iwp) :: ioff   !< offset index x-direction indicating location of soil grid point
    INTEGER(iwp) :: j      !< grid index y-direction
    INTEGER(iwp) :: joff   !< offset index x-direction indicating location of soil grid point
    INTEGER(iwp) :: k      !< grid index z-direction
    INTEGER(iwp) :: koff   !< offset index x-direction indicating location of soil grid point
    INTEGER(iwp) :: l      !< running index surface-orientation
    INTEGER(iwp) :: m      !< running index surface elements

    koff = surf_lsm_h%koff
    DO  m = 1, surf_lsm_h%ns
       i = surf_lsm_h%i(m)
       j = surf_lsm_h%j(m)
       k = surf_lsm_h%k(m)
       pt(k+koff,j,i) = pt(k,j,i)
    ENDDO

    DO  l = 0, 3
       ioff = surf_lsm_v(l)%ioff
       joff = surf_lsm_v(l)%joff
       DO  m = 1, surf_lsm_v(l)%ns
          i = surf_lsm_v(l)%i(m)
          j = surf_lsm_v(l)%j(m)
          k = surf_lsm_v(l)%k(m)
          pt(k,j+joff,i+ioff) = pt(k,j,i)
       ENDDO
    ENDDO
!
!-- In case of humidity, set boundary conditions also for q and vpt.
    IF ( humidity )  THEN
       koff = surf_lsm_h%koff
       DO  m = 1, surf_lsm_h%ns
          i = surf_lsm_h%i(m)
          j = surf_lsm_h%j(m)
          k = surf_lsm_h%k(m)
          q(k+koff,j,i)   = q(k,j,i)
          vpt(k+koff,j,i) = vpt(k,j,i)
       ENDDO

       DO  l = 0, 3
          ioff = surf_lsm_v(l)%ioff
          joff = surf_lsm_v(l)%joff
          DO  m = 1, surf_lsm_v(l)%ns
             i = surf_lsm_v(l)%i(m)
             j = surf_lsm_v(l)%j(m)
             k = surf_lsm_v(l)%k(m)
             q(k,j+joff,i+ioff)   = q(k,j,i)
             vpt(k,j+joff,i+ioff) = vpt(k,j,i)
          ENDDO
       ENDDO
    ENDIF

 END SUBROUTINE lsm_boundary_condition

!------------------------------------------------------------------------------!
! 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( 'lsm_check_data_output', '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( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
          ENDIF
          unit = 'K'   
             
       CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'm_liq*',                 &
              'qsws_liq*', 'qsws_soil*', 'qsws_veg*', 'r_s*' )
          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( 'lsm_check_data_output', 'PA0111', 1, 2, 0, 6, 0 )
          ENDIF
          IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                              'res land_surface = .TRUE.'
             CALL message( 'lsm_check_data_output', 'PA0404', 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( 'lsm_check_data_output', '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( 'lsm_check_data_output', '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( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
          ENDIF
          IF ( TRIM( var ) == 'm_liq*'  .AND.  .NOT.  land_surface )  THEN
             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                              'res land_surface = .TRUE.'
             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
          ENDIF
          IF ( TRIM( var ) == 'qsws_liq*'  .AND.  .NOT. land_surface )         &
          THEN
             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                              'res land_surface = .TRUE.'
             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
          ENDIF
          IF ( TRIM( var ) == 'qsws_soil*'  .AND.  .NOT.  land_surface )       &
          THEN
             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                              'res land_surface = .TRUE.'
             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
          ENDIF
          IF ( TRIM( var ) == 'qsws_veg*'  .AND.  .NOT. land_surface )         &
          THEN
             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
                              'res land_surface = .TRUE.'
             CALL message( 'lsm_check_data_output', '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( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( TRIM( var ) == 'lai*'   )      unit = 'none' 
          IF ( TRIM( var ) == 'c_liq*' )      unit = 'none'
          IF ( TRIM( var ) == 'c_soil*')      unit = 'none'
          IF ( TRIM( var ) == 'c_veg*' )      unit = 'none'
          IF ( TRIM( var ) == 'm_liq*'     )  unit = 'm'
          IF ( TRIM( var ) == 'qsws_liq*'  )  unit = 'W/m2'
          IF ( TRIM( var ) == 'qsws_soil*' )  unit = 'W/m2'
          IF ( TRIM( var ) == 'qsws_veg*'  )  unit = 'W/m2'
          IF ( TRIM( var ) == 'r_s*')         unit = 's/m' 
             
       CASE DEFAULT
          unit = 'illegal'

    END SELECT


 END SUBROUTINE lsm_check_data_output



!------------------------------------------------------------------------------!
! Description:
! ------------
!> Check data output of profiles for land surface model
!------------------------------------------------------------------------------!
 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) ::  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( 'lsm_check_data_output_pr', 'PA0402', 1, 2, 0, 6, 0 )
          ELSE
             dopr_index(var_count) = 89
             dopr_unit     = 'K'
             hom(0:nzs-1,2,89,:)  = SPREAD( - zs(nzb_soil:nzt_soil), 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(nzb_soil:nzt_soil), 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( 'lsm_check_data_output_pr', 'PA0402', 1, 2, 0, 6, 0 )
          ELSE
             dopr_index(var_count) = 91
             dopr_unit     = 'm3/m3'
             hom(0:nzs-1,2,91,:)  = SPREAD( - zs(nzb_soil:nzt_soil), 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(nzb_soil:nzt_soil), 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
                     
    
    IMPLICIT NONE

    INTEGER(iwp) ::  i        !< running index, x-dimension
    INTEGER(iwp) ::  j        !< running index, y-dimension
    INTEGER(iwp) ::  k        !< running index, z-dimension

!
!-- Check for a valid setting of surface_type. The default value is 'netcdf'.
!-- In that case, the surface types are read from NetCDF file
    IF ( TRIM( surface_type ) /= 'vegetation'  .AND.                           &
         TRIM( surface_type ) /= 'pavement'    .AND.                           &
         TRIM( surface_type ) /= 'water'       .AND.                           &
         TRIM( surface_type ) /= 'netcdf' )  THEN  
       message_string = 'unknown surface type: surface_type = "' //            &
                        TRIM( surface_type ) // '"'
       CALL message( 'lsm_check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
    ENDIF

!
!-- 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( 'lsm_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( 'lsm_check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
    ENDIF
    
    IF (  .NOT.  radiation )  THEN
       message_string = 'lsm requires '//                                      &
                        'the radiation model to be switched on'
       CALL message( 'lsm_check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
    ENDIF
!
!-- Check if soil types are set within a valid range.
    IF ( TRIM( surface_type ) == 'vegetation'  .OR.                            &
         TRIM( surface_type ) == 'pavement' )  THEN
       IF ( soil_type < LBOUND( soil_pars, 2 )  .AND.                          &
            soil_type > UBOUND( soil_pars, 2 ) )  THEN
          WRITE( message_string, * ) 'soil_type = ', soil_type, ' is out ' //  &
                                     'of the valid range'
          CALL message( 'lsm_check_parameters', 'PA0452', 2, 2, 0, 6, 0 )
       ENDIF
    ENDIF
    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
       IF ( soil_type_f%from_file )  THEN
          DO  i = nxl, nxr
             DO  j = nys, nyn
                IF ( soil_type_f%var_2d(j,i) /= soil_type_f%fill  .AND.        &
                     ( soil_type_f%var_2d(j,i) < LBOUND( soil_pars, 2 )  .OR.  &
                       soil_type_f%var_2d(j,i) > UBOUND( soil_pars, 2 ) ) )  THEN
                   WRITE( message_string, * ) 'soil_type = is out  of ' //     &
                                        'the valid range at (j,i) = ', j, i
                   CALL message( 'lsm_check_parameters', 'PA0452',             &
                                  2, 2, myid, 6, 0 )
                ENDIF
             ENDDO
          ENDDO
       ENDIF
    ENDIF
!
!-- Check if vegetation types are set within a valid range.    
    IF ( TRIM( surface_type ) == 'vegetation' )  THEN
       IF ( vegetation_type < LBOUND( vegetation_pars, 2 )  .AND.              &
            vegetation_type > UBOUND( vegetation_pars, 2 ) )  THEN
          WRITE( message_string, * ) 'vegetation_type = ', vegetation_type,    &
                                     ' is out of the valid range'
          CALL message( 'lsm_check_parameters', 'PA0526', 2, 2, 0, 6, 0 )
       ENDIF
    ENDIF
    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
       IF ( vegetation_type_f%from_file )  THEN
          DO  i = nxl, nxr
             DO  j = nys, nyn
                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill  .AND.&
              ( vegetation_type_f%var(j,i) < LBOUND( vegetation_pars, 2 )  .OR.&
                vegetation_type_f%var(j,i) > UBOUND( vegetation_pars, 2 ) ) )  &
                THEN
                   WRITE( message_string, * ) 'vegetation_type = is out of ' //&
                                        'the valid range at (j,i) = ', j, i
                   CALL message( 'lsm_check_parameters', 'PA0526',             &
                                  2, 2, myid, 6, 0 ) 
                ENDIF
             ENDDO
          ENDDO
       ENDIF
    ENDIF
!
!-- Check if pavement types are set within a valid range.    
    IF ( TRIM( surface_type ) == 'pavement' )  THEN
       IF ( pavement_type < LBOUND( pavement_pars, 2 )  .AND.                  &
            pavement_type > UBOUND( pavement_pars, 2 ) )  THEN
          WRITE( message_string, * ) 'pavement_type = ', pavement_type,        &
                                     ' is out of the valid range'
          CALL message( 'lsm_check_parameters', 'PA0527', 2, 2, 0, 6, 0 )
       ENDIF
    ENDIF
    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
       IF ( pavement_type_f%from_file )  THEN
          DO  i = nxl, nxr
             DO  j = nys, nyn
                IF ( pavement_type_f%var(j,i) /= pavement_type_f%fill  .AND.   &
              ( pavement_type_f%var(j,i) < LBOUND( pavement_pars, 2 )  .OR.    &
                pavement_type_f%var(j,i) > UBOUND( pavement_pars, 2 ) ) )  THEN
                   WRITE( message_string, * ) 'pavement_type = is out of ' //  &
                                        'the valid range at (j,i) = ', j, i
                   CALL message( 'lsm_check_parameters', 'PA0527',             &
                                  2, 2, myid, 6, 0 ) 
                ENDIF
             ENDDO
          ENDDO
       ENDIF
    ENDIF
!
!-- Check if water types are set within a valid range.    
    IF ( TRIM( surface_type ) == 'water' )  THEN
       IF ( water_type < LBOUND( water_pars, 2 )  .AND.                        &
            water_type > UBOUND( water_pars, 2 ) )  THEN
          WRITE( message_string, * ) 'water_type = ', water_type,              &
                                     ' is out of the valid range'
          CALL message( 'lsm_check_parameters', 'PA0528', 2, 2, 0, 6, 0 )
       ENDIF
    ENDIF
    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
       IF ( water_type_f%from_file )  THEN
          DO  i = nxl, nxr
             DO  j = nys, nyn
                IF ( water_type_f%var(j,i) /= water_type_f%fill  .AND.         &
              ( water_type_f%var(j,i) < LBOUND( water_pars, 2 )  .OR.          &
                water_type_f%var(j,i) > UBOUND( water_pars, 2 ) ) )  THEN
                   WRITE( message_string, * ) 'water_type = is out  of ' //    &
                                        'the valid range at (j,i) = ', j, i
                   CALL message( 'lsm_check_parameters', 'PA0528',             &
                                 2, 2, myid, 6, 0 ) 
                ENDIF
             ENDDO
          ENDDO
       ENDIF
    ENDIF
!
!-- Check further settings for consistency.
    IF ( TRIM( surface_type ) == 'vegetation' )  THEN
    
       IF ( vegetation_type == 0 )  THEN
          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
             message_string = 'vegetation_type = 0 (user defined)'//           &
                              'requires setting of min_canopy_resistance'//    &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( leaf_area_index == 9999999.9_wp )  THEN
             message_string = 'vegetation_type = 0 (user_defined)'//           &
                              'requires setting of leaf_area_index'//          &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( vegetation_coverage == 9999999.9_wp )  THEN
             message_string = 'vegetation_type = 0 (user_defined)'//           &
                              'requires setting of vegetation_coverage'//      &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
             message_string = 'vegetation_type = 0 (user_defined)'//           &
                              'requires setting of'//                          &
                              'canopy_resistance_coefficient /= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
             message_string = 'vegetation_type = 0 (user_defined)'//           &
                              'requires setting of lambda_surface_stable'//    &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
             message_string = 'vegetation_type = 0 (user_defined)'//           &
                              'requires setting of lambda_surface_unstable'//  &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
             message_string = 'vegetation_type = 0 (user_defined)'//           &
                              'requires setting of f_shortwave_incoming'//     &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( z0_vegetation == 9999999.9_wp )  THEN
             message_string = 'vegetation_type = 0 (user_defined)'//           &
                              'requires setting of z0_vegetation'//            &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( z0h_vegetation == 9999999.9_wp )  THEN
             message_string = 'vegetation_type = 0 (user_defined)'//           &
                              'requires setting of z0h_vegetation'//           &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
          ENDIF
       ENDIF

       IF ( vegetation_type == 1 )  THEN
          IF ( vegetation_coverage /= 9999999.9_wp  .AND.  vegetation_coverage &
               /= 0.0_wp )  THEN
             message_string = 'vegetation_type = 1 (bare soil)'//              &
                              ' requires vegetation_coverage = 0'
             CALL message( 'lsm_check_parameters', 'PA0294', 1, 2, 0, 6, 0 )
          ENDIF
       ENDIF
  
    ENDIF
    
    IF ( TRIM( surface_type ) == 'water' )  THEN

       IF ( water_type == 0 )  THEN  
       
          IF ( z0_water == 9999999.9_wp )  THEN
             message_string = 'water_type = 0 (user_defined)'//                &
                              'requires setting of z0_water'//                 &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
          ENDIF

          IF ( z0h_water == 9999999.9_wp )  THEN
             message_string = 'water_type = 0 (user_defined)'//                &
                              'requires setting of z0h_water'//                &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
          ENDIF
          
          IF ( water_temperature == 9999999.9_wp )  THEN
             message_string = 'water_type = 0 (user_defined)'//                &
                              'requires setting of water_temperature'//        &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0379', 1, 2, 0, 6, 0 )
          ENDIF       
          
       ENDIF
       
    ENDIF
    
    IF ( TRIM( surface_type ) == 'pavement' )  THEN

       IF ( ANY( dz_soil /= 9999999.9_wp )  .AND.  pavement_type /= 0 )  THEN
          message_string = 'non-default setting of dz_soil '//                  &
                           'does not allow to use pavement_type /= 0)'
          CALL message( 'lsm_check_parameters', 'PA0341', 1, 2, 0, 6, 0 )
       ENDIF

       IF ( pavement_type == 0 )  THEN  
       
          IF ( z0_pavement == 9999999.9_wp )  THEN
             message_string = 'pavement_type = 0 (user_defined)'//             &
                              'requires setting of z0_pavement'//              &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0352', 1, 2, 0, 6, 0 )
          ENDIF
          
          IF ( z0h_pavement == 9999999.9_wp )  THEN
             message_string = 'pavement_type = 0 (user_defined)'//             &
                              'requires setting of z0h_pavement'//             &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0353', 1, 2, 0, 6, 0 )
          ENDIF
          
          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
             message_string = 'pavement_type = 0 (user_defined)'//             &
                              'requires setting of pavement_heat_conduct'//    &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0342', 1, 2, 0, 6, 0 )
          ENDIF 
          
           IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
             message_string = 'pavement_type = 0 (user_defined)'//             &
                              'requires setting of pavement_heat_capacity'//   &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
          ENDIF 

          IF ( pavement_depth_level == 0 )  THEN
             message_string = 'pavement_type = 0 (user_defined)'//             &
                              'requires setting of pavement_depth_level'//     &
                              '/= 0'
             CALL message( 'lsm_check_parameters', 'PA0474', 1, 2, 0, 6, 0 )
          ENDIF 

       ENDIF 
    
    ENDIF

    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
!
!--    MS: Some problme here, after calling message everythings stucks at
!--        MPI_FINALIZE call. 
       IF ( ANY( pavement_type_f%var /= pavement_type_f%fill )  .AND.           &
            ANY( dz_soil /= 9999999.9_wp ) )  THEN
          message_string = 'pavement-surfaces are not allowed in ' //           &
                           'combination with a non-default setting of dz_soil'
          CALL message( 'lsm_check_parameters', 'PA0999', 2, 2, 0, 6, 0 )
       ENDIF
    ENDIF
    
!
!-- Temporary message as long as NetCDF input is not available
    IF ( TRIM( surface_type ) == 'netcdf'  .AND.  .NOT. input_pids_static )   &
    THEN
       message_string = 'surface_type = netcdf requires static input file.'
       CALL message( 'lsm_check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
    ENDIF

    IF ( soil_type == 0  .AND.  .NOT. input_pids_static )  THEN

       IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
          message_string = 'soil_type = 0 (user_defined)'//                    &
                           'requires setting of alpha_vangenuchten'//          &
                           '/= 9999999.9'
          CALL message( 'lsm_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( 'lsm_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( 'lsm_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( 'lsm_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( 'lsm_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( 'lsm_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( 'lsm_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( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
       ENDIF

    ENDIF


!!! these checks are not needed for water surfaces??

!
!-- Determine number of soil layers to be used and check whether an appropriate
!-- root fraction is prescribed
    nzb_soil = 0
    nzt_soil = -1
    IF ( ALL( dz_soil == 9999999.9_wp ) )  THEN
       nzt_soil = 7
       dz_soil(nzb_soil:nzt_soil) = dz_soil_default
    ELSE
       DO k = 0, 19
          IF ( dz_soil(k) /= 9999999.9_wp )  THEN
             nzt_soil = nzt_soil + 1
          ENDIF
       ENDDO    
    ENDIF
    nzs = nzt_soil + 1

!
!-- Check whether valid soil temperatures are prescribed
    IF ( .NOT. input_pids_dynamic )  THEN
       IF ( COUNT( soil_temperature /= 9999999.9_wp ) /= nzs )  THEN
          WRITE( message_string, * )                                           &
                                  'number of soil layers (', nzs, ') does not',&
                                  ' match to the number of layers specified',  &
                                  ' in soil_temperature (', COUNT(             &
                                  soil_temperature /= 9999999.9_wp ), ')'
             CALL message( 'lsm_check_parameters', 'PA0471', 1, 2, 0, 6, 0 )
       ENDIF

       IF ( deep_soil_temperature == 9999999.9_wp ) THEN
             message_string = 'deep_soil_temperature is not set but must be'// &
                              '/= 9999999.9'
             CALL message( 'lsm_check_parameters', 'PA0472', 1, 2, 0, 6, 0 )
       ENDIF
    ENDIF

!
!-- Check whether the sum of all root fractions equals one
    IF ( vegetation_type == 0 )  THEN
       IF ( SUM( root_fraction(nzb_soil:nzt_soil) ) /= 1.0_wp )  THEN
          message_string = 'vegetation_type = 0 (user_defined)'//              &
                           'requires setting of root_fraction'//               &
                           '/= 9999999.9 and SUM(root_fraction) = 1'
          CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
       ENDIF
    ENDIF     
!
!-- Calculate grid spacings. Temperature and moisture are defined at
!-- the center of the soil layers, whereas gradients/fluxes are 
!-- defined at the edges (_layer)
!
!-- Allocate global 1D arrays
    ALLOCATE ( ddz_soil_center(nzb_soil:nzt_soil) )
    ALLOCATE ( ddz_soil(nzb_soil:nzt_soil+1) )
    ALLOCATE ( dz_soil_center(nzb_soil:nzt_soil) )
    ALLOCATE ( zs(nzb_soil:nzt_soil+1) )


    zs(nzb_soil) = 0.5_wp * dz_soil(nzb_soil)
    zs_layer(nzb_soil) = dz_soil(nzb_soil)

    DO  k = nzb_soil+1, nzt_soil
       zs_layer(k) = zs_layer(k-1) + dz_soil(k)
       zs(k) = (zs_layer(k) +  zs_layer(k-1)) * 0.5_wp
    ENDDO

    dz_soil(nzt_soil+1) = zs_layer(nzt_soil) + dz_soil(nzt_soil)
    zs(nzt_soil+1) = zs_layer(nzt_soil) + 0.5_wp * dz_soil(nzt_soil)
 
    DO  k = nzb_soil, nzt_soil-1
       dz_soil_center(k) = zs(k+1) - zs(k)
       IF ( dz_soil_center(k) <= 0.0_wp )  THEN
          message_string = 'invalid soil layer configuration found ' //        &
                           '(dz_soil_center(k) <= 0.0)'
          CALL message( 'lsm_rrd_local', 'PA0140', 1, 2, 0, 6, 0 )
       ENDIF  
    ENDDO
  
    dz_soil_center(nzt_soil) = zs_layer(k-1) + dz_soil(k) - zs(nzt_soil)
       
    ddz_soil_center = 1.0_wp / dz_soil_center
    ddz_soil(nzb_soil:nzt_soil) = 1.0_wp / dz_soil(nzb_soil:nzt_soil)



 END SUBROUTINE lsm_check_parameters
 
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Solver for the energy balance at the surface.
!------------------------------------------------------------------------------!
 SUBROUTINE lsm_energy_balance( horizontal, l )

    USE pegrid
    USE radiation_model_mod,  ONLY:  rad_lw_out

    IMPLICIT NONE

    INTEGER(iwp) ::  i         !< running index
    INTEGER(iwp) ::  i_off     !< offset to determine index of surface element, seen from atmospheric grid point, for x
    INTEGER(iwp) ::  j         !< running index
    INTEGER(iwp) ::  j_off     !< offset to determine index of surface element, seen from atmospheric grid point, for y
    INTEGER(iwp) ::  k         !< running index
    INTEGER(iwp) ::  k_off     !< offset to determine index of surface element, seen from atmospheric grid point, for z
    INTEGER(iwp) ::  ks        !< running index
    INTEGER(iwp) ::  l         !< surface-facing index
    INTEGER(iwp) ::  m         !< running index concerning wall elements

    LOGICAL      ::  horizontal !< Flag indicating horizontal or vertical surfaces

    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
                f_qsws_veg,  & !< factor for qsws_veg
                f_qsws_soil, & !< factor for qsws_soil
                f_qsws_liq,  & !< factor for qsws_liq
                f_shf,       & !< factor for shf
                lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
                lambda_surface, & !< Current value of lambda_surface (W/m2/K)
                m_liq_max      !< maxmimum value of the liq. water reservoir

    TYPE(surf_type_lsm), POINTER ::  surf_t_surface
    TYPE(surf_type_lsm), POINTER ::  surf_t_surface_p
    TYPE(surf_type_lsm), POINTER ::  surf_tt_surface_m
    TYPE(surf_type_lsm), POINTER ::  surf_m_liq
    TYPE(surf_type_lsm), POINTER ::  surf_m_liq_p
    TYPE(surf_type_lsm), POINTER ::  surf_tm_liq_m

    TYPE(surf_type_lsm), POINTER ::  surf_m_soil
    TYPE(surf_type_lsm), POINTER ::  surf_t_soil

    TYPE(surf_type), POINTER  ::  surf  !< surface-date type variable


    IF ( debug_output_timestep )  THEN
       WRITE( debug_string, * ) 'lsm_energy_balance', horizontal, l
       CALL debug_message( debug_string, 'start' )
    ENDIF

    IF ( horizontal )  THEN
       surf              => surf_lsm_h

       surf_t_surface    => t_surface_h
       surf_t_surface_p  => t_surface_h_p
       surf_tt_surface_m => tt_surface_h_m
       surf_m_liq        => m_liq_h
       surf_m_liq_p      => m_liq_h_p
       surf_tm_liq_m     => tm_liq_h_m
       surf_m_soil       => m_soil_h
       surf_t_soil       => t_soil_h
    ELSE
       surf              => surf_lsm_v(l)

       surf_t_surface    => t_surface_v(l)
       surf_t_surface_p  => t_surface_v_p(l)
       surf_tt_surface_m => tt_surface_v_m(l)
       surf_m_liq        => m_liq_v(l)
       surf_m_liq_p      => m_liq_v_p(l)
       surf_tm_liq_m     => tm_liq_v_m(l)
       surf_m_soil       => m_soil_v(l)
       surf_t_soil       => t_soil_v(l)
    ENDIF

!
!-- Index offset of surface element point with respect to adjoining 
!-- atmospheric grid point 
    k_off = surf%koff
    j_off = surf%joff
    i_off = surf%ioff

    !$OMP PARALLEL PRIVATE (m, i, j, k, lambda_h_sat, ke, lambda_soil, lambda_surface,             &
    !$OMP&                  c_surface_tmp, f1,m_total, f2, e_s, e, f3, m_min, m_liq_max, q_s,      &
    !$OMP&                  f_qsws_veg, f_qsws_soil, f_qsws_liq, f_shf, f_qsws, e_s_dt, dq_s_dt,   &
    !$OMP&                  coef_1, coef_2, tend)
    !$OMP DO SCHEDULE (STATIC)
    DO  m = 1, surf%ns

       i   = surf%i(m)            
       j   = surf%j(m)
       k   = surf%k(m)

!
!--    Define heat conductivity between surface and soil depending on surface
!--    type. For vegetation, a skin layer parameterization is used. The new
!--    parameterization uses a combination of two conductivities: a constant
!--    conductivity for the skin layer, and a conductivity according to the
!--    uppermost soil layer. For bare soil and pavements, no skin layer is 
!--    applied. In these cases, the temperature is assumed to be constant 
!--    between the surface and the first soil layer. The heat conductivity is 
!--    then derived from the soil/pavement properties. 
!--    For water surfaces, the conductivity is already set to 1E10.
!--    Moreover, the heat capacity is set. For bare soil the heat capacity is
!--    the capacity of the uppermost soil layer, for pavement it is that of
!--    the material involved.

!
!--    for vegetation type surfaces, the thermal conductivity of the soil is
!--    needed

       IF ( surf%vegetation_surface(m) )  THEN

          lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(nzb_soil,m)) *      &
                         lambda_h_water ** surf_m_soil%var_2d(nzb_soil,m)
                         
          ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(nzb_soil,m) /   &
                                                     surf%m_sat(nzb_soil,m) ) )                   
                         
          lambda_soil = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry )   &
                           * ddz_soil(nzb_soil) * 2.0_wp

!
!--       When bare soil is set without a thermal conductivity (no skin layer),
!--       a heat capacity is that of the soil layer, otherwise it is a
!--       combination of the conductivities from the skin and the soil layer
          IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN 
            surf%c_surface(m) = (rho_c_soil * (1.0_wp - surf%m_sat(nzb_soil,m))&
                              + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) ) &
                              * dz_soil(nzb_soil) * 0.5_wp    
            lambda_surface = lambda_soil

          ELSE IF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m))&
          THEN
             lambda_surface = surf%lambda_surface_s(m) * lambda_soil           &
                              / ( surf%lambda_surface_s(m) + lambda_soil )
          ELSE

             lambda_surface = surf%lambda_surface_u(m) * lambda_soil           &
                              / ( surf%lambda_surface_u(m) + lambda_soil )
          ENDIF
       ELSE
          lambda_surface = surf%lambda_surface_s(m)
       ENDIF

!
!--    Set heat capacity of the skin/surface. It is ususally zero when a skin
!--    layer is used, and non-zero otherwise.
       c_surface_tmp = surf%c_surface(m) 

!
!--    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 ( bulk_cloud_model )  THEN
!           pt1 = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
!           qv1 = q(k,j,i) - ql(k,j,i)
!        ELSEIF ( cloud_droplets ) THEN
!           pt1 = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
!           qv1 = q(k,j,i) 
!        ELSE
!           pt1 = pt(k,j,i)
!           IF ( humidity )  THEN
!              qv1 = q(k,j,i)
!           ELSE
!              qv1 = 0.0_wp
!           ENDIF
!        ENDIF
!
!--     Calculation of r_a for vertical surfaces
!--
!--     heat transfer coefficient for forced convection along vertical walls
!--     follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
!--            
!--       H = httc (Tsfc - Tair)
!--       httc = rw * (11.8 + 4.2 * Ueff) - 4.0
!--            
!--             rw: wall patch roughness relative to 1.0 for concrete
!--             Ueff: effective wind speed
!--             - 4.0 is a reduction of Rowley et al (1930) formulation based on
!--             Cole and Sturrock (1977)
!--           
!--             Ucan: Canyon wind speed
!--             wstar: convective velocity
!--             Qs: surface heat flux
!--             zH: height of the convective layer
!--             wstar = (g/Tcan*Qs*zH)**(1./3.)
                
!--    Effective velocity components must always 
!--    be defined at scalar grid point. The wall normal component is 
!--    obtained by simple linear interpolation. ( An alternative would
!--    be an logarithmic interpolation. )
!--    A roughness lenght of 0.001 is assumed for concrete (the inverse,
!--    1000 is used in the nominator for scaling)
!--    To do: detailed investigation which approach gives more reliable results!
!--    Please note, in case of very small friction velocity, e.g. in little 
!--    holes, the resistance can become negative. For this reason, limit r_a 
!--    to positive values. 
       IF ( horizontal  .OR.  .NOT. aero_resist_kray )  THEN
          surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) /            &
                             ( surf%ts(m) * surf%us(m) + 1.0E-20_wp ) )
       ELSE
          surf%r_a(m) = rho_cp / ( surf%z0(m) * 1000.0_wp                      &
                        * ( 11.8_wp + 4.2_wp *                                 &
                        SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
                                   ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
                                   ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2,  &
                              0.01_wp ) )                                      &
                           )  - 4.0_wp  ) 
       ENDIF
!
!--    Make sure that the resistance does not drop to zero for neutral 
!--    stratification. Also, set a maximum resistance to avoid the breakdown of
!--    MOST for locations with zero wind speed
       IF ( surf%r_a(m) <   1.0_wp )  surf%r_a(m) =   1.0_wp
       IF ( surf%r_a(m) > 300.0_wp )  surf%r_a(m) = 300.0_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)
       f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) /          &
                        (0.81_wp * (0.004_wp * surf%rad_sw_in(m)               &
                         + 1.0_wp)) )

!
!--    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 + surf%root_fr(ks,m)                              &
                     * MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(ks,m) )
       ENDDO 

!
!--    The calculation of f2 is based on only one wilting point value for all
!--    soil layers. The value at k=nzb_soil is used here as a proxy but might
!--    need refinement in the future.
       IF ( m_total > surf%m_wilt(nzb_soil,m)  .AND.                           &
            m_total < surf%m_fc(nzb_soil,m) )  THEN
          f2 = ( m_total - surf%m_wilt(nzb_soil,m) ) /                         &
               ( surf%m_fc(nzb_soil,m) - surf%m_wilt(nzb_soil,m) )
       ELSEIF ( m_total >= surf%m_fc(nzb_soil,m) )  THEN
          f2 = 1.0_wp
       ELSE
          f2 = 1.0E-20_wp
       ENDIF

!
!--    Calculate water vapour pressure at saturation and convert to hPa
!--    The magnus formula is limited to temperatures up to 333.15 K to
!--    avoid negative values of q_s
       e_s = 0.01_wp * magnus( MIN(surf_t_surface%var_1d(m), 333.15_wp) )

!
!--    f3: correction for vapour pressure deficit
       IF ( surf%g_d(m) /= 0.0_wp )  THEN
!
!--       Calculate vapour pressure
          e  = surf%qv1(m) * surface_pressure / ( surf%qv1(m) + rd_d_rv )
          f3 = EXP ( - surf%g_d(m) * (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
       surf%r_canopy(m) = surf%r_canopy_min(m) /                               &
                              ( surf%lai(m) * f1 * f2 * f3 + 1.0E-20_wp )
!
!--    Third step: calculate bare soil resistance r_soil.
       m_min = surf%c_veg(m) * surf%m_wilt(nzb_soil,m) +                       &
                         ( 1.0_wp - surf%c_veg(m) ) * surf%m_res(nzb_soil,m)


       f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) /                       &
            ( surf%m_fc(nzb_soil,m) - m_min )
       f2 = MAX( f2, 1.0E-20_wp )
       f2 = MIN( f2, 1.0_wp     )

       surf%r_soil(m) = surf%r_soil_min(m) / 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 ( surf%pavement_surface(m) )  THEN
          m_liq_max = m_max_depth * 5.0_wp
          surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq%var_1d(m) / m_liq_max)**0.67 )
       ELSE
          m_liq_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m)              &
                      + ( 1.0_wp - surf%c_veg(m) ) )
          surf%c_liq(m) = MIN( 1.0_wp, surf_m_liq%var_1d(m) / m_liq_max )
       ENDIF
!
!--    Calculate saturation water vapor mixing ratio
       q_s = rd_d_rv * e_s / ( surface_pressure - e_s )
!
!--    In case of dewfall, set evapotranspiration to zero
!--    All super-saturated water is then removed from the air
       IF ( humidity  .AND.  q_s <= surf%qv1(m) )  THEN
          surf%r_canopy(m) = 0.0_wp
          surf%r_soil(m)   = 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 ( surf%water_surface(m) )  THEN
          f_qsws_veg  = 0.0_wp
          f_qsws_soil = 0.0_wp
          f_qsws_liq  = rho_lv / surf%r_a(m)
       ELSEIF ( surf%pavement_surface(m) )  THEN
          f_qsws_veg  = 0.0_wp
          f_qsws_soil = 0.0_wp
          f_qsws_liq  = rho_lv * surf%c_liq(m) / surf%r_a(m)
       ELSE
          f_qsws_veg  = rho_lv * surf%c_veg(m) *                               &
                            ( 1.0_wp        - surf%c_liq(m)    ) /             &
                            ( surf%r_a(m) + surf%r_canopy(m) )
          f_qsws_soil = rho_lv * (1.0_wp    - surf%c_veg(m)    ) /             &
                            ( surf%r_a(m) + surf%r_soil(m)   )
          f_qsws_liq  = rho_lv * surf%c_veg(m) * surf%c_liq(m)   /             &
                              surf%r_a(m)
       ENDIF

       f_shf  = rho_cp / surf%r_a(m)
       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.62_wp / ( surf_t_surface%var_1d(m) - 29.65_wp) -   &
                        17.62_wp*( surf_t_surface%var_1d(m) - 273.15_wp)      &
                       / ( surf_t_surface%var_1d(m) - 29.65_wp)**2 )

       dq_s_dt = rd_d_rv * e_s_dt / ( surface_pressure - e_s_dt )
!
!--    Calculate net radiation radiation without longwave outgoing flux because
!--    it has a dependency on surface temperature and thus enters the prognostic 
!--    equations directly
       surf%rad_net_l(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)              &
                           + surf%rad_lw_in(m)
!
!--    Calculate new skin temperature
       IF ( humidity )  THEN
!
!--       Numerator of the prognostic equation
          coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m)             &
                   * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)             &
                   + f_shf * surf%pt1(m) + f_qsws * ( surf%qv1(m) - q_s        &
                   + dq_s_dt * surf_t_surface%var_1d(m) ) + lambda_surface     &
                   * surf_t_soil%var_2d(nzb_soil,m)

!
!--       Denominator of the prognostic equation
          coef_2 = surf%rad_lw_out_change_0(m) + f_qsws * dq_s_dt              &
                   + lambda_surface + f_shf / exner(nzb)
       ELSE
!
!--       Numerator of the prognostic equation
          coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m)             &
                   * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)             &
                   + f_shf * surf%pt1(m)  + lambda_surface                     &
                   * surf_t_soil%var_2d(nzb_soil,m)
!
!--       Denominator of the prognostic equation
          coef_2 = surf%rad_lw_out_change_0(m) + lambda_surface + f_shf / exner(nzb)

       ENDIF

       tend = 0.0_wp

!
!--    Implicit solution when the surface layer has no heat capacity,
!--    otherwise use RK3 scheme.
       surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *&
                          surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2&
                                             * dt_3d * tsc(2) ) 

!
!--    Add RK3 term
       IF ( c_surface_tmp /= 0.0_wp )  THEN

          surf_t_surface_p%var_1d(m) = surf_t_surface_p%var_1d(m) + dt_3d *    &
                                       tsc(3) * surf_tt_surface_m%var_1d(m)

!
!--       Calculate true tendency
          tend = ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) -     &
                   dt_3d * tsc(3) * surf_tt_surface_m%var_1d(m)) / (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
                surf_tt_surface_m%var_1d(m) = tend
             ELSEIF ( intermediate_timestep_count <                            &
                      intermediate_timestep_count_max )  THEN
                surf_tt_surface_m%var_1d(m) = -9.5625_wp * tend +              &
                                               5.3125_wp * surf_tt_surface_m%var_1d(m)
             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( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )       &
            > 0.2_wp  .AND. &
            unscheduled_radiation_calls )  THEN
          force_radiation_call_l = .TRUE.
       ENDIF

       surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exner(nzb)

!
!--    Calculate fluxes
       surf%rad_net_l(m) = surf%rad_net_l(m) +                                 &
                            surf%rad_lw_out_change_0(m)                        &
                          * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)      &
                          - surf%rad_lw_out_change_0(m) * surf_t_surface_p%var_1d(m)

       surf%rad_net(m) = surf%rad_net_l(m)
       surf%rad_lw_out(m) = surf%rad_lw_out(m) + surf%rad_lw_out_change_0(m) * &
                     ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )

       surf%ghf(m) = lambda_surface * ( surf_t_surface_p%var_1d(m)             &
                                             - surf_t_soil%var_2d(nzb_soil,m) )

       surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / c_p

!
! update the 3d field of rad_lw_out array to have consistent output
       IF ( horizontal ) THEN
          IF ( radiation_scheme == 'rrtmg' ) THEN
             rad_lw_out(k+k_off,j+j_off,i+i_off) = surf%rad_lw_out(m)
          ELSE
             rad_lw_out(0,j+j_off,i+i_off) = surf%rad_lw_out(m)
          ENDIF
       ENDIF

       IF ( humidity )  THEN
          surf%qsws(m)  = - f_qsws * ( surf%qv1(m) - q_s + dq_s_dt             &
                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
                            surf_t_surface_p%var_1d(m) )

          surf%qsws_veg(m)  = - f_qsws_veg  * ( surf%qv1(m) - q_s              &
                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
                              * surf_t_surface_p%var_1d(m) )

          surf%qsws_soil(m) = - f_qsws_soil * ( surf%qv1(m) - q_s              &
                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
                              * surf_t_surface_p%var_1d(m) )

          surf%qsws_liq(m)  = - f_qsws_liq  * ( surf%qv1(m) - q_s              &
                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
                              * surf_t_surface_p%var_1d(m) )
       ENDIF

!
!--    Calculate the true surface resistance. ABS is used here to avoid negative
!--    values that can occur for very small fluxes due to the artifical addition
!--    of 1.0E-20.
       IF ( .NOT.  humidity )  THEN
          surf%r_s(m) = 1.0E10_wp
       ELSE
          surf%r_s(m) = ABS(rho_lv / (f_qsws + 1.0E-20_wp) - surf%r_a(m))
       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
!--       and qsws_soil 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 explicitly removed 
!--          (as fictive runoff by drainage systems)
             IF ( surf%pavement_surface(m) )  THEN
                IF ( surf_m_liq%var_1d(m) < m_liq_max )  THEN
                   surf%qsws_liq(m) = surf%qsws_liq(m)                         &
                                 + prr(k+k_off,j+j_off,i+i_off)                &
                                 * hyrho(k+k_off)                              &
                                 * 0.001_wp * rho_l * l_v
                ENDIF         
             ELSE
                IF ( surf_m_liq%var_1d(m) < m_liq_max )  THEN
                   surf%qsws_liq(m) = surf%qsws_liq(m)                         &
                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
                                 * hyrho(k+k_off)                              &
                                 * 0.001_wp * rho_l * l_v
                   surf%qsws_soil(m) = surf%qsws_soil(m) + ( 1.0_wp -          &
                                 surf%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
                                 * hyrho(k+k_off)                              &
                                 * 0.001_wp * rho_l * l_v                                 
                ELSE

!--                Add precipitation to bare soil according to the bare soil
!--                coverage.
                   surf%qsws_soil(m) = surf%qsws_soil(m)                       &
                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
                                 * hyrho(k+k_off)                              &
                                 * 0.001_wp * rho_l * l_v

                ENDIF
             ENDIF

          ENDIF

!
!--       If the air is saturated, check the reservoir water level
          IF ( surf%qsws(m) < 0.0_wp )  THEN
!
!--          Check if reservoir is full (avoid values > m_liq_max)
!--          In that case, qsws_liq goes to qsws_soil for pervious surfaces. In
!--          this case qsws_veg is zero anyway (because c_liq = 1),       
!--          so that tend is zero and no further check is needed
             IF ( surf_m_liq%var_1d(m) == m_liq_max )  THEN
                IF ( .NOT. surf%pavement_surface(m))  THEN
                   surf%qsws_soil(m) = surf%qsws_soil(m) + surf%qsws_liq(m)
                ENDIF
                surf%qsws_liq(m)  = 0.0_wp
             ENDIF

!
!--          In case qsws_veg becomes negative (unphysical behavior), 
!--          let the water enter the liquid water reservoir as dew on the
!--          plant
             IF ( surf%qsws_veg(m) < 0.0_wp )  THEN
                surf%qsws_liq(m) = surf%qsws_liq(m) + surf%qsws_veg(m)
                surf%qsws_veg(m) = 0.0_wp
             ENDIF
          ENDIF                    
 
          surf%qsws(m) = surf%qsws(m) / l_v
 
          tend = - surf%qsws_liq(m) * drho_l_lv
          surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d *              &
                                        ( tsc(2) * tend +                      &
                                          tsc(3) * surf_tm_liq_m%var_1d(m) )
!
!--       Check if reservoir is overfull -> reduce to maximum
!--       (conservation of water is violated here)
          surf_m_liq_p%var_1d(m) = MIN( surf_m_liq_p%var_1d(m),m_liq_max )

!
!--       Check if reservoir is empty (avoid values < 0.0)
!--       (conservation of water is violated here)
          surf_m_liq_p%var_1d(m) = MAX( surf_m_liq_p%var_1d(m), 0.0_wp )
!
!--       Calculate m_liq tendencies for the next Runge-Kutta step
          IF ( timestep_scheme(1:5) == 'runge' )  THEN
             IF ( intermediate_timestep_count == 1 )  THEN
                surf_tm_liq_m%var_1d(m) = tend
             ELSEIF ( intermediate_timestep_count <                            &
                      intermediate_timestep_count_max )  THEN
                surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend +                  &
                                           5.3125_wp * surf_tm_liq_m%var_1d(m)
             ENDIF
          ENDIF

       ENDIF

    ENDDO
    !$OMP END PARALLEL

!
!-- 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 water vapor mixing ratio
    IF ( humidity )  THEN
       CALL calc_q_surface
    ENDIF

!
!-- Calculate new roughness lengths (for water surfaces only)
    IF ( horizontal  .AND.  .NOT. constant_roughness )  CALL calc_z0_water_surface
    
    IF ( debug_output_timestep )  THEN
       WRITE( debug_string, * ) 'lsm_energy_balance', horizontal, l
       CALL debug_message( debug_string, 'end' )
    ENDIF

    CONTAINS
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculation of mixing ratio of the skin layer (surface). It is assumend
!> that the skin is always saturated.
!------------------------------------------------------------------------------!
    SUBROUTINE calc_q_surface

       IMPLICIT NONE

       REAL(wp) ::  e_s           !< saturation water vapor pressure
       REAL(wp) ::  q_s           !< saturation mixing ratio
       REAL(wp) ::  resistance    !< aerodynamic and soil resistance term


       !$OMP PARALLEL PRIVATE (m, i, j, k, e_s, q_s, resistance)
       !$OMP DO SCHEDULE (STATIC)
       DO  m = 1, surf%ns

          i   = surf%i(m)            
          j   = surf%j(m)
          k   = surf%k(m)
!
!--       Calculate water vapour pressure at saturation and convert to hPa
          e_s = 0.01_wp * magnus( MIN(surf_t_surface_p%var_1d(m), 333.15_wp) )                   

!
!--       Calculate mixing ratio at saturation
          q_s = rd_d_rv * e_s / ( surface_pressure - e_s )

          resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) + 1E-5_wp )

!
!--       Calculate mixing ratio at surface
          IF ( bulk_cloud_model )  THEN
             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
                                        ( 1.0_wp - resistance ) *              &
                                        ( q(k,j,i) - ql(k,j,i) )
          ELSE
             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
                                        ( 1.0_wp - resistance ) *              &
                                          q(k,j,i)
          ENDIF
          
          surf%q_surface(m) = q(k+k_off,j+j_off,i+i_off)
!
!--       Update virtual potential temperature
          surf%vpt_surface(m) = surf%pt_surface(m) *                           &
                                  ( 1.0_wp + 0.61_wp * surf%q_surface(m) )

        
                     
       ENDDO
       !$OMP END PARALLEL
  
    END SUBROUTINE calc_q_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

       IF ( vegetation_type_f%from_file )  THEN
          WRITE( io, 5 )
       ELSE
          WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ),        &
                         TRIM (soil_type_name(soil_type) )
       ENDIF
       WRITE( io, 6 ) 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 ('    --> Land surface type  : read from file' /                    &
            '    --> Soil porosity type : read from file' )
6   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
            '       Height:        ',A,'  m'/                                  &
            '       Temperature:   ',A,'  K'/                                  &
            '       Moisture:      ',A,'  m**3/m**3'/                          &
            '       Root fraction: ',A,'  '/                                   &
            '       Grid point:    ',A)


    END SUBROUTINE lsm_header


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Initialization of the land surface model
!------------------------------------------------------------------------------!
    SUBROUTINE lsm_init

       USE control_parameters,                                                 &
           ONLY:  message_string

       USE indices,                                                            &
           ONLY:  nx, ny, topo_min_level

       USE pmc_interface,                                                      &
           ONLY:  nested_run
    
       IMPLICIT NONE

       LOGICAL      ::  init_msoil_from_parent   !< flag controlling initialization of soil moisture in nested child domains
       LOGICAL      ::  init_tsoil_from_parent   !< flag controlling initialization of soil temperature in nested child domains

       INTEGER(iwp) ::  i                       !< running index
       INTEGER(iwp) ::  j                       !< running index
       INTEGER(iwp) ::  k                       !< running index
       INTEGER(iwp) ::  kn                      !< running index
       INTEGER(iwp) ::  ko                      !< running index
       INTEGER(iwp) ::  kroot                   !< running index
       INTEGER(iwp) ::  kzs                     !< running index
       INTEGER(iwp) ::  l                       !< running index surface facing
       INTEGER(iwp) ::  m                       !< running index
       INTEGER(iwp) ::  st                      !< soil-type index
       INTEGER(iwp) ::  n_soil_layers_total     !< temperature variable, stores the total number of soil layers + 4 

       REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_soil_init !< temporary array used for averaging soil profiles

       IF ( debug_output )  CALL debug_message( 'lsm_init', 'start' )
!
!--    If no cloud physics is used, rho_surface has not been calculated before
       IF (  .NOT.  bulk_cloud_model  .AND.  .NOT.  cloud_droplets )  THEN
          CALL calc_mean_profile( pt, 4 )
          rho_surface = hyp(nzb) / ( r_d * hom(topo_min_level+1,1,4,0) * exner(nzb) )
       ENDIF

!
!--    Calculate frequently used parameters
       rho_cp    = c_p * rho_surface
       rho_lv    = rho_surface * l_v
       drho_l_lv = 1.0_wp / (rho_l * l_v)

!
!--    Set initial values for prognostic quantities
!--    Horizontal surfaces
       tt_surface_h_m%var_1d = 0.0_wp
       tt_soil_h_m%var_2d    = 0.0_wp
       tm_soil_h_m%var_2d    = 0.0_wp
       tm_liq_h_m%var_1d     = 0.0_wp
       surf_lsm_h%c_liq      = 0.0_wp

       surf_lsm_h%ghf = 0.0_wp

       surf_lsm_h%qsws_liq  = 0.0_wp
       surf_lsm_h%qsws_soil = 0.0_wp
       surf_lsm_h%qsws_veg  = 0.0_wp

       surf_lsm_h%r_a        = 50.0_wp
       surf_lsm_h%r_s        = 50.0_wp
       surf_lsm_h%r_canopy   = 0.0_wp
       surf_lsm_h%r_soil     = 0.0_wp
!
!--    Do the same for vertical surfaces
       DO  l = 0, 3
          tt_surface_v_m(l)%var_1d = 0.0_wp
          tt_soil_v_m(l)%var_2d    = 0.0_wp
          tm_soil_v_m(l)%var_2d    = 0.0_wp
          tm_liq_v_m(l)%var_1d     = 0.0_wp
          surf_lsm_v(l)%c_liq      = 0.0_wp

          surf_lsm_v(l)%ghf = 0.0_wp

          surf_lsm_v(l)%qsws_liq  = 0.0_wp
          surf_lsm_v(l)%qsws_soil = 0.0_wp
          surf_lsm_v(l)%qsws_veg  = 0.0_wp

          surf_lsm_v(l)%r_a        = 50.0_wp
          surf_lsm_v(l)%r_s        = 50.0_wp
          surf_lsm_v(l)%r_canopy   = 0.0_wp
          surf_lsm_v(l)%r_soil     = 0.0_wp
       ENDDO

!
!--    Set initial values for prognostic soil quantities
       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
          t_soil_h%var_2d = 0.0_wp
          m_soil_h%var_2d = 0.0_wp
          m_liq_h%var_1d  = 0.0_wp

          DO  l = 0, 3
             t_soil_v(l)%var_2d = 0.0_wp
             m_soil_v(l)%var_2d = 0.0_wp
             m_liq_v(l)%var_1d  = 0.0_wp
          ENDDO
       ENDIF
!
!--    Allocate 3D soil model arrays
!--    First, for horizontal surfaces
       ALLOCATE ( surf_lsm_h%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
       ALLOCATE ( surf_lsm_h%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
       ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
       ALLOCATE ( surf_lsm_h%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns))
       ALLOCATE ( surf_lsm_h%l_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
       ALLOCATE ( surf_lsm_h%m_fc(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
       ALLOCATE ( surf_lsm_h%m_res(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
       ALLOCATE ( surf_lsm_h%m_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
       ALLOCATE ( surf_lsm_h%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_h%ns)      )
       ALLOCATE ( surf_lsm_h%n_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
       ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
       ALLOCATE ( surf_lsm_h%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
       ALLOCATE ( surf_lsm_h%root_fr(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
    
       surf_lsm_h%lambda_h     = 0.0_wp
!
!--    If required, allocate humidity-related variables for the soil model
       IF ( humidity )  THEN
          ALLOCATE ( surf_lsm_h%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
          ALLOCATE ( surf_lsm_h%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )  

          surf_lsm_h%lambda_w = 0.0_wp 
       ENDIF
!
!--    For vertical surfaces
       DO  l = 0, 3
          ALLOCATE ( surf_lsm_v(l)%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
          ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
          ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
          ALLOCATE ( surf_lsm_v(l)%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns))
          ALLOCATE ( surf_lsm_v(l)%l_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
          ALLOCATE ( surf_lsm_v(l)%m_fc(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
          ALLOCATE ( surf_lsm_v(l)%m_res(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
          ALLOCATE ( surf_lsm_v(l)%m_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
          ALLOCATE ( surf_lsm_v(l)%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)      )
          ALLOCATE ( surf_lsm_v(l)%n_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
          ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )  
          ALLOCATE ( surf_lsm_v(l)%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )  
          ALLOCATE ( surf_lsm_v(l)%root_fr(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )

          surf_lsm_v(l)%lambda_h     = 0.0_wp 
          
!
!--       If required, allocate humidity-related variables for the soil model
          IF ( humidity )  THEN
             ALLOCATE ( surf_lsm_v(l)%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
             ALLOCATE ( surf_lsm_v(l)%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )  

             surf_lsm_v(l)%lambda_w = 0.0_wp 
          ENDIF     
       ENDDO
!
!--    Allocate albedo type and emissivity for vegetation, water and pavement 
!--    fraction.
!--    Set default values at each surface element. 
       ALLOCATE ( surf_lsm_h%albedo_type(0:2,1:surf_lsm_h%ns) )
       ALLOCATE ( surf_lsm_h%emissivity(0:2,1:surf_lsm_h%ns) )
!
!--    Initialize albedo type according to its default type, in order to set values
!--    independent on default albedo_type in radiation model. 
       surf_lsm_h%albedo_type(ind_veg_wall,:)  =                               &
                             INT( vegetation_pars(ind_v_at,vegetation_type) )
       surf_lsm_h%albedo_type(ind_wat_win,:)   =                               &
                             INT( water_pars(ind_w_at,water_type)           )
       surf_lsm_h%albedo_type(ind_pav_green,:) =                               &
                             INT( pavement_pars(ind_p_at,pavement_type)     )
       surf_lsm_h%emissivity  = emissivity
       DO  l = 0, 3
          ALLOCATE ( surf_lsm_v(l)%albedo_type(0:2,1:surf_lsm_v(l)%ns) )
          ALLOCATE ( surf_lsm_v(l)%emissivity(0:2,1:surf_lsm_v(l)%ns)  )
!
!--       Initialize albedo type according to its default type, in order to 
!--       set values independent on default albedo_type in radiation model. 
          surf_lsm_v(l)%albedo_type(ind_veg_wall,:)  =                         &
                             INT( vegetation_pars(ind_v_at,vegetation_type) )
          surf_lsm_v(l)%albedo_type(ind_wat_win,:)   =                         &
                             INT( water_pars(ind_w_at,water_type)           )
          surf_lsm_v(l)%albedo_type(ind_pav_green,:) =                         &
                             INT( pavement_pars(ind_p_at,pavement_type)     )
          surf_lsm_v(l)%emissivity  = emissivity
       ENDDO
!
!--    Allocate arrays for relative surface fraction. 
!--    0 - vegetation fraction, 2 - water fraction, 1 - pavement fraction
       ALLOCATE( surf_lsm_h%frac(0:2,1:surf_lsm_h%ns) )
       surf_lsm_h%frac = 0.0_wp
       DO  l = 0, 3
          ALLOCATE( surf_lsm_v(l)%frac(0:2,1:surf_lsm_v(l)%ns) )
          surf_lsm_v(l)%frac = 0.0_wp
       ENDDO
!
!--    For vertical walls only - allocate special flag indicating if any building is on
!--    top of any natural surfaces. Used for initialization only. 
       DO  l = 0, 3
          ALLOCATE( surf_lsm_v(l)%building_covered(1:surf_lsm_v(l)%ns) )
       ENDDO
!
!--    Allocate arrays for the respective types and their names on the surface
!--    elements. This will be required to treat deposition of chemical species.
       ALLOCATE( surf_lsm_h%pavement_type(1:surf_lsm_h%ns)   )
       ALLOCATE( surf_lsm_h%vegetation_type(1:surf_lsm_h%ns) )
       ALLOCATE( surf_lsm_h%water_type(1:surf_lsm_h%ns)      )
       
       surf_lsm_h%pavement_type   = 0
       surf_lsm_h%vegetation_type = 0
       surf_lsm_h%water_type      = 0
       
       ALLOCATE( surf_lsm_h%pavement_type_name(1:surf_lsm_h%ns)   )
       ALLOCATE( surf_lsm_h%vegetation_type_name(1:surf_lsm_h%ns) )
       ALLOCATE( surf_lsm_h%water_type_name(1:surf_lsm_h%ns)      )
       
       surf_lsm_h%pavement_type_name   = 'none'
       surf_lsm_h%vegetation_type_name = 'none'
       surf_lsm_h%water_type_name      = 'none'
       
       DO  l = 0, 3
          ALLOCATE( surf_lsm_v(l)%pavement_type(1:surf_lsm_v(l)%ns)   )
          ALLOCATE( surf_lsm_v(l)%vegetation_type(1:surf_lsm_v(l)%ns) )
          ALLOCATE( surf_lsm_v(l)%water_type(1:surf_lsm_v(l)%ns)      )
          
          surf_lsm_v(l)%pavement_type   = 0
          surf_lsm_v(l)%vegetation_type = 0
          surf_lsm_v(l)%water_type      = 0
       
          ALLOCATE( surf_lsm_v(l)%pavement_type_name(1:surf_lsm_v(l)%ns)   )
          ALLOCATE( surf_lsm_v(l)%vegetation_type_name(1:surf_lsm_v(l)%ns) )
          ALLOCATE( surf_lsm_v(l)%water_type_name(1:surf_lsm_v(l)%ns)      )
       
          surf_lsm_v(l)%pavement_type_name   = 'none'
          surf_lsm_v(l)%vegetation_type_name = 'none'
          surf_lsm_v(l)%water_type_name      = 'none'       
       ENDDO
       
!
!--    Set flag parameter for the prescribed surface type depending on user 
!--    input. Set surface fraction to 1 for the respective type.
       SELECT CASE ( TRIM( surface_type ) )
          
          CASE ( 'vegetation' )
         
             surf_lsm_h%vegetation_surface = .TRUE.
             surf_lsm_h%frac(ind_veg_wall,:) = 1.0_wp
             DO  l = 0, 3
                surf_lsm_v(l)%vegetation_surface = .TRUE.
                surf_lsm_v(l)%frac(ind_veg_wall,:) = 1.0_wp
             ENDDO
   
          CASE ( 'water' )
             
             surf_lsm_h%water_surface = .TRUE.
             surf_lsm_h%frac(ind_wat_win,:) = 1.0_wp
!
!--          Note, vertical water surface does not really make sense.
             DO  l = 0, 3  
                surf_lsm_v(l)%water_surface   = .TRUE.
                surf_lsm_v(l)%frac(ind_wat_win,:) = 1.0_wp
             ENDDO

          CASE ( 'pavement' )
             
             surf_lsm_h%pavement_surface = .TRUE.
                surf_lsm_h%frac(ind_pav_green,:) = 1.0_wp
             DO  l = 0, 3
                surf_lsm_v(l)%pavement_surface   = .TRUE.
                surf_lsm_v(l)%frac(ind_pav_green,:) = 1.0_wp
             ENDDO

          CASE ( 'netcdf' )

             DO  m = 1, surf_lsm_h%ns
                i = surf_lsm_h%i(m)
                j = surf_lsm_h%j(m)
                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )    &
                   surf_lsm_h%vegetation_surface(m) = .TRUE.
                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )      &
                   surf_lsm_h%pavement_surface(m) = .TRUE.
                IF ( water_type_f%var(j,i)      /= water_type_f%fill )         &
                   surf_lsm_h%water_surface(m) = .TRUE.
!
!--             Check if at least one type is set.
                IF ( .NOT. surf_lsm_h%vegetation_surface(m)  .AND.             &
                     .NOT. surf_lsm_h%pavement_surface(m)    .AND.             &
                     .NOT. surf_lsm_h%water_surface(m) )  THEN
                   WRITE( message_string, * ) 'Horizontal surface element ' // &
                                       ' at i, j = ',  i, j,                   &
                                       ' is neither a vegetation, ' //         &
                                       'pavement, nor a water surface.'
                   CALL message( 'land_surface_model_mod', 'PA0619',          &
                                  2, 2, myid, 6, 0 )
                ENDIF
                
             ENDDO
!
!--          For vertical surfaces some special checks and treatment are
!--          required for correct initialization. 
             DO  l = 0, 3
                DO  m = 1, surf_lsm_v(l)%ns
!
!--                Only for vertical surfaces. Check if at the grid point where 
!--                the wall is defined (i+ioff, j+joff) is any building. 
!--                This case, no natural surfaces properties will be defined at
!--                at this grid point, leading to problems in the initialization.
!--                To overcome this, define a special flag which
!--                indicates that a building is defined at the wall grid point  
!--                and take the surface properties from the adjoining grid 
!--                point, i.e. without offset values.
!--                Further, there can occur a special case where elevation
!--                changes are larger than building heights. This case, (j,i)
!--                and (j+joff,i+ioff) grid points may be both covered by 
!--                buildings, but vertical, but vertically natural walls may 
!--                be located between the buildings. This case, it is not 
!--                guaranteed that information about natural surface types is 
!--                given, neither at (j,i) nor at (j+joff,i+ioff), again leading 
!--                to non-initialized surface properties. 
                   surf_lsm_v(l)%building_covered(m) = .FALSE.
!
!--                Wall grid point is building-covered. This case, set 
!--                flag indicating that surface properties are initialized 
!--                from neighboring reference grid point, which is not 
!--                building_covered. 
                   IF ( building_type_f%from_file )  THEN
                      i = surf_lsm_v(l)%i(m)
                      j = surf_lsm_v(l)%j(m)
                      IF ( building_type_f%var(j+surf_lsm_v(l)%joff,           &
                                               i+surf_lsm_v(l)%ioff) /=        &
                           building_type_f%fill )                              &
                         surf_lsm_v(l)%building_covered(m) = .TRUE.
!
!--                   Wall grid point as well as neighboring reference grid
!--                   point are both building-covered. This case, surface
!--                   properties are not necessarily defined (not covered by
!--                   checks for static input file) at this surface. Hence, 
!--                   initialize surface properties by simply setting 
!--                   vegetation_type_f to bare-soil bulk parametrization. 
!--                   soil_type_f as well as surface_fractions_f will be set 
!--                   also.                      
                      IF ( building_type_f%var(j+surf_lsm_v(l)%joff,           &
                                               i+surf_lsm_v(l)%ioff) /=        &
                           building_type_f%fill  .AND.                         &
                           building_type_f%var(j,i) /= building_type_f%fill )  &
                      THEN
                         vegetation_type_f%var(j,i)                 = 1 ! bare soil
                         soil_type_f%var_2d(j,i)                    = 1 
!                         
!--                      Set surface_fraction if provided in static input, 
!--                      else, in case no tiles are used, this will be done
!--                      on basis of the prescribed types (vegetation/pavement/
!--                      water_type). 
                         IF ( surface_fraction_f%from_file )  THEN
                            surface_fraction_f%frac(ind_veg_wall,j,i)  = 1.0_wp
                            surface_fraction_f%frac(ind_pav_green,j,i) = 0.0_wp
                            surface_fraction_f%frac(ind_wat_win,j,i)   = 0.0_wp
                         ENDIF
                      ENDIF
                      
                   ENDIF
!
!--                Normally proceed with setting surface types.
                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
                                            surf_lsm_v(l)%building_covered(m) )
                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
                                            surf_lsm_v(l)%building_covered(m) )
                   IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
                      surf_lsm_v(l)%vegetation_surface(m) = .TRUE.
                   IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )   &
                      surf_lsm_v(l)%pavement_surface(m) = .TRUE.
                   IF ( water_type_f%var(j,i)      /= water_type_f%fill )      &
                      surf_lsm_v(l)%water_surface(m) = .TRUE.
!
!--                Check if at least one type is set.
                   IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m)  .AND.       &
                        .NOT. surf_lsm_v(l)%pavement_surface(m)    .AND.       &
                        .NOT. surf_lsm_v(l)%water_surface(m) )  THEN
                      WRITE( message_string, * ) 'Vertical surface element ' //&
                                       ' at i, j = ',  i, j,                   &
                                       ' is neither a vegetation, ' //         &
                                       'pavement, nor a water surface.'
                      CALL message( 'land_surface_model_mod', 'PA0619',        &
                                     2, 2, myid, 6, 0 )
                   ENDIF
                ENDDO
             ENDDO

       END SELECT
!
!--    In case of netcdf input file, further initialize surface fractions. 
!--    At the moment only 1 surface is given at a location, so that the fraction
!--    is either 0 or 1. This will be revised later. If surface fraction 
!--    is not given in static input file, relative fractions will be derived
!--    from given surface type. In this case, only 1 type is given at a certain
!--    location (already checked).  
       IF ( input_pids_static  .AND.  surface_fraction_f%from_file )  THEN
          DO  m = 1, surf_lsm_h%ns
             i = surf_lsm_h%i(m)
             j = surf_lsm_h%j(m)
!
!--          0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction             
             surf_lsm_h%frac(ind_veg_wall,m)  =                                &
                                    surface_fraction_f%frac(ind_veg_wall,j,i)         
             surf_lsm_h%frac(ind_pav_green,m) =                                &
                                    surface_fraction_f%frac(ind_pav_green,j,i)        
             surf_lsm_h%frac(ind_wat_win,m)   =                                &
                                    surface_fraction_f%frac(ind_wat_win,j,i)

          ENDDO
          DO  l = 0, 3
             DO  m = 1, surf_lsm_v(l)%ns
                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
                                                surf_lsm_v(l)%building_covered(m) ) 
                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
                                                surf_lsm_v(l)%building_covered(m) ) 
!
!--             0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction        
                surf_lsm_v(l)%frac(ind_veg_wall,m)  =                          &
                                    surface_fraction_f%frac(ind_veg_wall,j,i)         
                surf_lsm_v(l)%frac(ind_pav_green,m) =                          &
                                    surface_fraction_f%frac(ind_pav_green,j,i)        
                surf_lsm_v(l)%frac(ind_wat_win,m)   =                          &
                                    surface_fraction_f%frac(ind_wat_win,j,i)

             ENDDO
          ENDDO
       ELSEIF ( input_pids_static  .AND.  .NOT. surface_fraction_f%from_file ) &
       THEN
          DO  m = 1, surf_lsm_h%ns
             i = surf_lsm_h%i(m)
             j = surf_lsm_h%j(m)

             IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )       &        
                surf_lsm_h%frac(ind_veg_wall,m)  = 1.0_wp
             IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill   )       &        
                surf_lsm_h%frac(ind_pav_green,m) = 1.0_wp 
             IF ( water_type_f%var(j,i)      /= water_type_f%fill      )       &        
                surf_lsm_h%frac(ind_wat_win,m)   = 1.0_wp       
          ENDDO
          DO  l = 0, 3
             DO  m = 1, surf_lsm_v(l)%ns
                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
                                                surf_lsm_v(l)%building_covered(m) ) 
                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
                                                surf_lsm_v(l)%building_covered(m) ) 
     
                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )    &        
                   surf_lsm_v(l)%frac(ind_veg_wall,m)  = 1.0_wp
                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill   )    &        
                   surf_lsm_v(l)%frac(ind_pav_green,m) = 1.0_wp 
                IF ( water_type_f%var(j,i)      /= water_type_f%fill      )    &        
                   surf_lsm_v(l)%frac(ind_wat_win,m)   = 1.0_wp     
             ENDDO
          ENDDO
       ENDIF
!
!--    Level 1, initialization of soil parameters.
!--    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 = soil_pars(4,soil_type)           
          ENDIF

          IF ( field_capacity == 9999999.9_wp )  THEN
             field_capacity = soil_pars(5,soil_type)           
          ENDIF

          IF ( wilting_point == 9999999.9_wp )  THEN
             wilting_point = soil_pars(6,soil_type)            
          ENDIF

          IF ( residual_moisture == 9999999.9_wp )  THEN
             residual_moisture = soil_pars(7,soil_type)       
          ENDIF

       ENDIF
!
!--    Map values to the respective 2D/3D arrays
!--    Horizontal surfaces
       surf_lsm_h%alpha_vg      = alpha_vangenuchten
       surf_lsm_h%l_vg          = l_vangenuchten
       surf_lsm_h%n_vg          = n_vangenuchten 
       surf_lsm_h%gamma_w_sat   = hydraulic_conductivity
       surf_lsm_h%m_sat         = saturation_moisture
       surf_lsm_h%m_fc          = field_capacity
       surf_lsm_h%m_wilt        = wilting_point
       surf_lsm_h%m_res         = residual_moisture
       surf_lsm_h%r_soil_min    = min_soil_resistance
!
!--    Vertical surfaces
       DO  l = 0, 3
          surf_lsm_v(l)%alpha_vg      = alpha_vangenuchten
          surf_lsm_v(l)%l_vg          = l_vangenuchten
          surf_lsm_v(l)%n_vg          = n_vangenuchten 
          surf_lsm_v(l)%gamma_w_sat   = hydraulic_conductivity
          surf_lsm_v(l)%m_sat         = saturation_moisture
          surf_lsm_v(l)%m_fc          = field_capacity
          surf_lsm_v(l)%m_wilt        = wilting_point
          surf_lsm_v(l)%m_res         = residual_moisture
          surf_lsm_v(l)%r_soil_min    = min_soil_resistance
       ENDDO
!
!--    Level 2, initialization of soil parameters via soil_type read from file.
!--    Soil parameters are initialized for each (y,x)-grid point 
!--    individually using default paramter settings according to the given
!--    soil type.
       IF ( soil_type_f%from_file )  THEN
!
!--       Level of detail = 1, i.e. a homogeneous soil distribution along the 
!--       vertical dimension is assumed.
          IF ( soil_type_f%lod == 1 )  THEN
!
!--          Horizontal surfaces
             DO  m = 1, surf_lsm_h%ns
                i = surf_lsm_h%i(m)
                j = surf_lsm_h%j(m)
             
                st = soil_type_f%var_2d(j,i)
                IF ( st /= soil_type_f%fill )  THEN
                   surf_lsm_h%alpha_vg(:,m)    = soil_pars(0,st)
                   surf_lsm_h%l_vg(:,m)        = soil_pars(1,st)
                   surf_lsm_h%n_vg(:,m)        = soil_pars(2,st)
                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars(3,st)
                   surf_lsm_h%m_sat(:,m)       = soil_pars(4,st)
                   surf_lsm_h%m_fc(:,m)        = soil_pars(5,st)
                   surf_lsm_h%m_wilt(:,m)      = soil_pars(6,st)
                   surf_lsm_h%m_res(:,m)       = soil_pars(7,st)
                ENDIF
             ENDDO
!
!--          Vertical surfaces ( assumes the soil type given at respective (x,y) 
             DO  l = 0, 3
                DO  m = 1, surf_lsm_v(l)%ns
                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
                                                   surf_lsm_v(l)%building_covered(m) ) 
                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
                                                   surf_lsm_v(l)%building_covered(m) ) 

                   st = soil_type_f%var_2d(j,i)
                   IF ( st /= soil_type_f%fill )  THEN
                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars(0,st)
                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars(1,st)
                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars(2,st)
                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars(3,st)
                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars(4,st)
                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars(5,st)
                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars(6,st)
                      surf_lsm_v(l)%m_res(:,m)       = soil_pars(7,st)
                   ENDIF
                ENDDO
             ENDDO
!
!--       Level of detail = 2, i.e. soil type and thus the soil parameters
!--       can be heterogeneous along the vertical dimension. 
          ELSE
!
!--          Horizontal surfaces
             DO  m = 1, surf_lsm_h%ns
                i = surf_lsm_h%i(m)
                j = surf_lsm_h%j(m)
             
                DO  k = nzb_soil, nzt_soil
                   st = soil_type_f%var_3d(k,j,i)
                   IF ( st /= soil_type_f%fill )  THEN
                      surf_lsm_h%alpha_vg(k,m)    = soil_pars(0,st)
                      surf_lsm_h%l_vg(k,m)        = soil_pars(1,st)
                      surf_lsm_h%n_vg(k,m)        = soil_pars(2,st)
                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars(3,st)
                      surf_lsm_h%m_sat(k,m)       = soil_pars(4,st)
                      surf_lsm_h%m_fc(k,m)        = soil_pars(5,st)
                      surf_lsm_h%m_wilt(k,m)      = soil_pars(6,st)
                      surf_lsm_h%m_res(k,m)       = soil_pars(7,st)
                   ENDIF
                ENDDO
             ENDDO
!
!--          Vertical surfaces ( assumes the soil type given at respective (x,y) 
             DO  l = 0, 3
                DO  m = 1, surf_lsm_v(l)%ns
                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
                                                   surf_lsm_v(l)%building_covered(m) ) 
                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
                                                   surf_lsm_v(l)%building_covered(m) ) 

                   DO  k = nzb_soil, nzt_soil
                      st = soil_type_f%var_3d(k,j,i)
                      IF ( st /= soil_type_f%fill )  THEN
                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars(0,st)
                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars(1,st)
                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars(2,st)
                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars(3,st)
                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars(4,st)
                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars(5,st)
                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars(6,st)
                         surf_lsm_v(l)%m_res(k,m)       = soil_pars(7,st)
                      ENDIF
                   ENDDO
                ENDDO
             ENDDO
          ENDIF
       ENDIF
!
!--    Level 3, initialization of single soil parameters at single z,x,y 
!--    position via soil_pars read from file.
       IF ( soil_pars_f%from_file )  THEN
!
!--       Level of detail = 1, i.e. a homogeneous vertical distribution of soil 
!--       parameters is assumed.
!--       Horizontal surfaces
          IF ( soil_pars_f%lod == 1 )  THEN
!
!--          Horizontal surfaces 
             DO  m = 1, surf_lsm_h%ns
                i = surf_lsm_h%i(m)
                j = surf_lsm_h%j(m)

                IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )              &
                   surf_lsm_h%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
                IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )              &
                   surf_lsm_h%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
                IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )              &
                   surf_lsm_h%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
                IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )              &
                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
                IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )              &
                   surf_lsm_h%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
                IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )              &
                   surf_lsm_h%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
                IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )              &
                   surf_lsm_h%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
                IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )              &
                   surf_lsm_h%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)

             ENDDO
!
!--          Vertical surfaces 
             DO  l = 0, 3
                DO  m = 1, surf_lsm_v(l)%ns
                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
                                                   surf_lsm_v(l)%building_covered(m) ) 
                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
                                                   surf_lsm_v(l)%building_covered(m) ) 

                   IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )           &
                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
                   IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )           &
                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
                   IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )           &
                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
                   IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )           &
                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
                   IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )           &
                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
                   IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )           &
                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
                   IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )           &
                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
                   IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )           &
                      surf_lsm_v(l)%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)

                ENDDO
             ENDDO
!
!--       Level of detail = 2, i.e. soil parameters can be set at each soil 
!--       layer individually. 
          ELSE
!
!--          Horizontal surfaces
             DO  m = 1, surf_lsm_h%ns
                i = surf_lsm_h%i(m)
                j = surf_lsm_h%j(m)

                DO  k = nzb_soil, nzt_soil
                   IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
                      surf_lsm_h%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
                   IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
                      surf_lsm_h%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
                   IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
                      surf_lsm_h%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
                   IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
                   IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
                      surf_lsm_h%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
                   IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
                      surf_lsm_h%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
                   IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
                      surf_lsm_h%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
                   IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
                      surf_lsm_h%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
                ENDDO

             ENDDO
!
!--          Vertical surfaces 
             DO  l = 0, 3
                DO  m = 1, surf_lsm_v(l)%ns
                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
                                                   surf_lsm_v(l)%building_covered(m) ) 
                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
                                                   surf_lsm_v(l)%building_covered(m) ) 

                   DO  k = nzb_soil, nzt_soil
                      IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
                      IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
                      IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
                      IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
                      IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
                      IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
                      IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
                      IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
                         surf_lsm_v(l)%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
                   ENDDO

                ENDDO
             ENDDO

          ENDIF
       ENDIF

!
!--    Level 1, initialization of vegetation parameters. A horizontally 
!--    homogeneous distribution is assumed here. 
       IF ( vegetation_type /= 0 )  THEN

          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
             min_canopy_resistance = vegetation_pars(ind_v_rc_min,vegetation_type)
          ENDIF

          IF ( leaf_area_index == 9999999.9_wp )  THEN
             leaf_area_index = vegetation_pars(ind_v_rc_lai,vegetation_type)         
          ENDIF

          IF ( vegetation_coverage == 9999999.9_wp )  THEN
             vegetation_coverage = vegetation_pars(ind_v_c_veg,vegetation_type)     
          ENDIF

          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
              canopy_resistance_coefficient= vegetation_pars(ind_v_gd,vegetation_type)     
          ENDIF

          IF ( z0_vegetation == 9999999.9_wp )  THEN
             z0_vegetation  = vegetation_pars(ind_v_z0,vegetation_type) 
          ENDIF

          IF ( z0h_vegetation == 9999999.9_wp )  THEN
             z0h_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
          ENDIF 
          
          IF ( z0q_vegetation == 9999999.9_wp )  THEN
             z0q_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
          ENDIF
          
          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
             lambda_surface_stable = vegetation_pars(ind_v_lambda_s,vegetation_type) 
          ENDIF

          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
             lambda_surface_unstable = vegetation_pars(ind_v_lambda_u,vegetation_type)           
          ENDIF

          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
             f_shortwave_incoming = vegetation_pars(ind_v_f_sw_in,vegetation_type)        
          ENDIF

          IF ( c_surface == 9999999.9_wp )  THEN
             c_surface = vegetation_pars(ind_v_c_surf,vegetation_type)        
          ENDIF

          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
             albedo_type = INT(vegetation_pars(ind_v_at,vegetation_type))       
          ENDIF
   
          IF ( emissivity == 9999999.9_wp )  THEN
             emissivity = vegetation_pars(ind_v_emis,vegetation_type)     
          ENDIF

       ENDIF
!
!--    Map values onto horizontal elemements
       DO  m = 1, surf_lsm_h%ns
          IF ( surf_lsm_h%vegetation_surface(m) )  THEN
             surf_lsm_h%r_canopy_min(m)     = min_canopy_resistance
             surf_lsm_h%lai(m)              = leaf_area_index
             surf_lsm_h%c_veg(m)            = vegetation_coverage
             surf_lsm_h%g_d(m)              = canopy_resistance_coefficient
             surf_lsm_h%z0(m)               = z0_vegetation
             surf_lsm_h%z0h(m)              = z0h_vegetation
             surf_lsm_h%z0q(m)              = z0q_vegetation
             surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
             surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
             surf_lsm_h%f_sw_in(m)          = f_shortwave_incoming
             surf_lsm_h%c_surface(m)        = c_surface
             surf_lsm_h%albedo_type(ind_veg_wall,m) = albedo_type
             surf_lsm_h%emissivity(ind_veg_wall,m)  = emissivity
             
             surf_lsm_h%vegetation_type(m)      = vegetation_type
             surf_lsm_h%vegetation_type_name(m) = vegetation_type_name(vegetation_type)
          ELSE
             surf_lsm_h%lai(m)   = 0.0_wp
             surf_lsm_h%c_veg(m) = 0.0_wp
             surf_lsm_h%g_d(m)   = 0.0_wp
          ENDIF
  
       ENDDO
!
!--    Map values onto vertical elements, even though this does not make
!--    much sense.
       DO  l = 0, 3
          DO  m = 1, surf_lsm_v(l)%ns
             IF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
                surf_lsm_v(l)%r_canopy_min(m)     = min_canopy_resistance
                surf_lsm_v(l)%lai(m)              = leaf_area_index
                surf_lsm_v(l)%c_veg(m)            = vegetation_coverage
                surf_lsm_v(l)%g_d(m)              = canopy_resistance_coefficient
                surf_lsm_v(l)%z0(m)               = z0_vegetation
                surf_lsm_v(l)%z0h(m)              = z0h_vegetation
                surf_lsm_v(l)%z0q(m)              = z0q_vegetation
                surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
                surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
                surf_lsm_v(l)%f_sw_in(m)          = f_shortwave_incoming
                surf_lsm_v(l)%c_surface(m)        = c_surface
                surf_lsm_v(l)%albedo_type(ind_veg_wall,m) = albedo_type
                surf_lsm_v(l)%emissivity(ind_veg_wall,m)  = emissivity
                
                surf_lsm_v(l)%vegetation_type(m)      = vegetation_type
                surf_lsm_v(l)%vegetation_type_name(m) = vegetation_type_name(vegetation_type)
             ELSE
                surf_lsm_v(l)%lai(m)   = 0.0_wp
                surf_lsm_v(l)%c_veg(m) = 0.0_wp
                surf_lsm_v(l)%g_d(m)   = 0.0_wp
             ENDIF
          ENDDO
       ENDDO

!
!--    Level 2, initialization of vegation parameters via vegetation_type read
!--    from file. Vegetation parameters are initialized for each (y,x)-grid point 
!--    individually using default paramter settings according to the given
!--    vegetation type.
       IF ( vegetation_type_f%from_file )  THEN
!
!--       Horizontal surfaces
          DO  m = 1, surf_lsm_h%ns
             i = surf_lsm_h%i(m)
             j = surf_lsm_h%j(m)
             
             st = vegetation_type_f%var(j,i)
             IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
                surf_lsm_h%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
                surf_lsm_h%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
                surf_lsm_h%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
                surf_lsm_h%g_d(m)              = vegetation_pars(ind_v_gd,st)
                surf_lsm_h%z0(m)               = vegetation_pars(ind_v_z0,st)
                surf_lsm_h%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
                surf_lsm_h%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
                surf_lsm_h%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
                surf_lsm_h%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
                surf_lsm_h%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
                surf_lsm_h%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
                surf_lsm_h%albedo_type(ind_veg_wall,m) = INT( vegetation_pars(ind_v_at,st) )
                surf_lsm_h%emissivity(ind_veg_wall,m)  = vegetation_pars(ind_v_emis,st)
                
                surf_lsm_h%vegetation_type(m)      = st
                surf_lsm_h%vegetation_type_name(m) = vegetation_type_name(st)
             ENDIF
          ENDDO
!
!--       Vertical surfaces
          DO  l = 0, 3
             DO  m = 1, surf_lsm_v(l)%ns
                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
                                                surf_lsm_v(l)%building_covered(m) ) 
                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
                                                surf_lsm_v(l)%building_covered(m) ) 
             
                st = vegetation_type_f%var(j,i)
                IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
                   surf_lsm_v(l)%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
                   surf_lsm_v(l)%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
                   surf_lsm_v(l)%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
                   surf_lsm_v(l)%g_d(m)              = vegetation_pars(ind_v_gd,st)
                   surf_lsm_v(l)%z0(m)               = vegetation_pars(ind_v_z0,st)
                   surf_lsm_v(l)%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
                   surf_lsm_v(l)%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
                   surf_lsm_v(l)%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
                   surf_lsm_v(l)%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
                   surf_lsm_v(l)%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
                   surf_lsm_v(l)%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
                   surf_lsm_v(l)%albedo_type(ind_veg_wall,m) = INT( vegetation_pars(ind_v_at,st) )
                   surf_lsm_v(l)%emissivity(ind_veg_wall,m)  = vegetation_pars(ind_v_emis,st)
                   
                   surf_lsm_v(l)%vegetation_type(m)      = st
                   surf_lsm_v(l)%vegetation_type_name(m) = vegetation_type_name(st)
                ENDIF
             ENDDO
          ENDDO
       ENDIF
!
!--    Level 3, initialization of vegation parameters at single (x,y) 
!--    position via vegetation_pars read from file.
       IF ( vegetation_pars_f%from_file )  THEN
!
!--       Horizontal surfaces 
          DO  m = 1, surf_lsm_h%ns

             i = surf_lsm_h%i(m)
             j = surf_lsm_h%j(m)
!
!--          If surface element is not a vegetation surface and any value in 
!--          vegetation_pars is given, neglect this information and give an
!--          informative message that this value will not be used.   
             IF ( .NOT. surf_lsm_h%vegetation_surface(m)  .AND.                &
                   ANY( vegetation_pars_f%pars_xy(:,j,i) /=                    &
                   vegetation_pars_f%fill ) )  THEN
                WRITE( message_string, * )                                     &
                                 'surface element at grid point (j,i) = (',    &
                                 j, i, ') is not a vegation surface, ',        &
                                 'so that information given in ',              &
                                 'vegetation_pars at this point is neglected.' 
                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
             ELSE

                IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=            &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%r_canopy_min(m)  =                               &
                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
                IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=            &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%lai(m)           =                               &
                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
                IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=             &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%c_veg(m)         =                               &
                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
                IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=                &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%g_d(m)           =                               &
                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
                IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=                &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%z0(m)            =                               &
                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
                IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=              &
                     vegetation_pars_f%fill )  THEN
                   surf_lsm_h%z0h(m)           =                               &
                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
                   surf_lsm_h%z0q(m)           =                               &
                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
                ENDIF
                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=          &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%lambda_surface_s(m) =                            &
                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=          &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%lambda_surface_u(m) =                            &
                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
                IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=           &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%f_sw_in(m)          =                            &
                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
                IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=            &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%c_surface(m)        =                            &
                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
                IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=                &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%albedo_type(ind_veg_wall,m) =                    &
                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
                IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=              &
                     vegetation_pars_f%fill )                                  &
                   surf_lsm_h%emissivity(ind_veg_wall,m)  =                    &
                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
             ENDIF
          ENDDO
!
!--       Vertical surfaces 
          DO  l = 0, 3
             DO  m = 1, surf_lsm_v(l)%ns
                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
                                                surf_lsm_v(l)%building_covered(m) ) 
                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
                                                surf_lsm_v(l)%building_covered(m) ) 
!
!--             If surface element is not a vegetation surface and any value in 
!--             vegetation_pars is given, neglect this information and give an
!--             informative message that this value will not be used.   
                IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m)  .AND.          &
                      ANY( vegetation_pars_f%pars_xy(:,j,i) /=                 &
                      vegetation_pars_f%fill ) )  THEN
                   WRITE( message_string, * )                                  &
                                 'surface element at grid point (j,i) = (',    &
                                 j, i, ') is not a vegation surface, ',        &
                                 'so that information given in ',              &
                                 'vegetation_pars at this point is neglected.' 
                   CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
                ELSE

                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=         &
                        vegetation_pars_f%fill )                               &
                      surf_lsm_v(l)%r_canopy_min(m)  =                         &
                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=         &
                        vegetation_pars_f%fill )                               &
                      surf_lsm_v(l)%lai(m)           =                         &
                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
                   IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=          &
                        vegetation_pars_f%fill )                               &
                      surf_lsm_v(l)%c_veg(m)         =                         &
                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
                   IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=             &
                        vegetation_pars_f%fill )                               &
                     surf_lsm_v(l)%g_d(m)            =                         &
                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
                   IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=             &
                        vegetation_pars_f%fill )                               &
                      surf_lsm_v(l)%z0(m)            =                         &
                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
                   IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=           &
                        vegetation_pars_f%fill )  THEN
                      surf_lsm_v(l)%z0h(m)           =                         &
                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
                      surf_lsm_v(l)%z0q(m)           =                         &
                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
                   ENDIF
                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=       &
                        vegetation_pars_f%fill )                               &
                      surf_lsm_v(l)%lambda_surface_s(m)  =                     &
                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=       &
                        vegetation_pars_f%fill )                               &
                      surf_lsm_v(l)%lambda_surface_u(m)  =                     &
                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
                   IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=        &
                        vegetation_pars_f%fill )                               &
                      surf_lsm_v(l)%f_sw_in(m)           =                     &
                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
                   IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=         &
                        vegetation_pars_f%fill )                               &
                      surf_lsm_v(l)%c_surface(m)         =                     &
                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
                   IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=             &
                        vegetation_pars_f%fill )                               &
                      surf_lsm_v(l)%albedo_type(ind_veg_wall,m) =              &
                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
                   IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=           &
                        vegetation_pars_f%fill )                               &
                      surf_lsm_v(l)%emissivity(ind_veg_wall,m)  =              &
                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
                ENDIF

             ENDDO
          ENDDO
       ENDIF 

!
!--    Level 1, initialization of water parameters. A horizontally 
!--    homogeneous distribution is assumed here. 
       IF ( water_type /= 0 )  THEN

          IF ( water_temperature == 9999999.9_wp )  THEN
             water_temperature = water_pars(ind_w_temp,water_type)       
          ENDIF

          IF ( z0_water == 9999999.9_wp )  THEN
             z0_water = water_pars(ind_w_z0,water_type)       
          ENDIF       

          IF ( z0h_water == 9999999.9_wp )  THEN
             z0h_water = water_pars(ind_w_z0h,water_type)       
          ENDIF  
          
          IF ( z0q_water == 9999999.9_wp )  THEN
             z0q_water = water_pars(ind_w_z0h,water_type)       
          ENDIF 

          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
             albedo_type = INT(water_pars(ind_w_at,water_type))        
          ENDIF
   
          IF ( emissivity == 9999999.9_wp )  THEN
             emissivity = water_pars(ind_w_emis,water_type)        
          ENDIF

       ENDIF 
!
!--    Map values onto horizontal elemements
       DO  m = 1, surf_lsm_h%ns
          IF ( surf_lsm_h%water_surface(m) )  THEN
             IF ( TRIM( initializing_actions ) /= 'read_restart_data' )        &
                t_soil_h%var_2d(:,m)        = water_temperature
             surf_lsm_h%z0(m)               = z0_water
             surf_lsm_h%z0h(m)              = z0h_water
             surf_lsm_h%z0q(m)              = z0q_water
             surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
             surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp               
             surf_lsm_h%c_surface(m)        = 0.0_wp
             surf_lsm_h%albedo_type(ind_wat_win,m) = albedo_type
             surf_lsm_h%emissivity(ind_wat_win,m)  = emissivity
             
             surf_lsm_h%water_type(m)      = water_type
             surf_lsm_h%water_type_name(m) = water_type_name(water_type)
          ENDIF
       ENDDO
!
!--    Map values onto vertical elements, even though this does not make
!--    much sense.
       DO  l = 0, 3
          DO  m = 1, surf_lsm_v(l)%ns
             IF ( surf_lsm_v(l)%water_surface(m) )  THEN
                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
                   t_soil_v(l)%var_2d(:,m)           = water_temperature
                surf_lsm_v(l)%z0(m)               = z0_water
                surf_lsm_v(l)%z0h(m)              = z0h_water
                surf_lsm_v(l)%z0q(m)              = z0q_water
                surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp
                surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp               
                surf_lsm_v(l)%c_surface(m)        = 0.0_wp
                surf_lsm_v(l)%albedo_type(ind_wat_win,m) = albedo_type
                surf_lsm_v(l)%emissivity(ind_wat_win,m)  = emissivity
                
                surf_lsm_v(l)%water_type(m)      = water_type
                surf_lsm_v(l)%water_type_name(m) = water_type_name(water_type)
             ENDIF 
          ENDDO
       ENDDO
!
!
!--    Level 2, initialization of water parameters via water_type read
!--    from file. Water surfaces are initialized for each (y,x)-grid point 
!--    individually using default paramter settings according to the given
!--    water type.
!--    Note, parameter 3/4 of water_pars are albedo and emissivity,
!--    whereas paramter 3/4 of water_pars_f are heat conductivities!
       IF ( water_type_f%from_file )  THEN
!
!--       Horizontal surfaces
          DO  m = 1, surf_lsm_h%ns
             i = surf_lsm_h%i(m)
             j = surf_lsm_h%j(m)
             
             st = water_type_f%var(j,i)
             IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
                   t_soil_h%var_2d(:,m) = water_pars(ind_w_temp,st)
                surf_lsm_h%z0(m)     = water_pars(ind_w_z0,st)
                surf_lsm_h%z0h(m)    = water_pars(ind_w_z0h,st)
                surf_lsm_h%z0q(m)    = water_pars(ind_w_z0h,st)
                surf_lsm_h%lambda_surface_s(m) = water_pars(ind_w_lambda_s,st)
                surf_lsm_h%lambda_surface_u(m) = water_pars(ind_w_lambda_u,st)              
                surf_lsm_h%c_surface(m)        = 0.0_wp
                surf_lsm_h%albedo_type(ind_wat_win,m) = INT( water_pars(ind_w_at,st) )
                surf_lsm_h%emissivity(ind_wat_win,m)  = water_pars(ind_w_emis,st)
                
                surf_lsm_h%water_type(m)      = st
                surf_lsm_h%water_type_name(m) = water_type_name(st)
             ENDIF
          ENDDO
!
!--       Vertical surfaces
          DO  l = 0, 3
             DO  m = 1, surf_lsm_v(l)%ns
                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
                                                surf_lsm_v(l)%building_covered(m) ) 
                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
                                                surf_lsm_v(l)%building_covered(m) ) 
             
                st = water_type_f%var(j,i)
                IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
                   IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  &
                      t_soil_v(l)%var_2d(:,m) = water_pars(ind_w_temp,st)
                   surf_lsm_v(l)%z0(m)     = water_pars(ind_w_z0,st)
                   surf_lsm_v(l)%z0h(m)    = water_pars(ind_w_z0h,st)
                   surf_lsm_v(l)%z0q(m)    = water_pars(ind_w_z0h,st)
                   surf_lsm_v(l)%lambda_surface_s(m) =                         &
                                                   water_pars(ind_w_lambda_s,st)
                   surf_lsm_v(l)%lambda_surface_u(m) =                         &
                                                   water_pars(ind_w_lambda_u,st)           
                   surf_lsm_v(l)%c_surface(m)     = 0.0_wp
                   surf_lsm_v(l)%albedo_type(ind_wat_win,m) =                  &
                                                  INT( water_pars(ind_w_at,st) )
                   surf_lsm_v(l)%emissivity(ind_wat_win,m)  =                  &
                                                  water_pars(ind_w_emis,st)
                                                  
                   surf_lsm_v(l)%water_type(m)      = st
                   surf_lsm_v(l)%water_type_name(m) = water_type_name(st)
                ENDIF
             ENDDO
          ENDDO
       ENDIF      

!
!--    Level 3, initialization of water parameters at single (x,y) 
!--    position via water_pars read from file.
       IF ( water_pars_f%from_file )  THEN
!
!--       Horizontal surfaces 
          DO  m = 1, surf_lsm_h%ns
             i = surf_lsm_h%i(m)
             j = surf_lsm_h%j(m)
!
!--          If surface element is not a water surface and any value in 
!--          water_pars is given, neglect this information and give an
!--          informative message that this value will not be used.   
             IF ( .NOT. surf_lsm_h%water_surface(m)  .AND.                     &
                   ANY( water_pars_f%pars_xy(:,j,i) /= water_pars_f%fill ) )  THEN
                WRITE( message_string, * )                                     &
                              'surface element at grid point (j,i) = (',       &
                              j, i, ') is not a water surface, ',              &
                              'so that information given in ',                 &
                              'water_pars at this point is neglected.' 
                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
             ELSE
                IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                   &
                     water_pars_f%fill  .AND.                                  &
                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
                      t_soil_h%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)

                IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /= water_pars_f%fill ) &
                   surf_lsm_h%z0(m)     = water_pars_f%pars_xy(ind_w_z0,j,i)

                IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /= water_pars_f%fill )&
                THEN
                   surf_lsm_h%z0h(