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

Last change on this file since 2705 was 2705, checked in by maronga, 4 years ago

minor bugfix for restarts

  • Property svn:keywords set to Id
File size: 293.9 KB
Line 
1!> @file land_surface_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: land_surface_model_mod.f90 2705 2017-12-18 11:26:23Z maronga $
27! Bugfix in binary output (wrong sequence)
28!
29! 2696 2017-12-14 17:12:51Z kanani
30! Bugfix: missing USE statement for calc_mean_profile
31! do not write surface temperatures onto pt array as this might cause
32! problems with nesting (MS)
33! Revised calculation of pt1 and qv1 (now done in surface_layer_fluxes). Bugfix
34! in calculation of surface density (cannot be done via an surface non-air
35! temperature) (BM)
36! Bugfix: g_d was NaN for non-vegetaed surface types (BM)
37! Bugfix initialization of c_veg and lai
38! Revise data output to enable _FillValues
39! Bugfix in calcultion of r_a and rad_net_l (MS)
40! Bugfix: rad_net is not updated in case of radiation_interaction and must thu
41! be calculated again from the radiative fluxes
42! Temporary fix for cases where no soil model is used on some PEs (BM)
43! Revised input and initialization of soil and surface paramters
44! pavement_depth is variable for each surface element
45! radiation quantities belong to surface type now
46! surface fractions initialized
47! Rename lsm_last_actions into lsm_write_restart_data (MS)
48!
49! 2608 2017-11-13 14:04:26Z schwenkel
50! Calculation of magnus equation in external module (diagnostic_quantities_mod).
51! Adjust calculation of vapor pressure and saturation mixing ratio that it is
52! consistent with formulations in other parts of PALM.
53!
54! 2575 2017-10-24 09:57:58Z maronga
55! Pavement parameterization revised
56!
57! 2573 2017-10-20 15:57:49Z scharf
58! bugfixes in last_actions
59!
60! 2548 2017-10-16 13:18:20Z suehring
61! extended by cloud_droplets option
62!
63! 2532 2017-10-11 16:00:46Z scharf
64! bugfixes in data_output_3d
65!
66! 2516 2017-10-04 11:03:04Z suehring
67! Remove tabs
68!
69! 2514 2017-10-04 09:52:37Z suehring
70! upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny
71! no output of ghost layer data
72!
73! 2504 2017-09-27 10:36:13Z maronga
74! Support roots and water under pavement. Added several pavement types.
75!
76! 2476 2017-09-18 07:54:32Z maronga
77! Bugfix for last commit
78!
79! 2475 2017-09-18 07:42:36Z maronga
80! Bugfix: setting of vegetation_pars for bare soil corrected.
81!
82! 2354 2017-08-17 10:49:36Z schwenkel
83! minor bugfixes
84!
85! 2340 2017-08-07 17:11:13Z maronga
86! Revised root_distribution tabel and implemented a pseudo-generic root fraction
87! calculation
88!
89! 2333 2017-08-04 09:08:26Z maronga
90! minor bugfixes
91!
92! 2332 2017-08-03 21:15:22Z maronga
93! bugfix in pavement_pars
94!
95! 2328 2017-08-03 12:34:22Z maronga
96! Revised skin layer concept.
97! Bugfix for runs with pavement surface and humidity
98! Revised some standard values in vegetation_pars
99! Added emissivity and default albedo_type as variable to tables
100! Changed default surface type to vegetation
101! Revised input of soil layer configuration
102!
103! 2307 2017-07-07 11:32:10Z suehring
104! Bugfix, variable names corrected
105!
106! 2299 2017-06-29 10:14:38Z maronga
107! Removed pt_p from USE statement. Adjusted call to lsm_soil_model to allow
108! spinups without soil moisture prediction
109!
110! 2298 2017-06-29 09:28:18Z raasch
111! type of write_binary changed from CHARACTER to LOGICAL
112!
113! 2296 2017-06-28 07:53:56Z maronga
114! Bugfix in calculation of bare soil heat capacity.
115! Bugfix in calculation of shf
116! Added support for spinups
117!
118! 2282 2017-06-13 11:38:46Z schwenkel
119! Bugfix for check of saturation moisture
120!
121! 2273 2017-06-09 12:46:06Z sward
122! Error number changed
123!
124! 2270 2017-06-09 12:18:47Z maronga
125! Revised parameterization of heat conductivity between skin layer and soil.
126! Temperature and moisture are now defined at the center of the layers.
127! Renamed veg_type to vegetation_type and pave_type to pavement_type_name
128! Renamed and reduced the number of look-up tables (vegetation_pars, soil_pars)
129! Revised land surface model initialization
130! Removed output of shf_eb and qsws_eb and removed _eb throughout code
131! Removed Clapp & Hornberger parameterization
132!
133! 2249 2017-06-06 13:58:01Z sward
134!
135! 2248 2017-06-06 13:52:54Z sward $
136! Error no changed
137!
138! 2246 2017-06-06 13:09:34Z sward
139! Error no changed
140!
141! Changed soil configuration to 8 layers. The number of soil layers is now
142! freely adjustable via the NAMELIST.
143!
144! 2237 2017-05-31 10:34:53Z suehring
145! Bugfix in write restart data
146!
147! 2233 2017-05-30 18:08:54Z suehring
148!
149! 2232 2017-05-30 17:47:52Z suehring
150! Adjustments to new topography and surface concept
151!   - now, also vertical walls are possible
152!   - for vertical walls, parametrization of r_a (aerodynamic resisistance) is
153!     implemented.
154!
155! Add check for soil moisture, it must not exceed its saturation value.
156!
157! 2149 2017-02-09 16:57:03Z scharf
158! Land surface parameters II corrected for vegetation_type 18 and 19
159!
160! 2031 2016-10-21 15:11:58Z knoop
161! renamed variable rho to rho_ocean
162!
163! 2000 2016-08-20 18:09:15Z knoop
164! Forced header and separation lines into 80 columns
165!
166! 1978 2016-07-29 12:08:31Z maronga
167! Bugfix: initial values of pave_surface and water_surface were not set.
168!
169! 1976 2016-07-27 13:28:04Z maronga
170! Parts of the code have been reformatted. Use of radiation model output is
171! generalized and simplified. Added more output quantities due to modularization
172!
173! 1972 2016-07-26 07:52:02Z maronga
174! Further modularization: output of cross sections and 3D data is now done in this
175! module. Moreover, restart data is written and read directly within this module.
176!
177!
178! 1966 2016-07-18 11:54:18Z maronga
179! Bugfix: calculation of m_total in soil model was not set to zero at model start
180!
181! 1949 2016-06-17 07:19:16Z maronga
182! Bugfix: calculation of qsws_soil_eb with precipitation = .TRUE. gave
183! qsws_soil_eb = 0 due to a typo
184!
185! 1856 2016-04-13 12:56:17Z maronga
186! Bugfix: for water surfaces, the initial water surface temperature is set equal
187! to the intital skin temperature. Moreover, the minimum value of r_a is now
188! 1.0 to avoid too large fluxes at the first model time step
189!
190! 1849 2016-04-08 11:33:18Z hoffmann
191! prr moved to arrays_3d
192!
193! 1826 2016-04-07 12:01:39Z maronga
194! Cleanup after modularization
195!
196! 1817 2016-04-06 15:44:20Z maronga
197! Added interface for lsm_init_arrays. Added subroutines for check_parameters,
198! header, and parin. Renamed some subroutines.
199!
200! 1788 2016-03-10 11:01:04Z maronga
201! Bugfix: calculate lambda_surface based on temperature gradient between skin
202! layer and soil layer instead of Obukhov length
203! Changed: moved calculation of surface specific humidity to energy balance solver
204! New: water surfaces are available by using a fixed sea surface temperature.
205! The roughness lengths are calculated dynamically using the Charnock
206! parameterization. This involves the new roughness length for moisture z0q.
207! New: modified solution of the energy balance solver and soil model for
208! paved surfaces (i.e. asphalt concrete).
209! Syntax layout improved.
210! Changed: parameter dewfall removed.
211!
212! 1783 2016-03-06 18:36:17Z raasch
213! netcdf variables moved to netcdf module
214!
215! 1757 2016-02-22 15:49:32Z maronga
216! Bugfix: set tm_soil_m to zero after allocation. Added parameter
217! unscheduled_radiation_calls to control calls of the radiation model based on
218! the skin temperature change during one time step (preliminary version). Set
219! qsws_soil_eb to zero at model start (previously set to qsws_eb). Removed MAX
220! function as it cannot be vectorized.
221!
222! 1709 2015-11-04 14:47:01Z maronga
223! Renamed pt_1 and qv_1 to pt1 and qv1.
224! Bugfix: set initial values for t_surface_p in case of restart runs
225! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky'
226! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky'
227! Added todo action
228!
229! 1697 2015-10-28 17:14:10Z raasch
230! bugfix: misplaced cpp-directive
231!
232! 1695 2015-10-27 10:03:11Z maronga
233! Bugfix: REAL constants provided with KIND-attribute in call of
234! Replaced rif with ol
235!
236! 1691 2015-10-26 16:17:44Z maronga
237! Added skip_time_do_lsm to allow for spin-ups without LSM. Various bugfixes:
238! Soil temperatures are now defined at the edges of the layers, calculation of
239! shb_eb corrected, prognostic equation for skin temperature corrected. Surface
240! fluxes are now directly transfered to atmosphere
241!
242! 1682 2015-10-07 23:56:08Z knoop
243! Code annotations made doxygen readable
244!
245! 1590 2015-05-08 13:56:27Z maronga
246! Bugfix: definition of character strings requires same length for all elements
247!
248! 1585 2015-04-30 07:05:52Z maronga
249! Modifications for RRTMG. Changed tables to PARAMETER type.
250!
251! 1571 2015-03-12 16:12:49Z maronga
252! Removed upper-case variable names. Corrected distribution of precipitation to
253! the liquid water reservoir and the bare soil fractions.
254!
255! 1555 2015-03-04 17:44:27Z maronga
256! Added output of r_a and r_s
257!
258! 1553 2015-03-03 17:33:54Z maronga
259! Improved better treatment of roughness lengths. Added default soil temperature
260! profile
261!
262! 1551 2015-03-03 14:18:16Z maronga
263! Flux calculation is now done in prandtl_fluxes. Added support for data output.
264! Vertical indices have been replaced. Restart runs are now possible. Some
265! variables have beem renamed. Bugfix in the prognostic equation for the surface
266! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
267! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
268! the hydraulic conductivity. Bugfix for root fraction and extraction
269! calculation
270!
271! intrinsic function MAX and MIN
272!
273! 1500 2014-12-03 17:42:41Z maronga
274! Corrected calculation of aerodynamic resistance (r_a).
275! Precipitation is now added to liquid water reservoir using LE_liq.
276! Added support for dry runs.
277!
278! 1496 2014-12-02 17:25:50Z maronga
279! Initial revision
280!
281!
282! Description:
283! ------------
284!> Land surface model, consisting of a solver for the energy balance at the
285!> surface and a multi layer soil scheme. The scheme is similar to the TESSEL
286!> scheme implemented in the ECMWF IFS model, with modifications according to
287!> H-TESSEL. The implementation is based on the formulation implemented in the
288!> DALES and UCLA-LES models.
289!>
290!> @todo Extensive verification energy-balance solver for vertical surfaces,
291!>       e.g. parametrization of r_a
292!> @todo Revise single land-surface processes for vertical surfaces, e.g.
293!>       treatment of humidity, etc.
294!> @todo Consider partial absorption of the net shortwave radiation by the
295!>       skin layer.
296!> @todo Improve surface water parameterization
297!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
298!>       nzt_soil=3)).
299!> @todo Implement surface runoff model (required when performing long-term LES
300!>       with considerable precipitation.
301!> @todo Revise calculation of f2 when wilting point is non-constant in the
302!>       soil
303!> @todo Allow for zero soil moisture (currently, it is set to wilting point)
304!> @note No time step criterion is required as long as the soil layers do not
305!>       become too thin.
306!> @todo Attention, pavement_subpars_1/2 are hardcoded to 8 levels, in case
307!>       more levels are used this may cause an potential bug
308!> @todo Routine calc_q_surface required?
309!------------------------------------------------------------------------------!
310 MODULE land_surface_model_mod
311 
312    USE arrays_3d,                                                             &
313        ONLY:  hyp, pt, prr, q, q_p, ql, vpt, u, v, w
314
315    USE calc_mean_profile_mod,                                                 &
316        ONLY:  calc_mean_profile
317
318    USE cloud_parameters,                                                      &
319        ONLY:  cp, hyrho, l_d_cp, l_d_r, l_v, pt_d_t, rho_l, r_d, r_v
320
321    USE control_parameters,                                                    &
322        ONLY:  cloud_droplets, cloud_physics, coupling_start_time, dt_3d,      &
323               end_time, humidity, intermediate_timestep_count,                &
324               initializing_actions, intermediate_timestep_count_max,          &
325               land_surface, max_masks, precipitation, pt_surface,             &
326               rho_surface, spinup, spinup_pt_mean, spinup_time,               &
327               surface_pressure, timestep_scheme, tsc,                         &
328               time_since_reference_point
329
330    USE indices,                                                               &
331        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb
332
333    USE netcdf_data_input_mod,                                                 &
334        ONLY :  building_type_f, init_3d, input_pids_static,                   &
335                netcdf_data_input_interpolate,                                 &
336                pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,  &
337                root_area_density_lsm_f, soil_pars_f, soil_type_f,             &
338                surface_fraction_f, vegetation_pars_f, vegetation_type_f,      &
339                water_pars_f, water_type_f
340
341    USE kinds
342
343    USE pegrid
344
345    USE radiation_model_mod,                                                   &
346        ONLY:  albedo, albedo_type, emissivity, force_radiation_call,          &
347               radiation_scheme, unscheduled_radiation_calls
348       
349    USE statistics,                                                            &
350        ONLY:  hom, statistic_regions
351
352    USE surface_mod,                                                           &
353        ONLY :  surf_lsm_h, surf_lsm_v, surf_type, surface_restore_elements
354
355    IMPLICIT NONE
356
357    TYPE surf_type_lsm
358       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  var_1d !< 1D prognostic variable
359       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2d !< 2D prognostic variable
360    END TYPE surf_type_lsm
361
362!
363!-- LSM model constants
364
365    REAL(wp), PARAMETER  ::                    &
366              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
367              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil (W/m/K) 
368              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix (W/m/K)
369              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water (W/m/K)
370              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
371              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil (J/m3/K)
372              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water (J/m3/K)
373              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
374
375
376    REAL(wp), DIMENSION(0:7), PARAMETER  :: dz_soil_default =                  & ! default soil layer configuration
377                                            (/ 0.01_wp, 0.02_wp, 0.04_wp,      &
378                                               0.06_wp, 0.14_wp, 0.26_wp,      &
379                                               0.54_wp, 1.86_wp/)
380
381    REAL(wp), DIMENSION(0:3), PARAMETER  :: dz_soil_ref =                      & ! reference four layer soil configuration used for estimating the root fractions
382                                            (/ 0.07_wp, 0.21_wp, 0.72_wp,      &
383                                               1.89_wp /)
384
385    REAL(wp), DIMENSION(0:3), PARAMETER  :: zs_ref =                           & ! reference four layer soil configuration used for estimating the root fractions
386                                            (/ 0.07_wp, 0.28_wp, 1.0_wp,       &
387                                               2.89_wp /)
388
389
390!
391!-- LSM variables
392    CHARACTER(10) :: surface_type = 'netcdf'      !< general classification. Allowed are:
393                                                  !< 'vegetation', 'pavement', ('building'),
394                                                  !< 'water', and 'netcdf'
395
396
397
398    INTEGER(iwp) :: nzb_soil = 0,             & !< bottom of the soil model (Earth's surface)
399                    nzt_soil = 7,             & !< top of the soil model
400                    nzt_pavement = 0,         & !< top of the pavement within the soil
401                    nzs = 8,                  & !< number of soil layers
402                    pavement_depth_level = 0, & !< default NAMELIST nzt_pavement
403                    pavement_type = 1,        & !< default NAMELIST pavement_type                 
404                    soil_type = 3,            & !< default NAMELIST soil_type
405                    vegetation_type = 2,      & !< default NAMELIST vegetation_type
406                    water_type = 1              !< default NAMELISt water_type
407                   
408   
409       
410    LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
411               constant_roughness = .FALSE.,     & !< use fixed/dynamic roughness lengths for water surfaces
412               force_radiation_call_l = .FALSE., & !< flag to force calling of radiation routine
413               aero_resist_kray = .TRUE.           !< flag to control parametrization of aerodynamic resistance at vertical surface elements
414
415!   value 9999999.9_wp -> generic available or user-defined value must be set
416!   otherwise -> no generic variable and user setting is optional
417    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
418                canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
419                c_surface = 9999999.9_wp,               & !< Surface (skin) heat capacity (J/m2/K)
420                drho_l_lv,                              & !< (rho_l * l_v)**-1
421                exn,                                    & !< value of the Exner function
422                e_s = 0.0_wp,                           & !< saturation water vapour pressure
423                field_capacity = 9999999.9_wp,          & !< NAMELIST m_fc
424                f_shortwave_incoming = 9999999.9_wp,    & !< NAMELIST f_sw_in
425                hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
426                ke = 0.0_wp,                            & !< Kersten number
427                lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil (W/m/K)
428                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s (W/m2/K)
429                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u (W/m2/K)
430                leaf_area_index = 9999999.9_wp,         & !< NAMELIST lai
431                l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
432                min_canopy_resistance = 9999999.9_wp,   & !< NAMELIST r_canopy_min
433                min_soil_resistance = 50.0_wp,          & !< NAMELIST r_soil_min
434                m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
435                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
436                pavement_heat_capacity = 9999999.9_wp,  & !< volumetric heat capacity of pavement (e.g. roads) (J/m3/K)
437                pavement_heat_conduct  = 9999999.9_wp,  & !< heat conductivity for pavements (e.g. roads) (W/m/K)
438                q_s = 0.0_wp,                           & !< saturation specific humidity
439                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
440                rho_cp,                                 & !< rho_surface * cp
441                rho_lv,                                 & !< rho_ocean * l_v
442                rd_d_rv,                                & !< r_d / r_v
443                saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
444                skip_time_do_lsm = 0.0_wp,              & !< LSM is not called before this time
445                vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
446                water_temperature = 9999999.9_wp,       & !< water temperature
447                wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
448                z0_vegetation  = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
449                z0h_vegetation = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
450                z0q_vegetation = 9999999.9_wp,          & !< NAMELIST z0q (lsm_par)
451                z0_pavement    = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
452                z0h_pavement   = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
453                z0q_pavement   = 9999999.9_wp,          & !< NAMELIST z0q (lsm_par)
454                z0_water       = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
455                z0h_water      = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
456                z0q_water      = 9999999.9_wp             !< NAMELIST z0q (lsm_par) 
457               
458               
459    REAL(wp), DIMENSION(:), ALLOCATABLE  :: ddz_soil_center, & !< 1/dz_soil_center
460                                            ddz_soil,        & !< 1/dz_soil
461                                            dz_soil_center,  & !< soil grid spacing (center-center)
462                                            zs,              & !< depth of the temperature/moisute levels
463                                            root_extr          !< root extraction
464
465
466                                           
467    REAL(wp), DIMENSION(0:20)  ::  root_fraction = 9999999.9_wp, & !< (NAMELIST) distribution of root surface area to the individual soil layers
468                                   soil_moisture = 0.0_wp,       & !< NAMELIST soil moisture content (m3/m3)
469                                   soil_temperature = 300.0_wp,  & !< NAMELIST soil temperature (K) +1
470                                   dz_soil  = 9999999.9_wp,      & !< (NAMELIST) soil layer depths (spacing)
471                                   zs_layer = 9999999.9_wp         !< soil layer depths (edge)
472                                 
473#if defined( __nopointer )
474    TYPE(surf_type_lsm), TARGET  ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
475                                     t_soil_h_p,  & !< Prog. soil temperature (K), horizontal surface elements
476                                     m_soil_h,    & !< Soil moisture (m3/m3), horizontal surface elements
477                                     m_soil_h_p     !< Prog. soil moisture (m3/m3), horizontal surface elements
478
479    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET  ::  &
480                                     t_soil_v,       & !< Soil temperature (K), vertical surface elements
481                                     t_soil_v_p,     & !< Prog. soil temperature (K), vertical surface elements
482                                     m_soil_v,       & !< Soil moisture (m3/m3), vertical surface elements
483                                     m_soil_v_p        !< Prog. soil moisture (m3/m3), vertical surface elements
484
485#else
486    TYPE(surf_type_lsm), POINTER ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
487                                     t_soil_h_p,  & !< Prog. soil temperature (K), horizontal surface elements
488                                     m_soil_h,    & !< Soil moisture (m3/m3), horizontal surface elements
489                                     m_soil_h_p     !< Prog. soil moisture (m3/m3), horizontal surface elements
490
491    TYPE(surf_type_lsm), TARGET  ::  t_soil_h_1,  & !<
492                                     t_soil_h_2,  & !<
493                                     m_soil_h_1,  & !<
494                                     m_soil_h_2     !<
495
496    TYPE(surf_type_lsm), DIMENSION(:), POINTER :: &
497                                     t_soil_v,    & !< Soil temperature (K), vertical surface elements
498                                     t_soil_v_p,  & !< Prog. soil temperature (K), vertical surface elements
499                                     m_soil_v,    & !< Soil moisture (m3/m3), vertical surface elements
500                                     m_soil_v_p     !< Prog. soil moisture (m3/m3), vertical surface elements   
501
502    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::&
503                                     t_soil_v_1,  & !<
504                                     t_soil_v_2,  & !<
505                                     m_soil_v_1,  & !<
506                                     m_soil_v_2     !<
507#endif   
508
509#if defined( __nopointer )
510    TYPE(surf_type_lsm), TARGET   ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements
511                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements
512                                      m_liq_h,        & !< liquid water reservoir (m), horizontal surface elements
513                                      m_liq_h_p         !< progn. liquid water reservoir (m), horizontal surface elements
514
515    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
516                                      t_surface_v,    & !< surface temperature (K), vertical surface elements
517                                      t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements
518                                      m_liq_v,        & !< liquid water reservoir (m), vertical surface elements
519                                      m_liq_v_p         !< progn. liquid water reservoir (m), vertical surface elements
520#else
521    TYPE(surf_type_lsm), POINTER  ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements
522                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements
523                                      m_liq_h,        & !< liquid water reservoir (m), horizontal surface elements
524                                      m_liq_h_p         !< progn. liquid water reservoir (m), horizontal surface elements
525
526    TYPE(surf_type_lsm), TARGET   ::  t_surface_h_1,  & !<
527                                      t_surface_h_2,  & !<
528                                      m_liq_h_1,      & !<
529                                      m_liq_h_2         !<
530
531    TYPE(surf_type_lsm), DIMENSION(:), POINTER  ::    &
532                                      t_surface_v,    & !< surface temperature (K), vertical surface elements
533                                      t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements
534                                      m_liq_v,        & !< liquid water reservoir (m), vertical surface elements
535                                      m_liq_v_p         !< progn. liquid water reservoir (m), vertical surface elements
536
537    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
538                                      t_surface_v_1,  & !<
539                                      t_surface_v_2,  & !<
540                                      m_liq_v_1,      & !<
541                                      m_liq_v_2         !<
542#endif
543
544#if defined( __nopointer )
545    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av
546#else
547    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av
548#endif
549
550#if defined( __nopointer )
551    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  t_soil_av, & !< Average of t_soil
552                                                        m_soil_av    !< Average of m_soil
553#else
554    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  t_soil_av, & !< Average of t_soil
555                                                        m_soil_av    !< Average of m_soil
556#endif
557
558    TYPE(surf_type_lsm), TARGET ::  tm_liq_h_m      !< liquid water reservoir tendency (m), horizontal surface elements
559    TYPE(surf_type_lsm), TARGET ::  tt_surface_h_m  !< surface temperature tendency (K), horizontal surface elements
560    TYPE(surf_type_lsm), TARGET ::  tt_soil_h_m     !< t_soil storage array, horizontal surface elements
561    TYPE(surf_type_lsm), TARGET ::  tm_soil_h_m     !< m_soil storage array, horizontal surface elements
562
563    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_liq_v_m      !< liquid water reservoir tendency (m), vertical surface elements
564    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_surface_v_m  !< surface temperature tendency (K), vertical surface elements
565    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_soil_v_m     !< t_soil storage array, vertical surface elements
566    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_soil_v_m     !< m_soil storage array, vertical surface elements
567
568!
569!-- Energy balance variables               
570    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
571              c_liq_av,         & !< average of c_liq
572              c_soil_av,        & !< average of c_soil
573              c_veg_av,         & !< average of c_veg
574              ghf_av,           & !< average of ghf
575              lai_av,           & !< average of lai
576              qsws_liq_av,      & !< average of qsws_liq
577              qsws_soil_av,     & !< average of qsws_soil
578              qsws_veg_av,      & !< average of qsws_veg
579              r_a_av,           & !< average of r_a
580              r_s_av              !< average of r_s
581                   
582
583!
584!-- Predefined Land surface classes (vegetation_type)
585    CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ &
586                                   'user defined              ',           & !  0
587                                   'bare soil                 ',           & !  1                           
588                                   'crops, mixed farming      ',           & !  2
589                                   'short grass               ',           & !  3
590                                   'evergreen needleleaf trees',           & !  4
591                                   'deciduous needleleaf trees',           & !  5
592                                   'evergreen broadleaf trees ',           & !  6
593                                   'deciduous broadleaf trees ',           & !  7
594                                   'tall grass                ',           & !  8
595                                   'desert                    ',           & !  9
596                                   'tundra                    ',           & ! 10
597                                   'irrigated crops           ',           & ! 11
598                                   'semidesert                ',           & ! 12
599                                   'ice caps and glaciers     ',           & ! 13
600                                   'bogs and marshes          ',           & ! 14
601                                   'evergreen shrubs          ',           & ! 15
602                                   'deciduous shrubs          ',           & ! 16
603                                   'mixed forest/woodland     ',           & ! 17
604                                   'interrupted forest        '            & ! 18
605                                                                 /)
606
607!
608!-- Soil model classes (soil_type)
609    CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ &
610                                   'user defined',                  & ! 0
611                                   'coarse      ',                  & ! 1
612                                   'medium      ',                  & ! 2
613                                   'medium-fine ',                  & ! 3
614                                   'fine        ',                  & ! 4
615                                   'very fine   ',                  & ! 5
616                                   'organic     '                   & ! 6
617                                                                 /)
618
619!
620!-- Pavement classes
621    CHARACTER(29), DIMENSION(0:15), PARAMETER :: pavement_type_name = (/ &
622                                   'user defined                 ', & ! 0
623                                   'asphalt/concrete mix         ', & ! 1
624                                   'asphalt (asphalt concrete)   ', & ! 2
625                                   'concrete (Portland concrete) ', & ! 3
626                                   'sett                         ', & ! 4
627                                   'paving stones                ', & ! 5
628                                   'cobblestone                  ', & ! 6
629                                   'metal                        ', & ! 7
630                                   'wood                         ', & ! 8
631                                   'gravel                       ', & ! 9
632                                   'fine gravel                  ', & ! 10
633                                   'pebblestone                  ', & ! 11
634                                   'woodchips                    ', & ! 12
635                                   'tartan (sports)              ', & ! 13
636                                   'artifical turf (sports)      ', & ! 14
637                                   'clay (sports)                '  & ! 15
638                                                                 /)                                                             
639                                                                 
640!
641!-- Water classes
642    CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ &
643                                   'user defined',                   & ! 0
644                                   'lake        ',                   & ! 1
645                                   'river       ',                   & ! 2
646                                   'ocean       ',                   & ! 3
647                                   'pond        ',                   & ! 4
648                                   'fountain    '                    & ! 5
649                                                                  /)                                                                 
650!
651!-- IDs for vegetation, pavement and water surfaces
652    INTEGER(iwp) ::  ind_veg = 0    !< index for vegetation surfaces, used to assess surface-fraction, albedo, etc.     
653    INTEGER(iwp) ::  ind_pav = 1    !< index for pavement surfaces, used to assess surface-fraction, albedo, etc.       
654    INTEGER(iwp) ::  ind_wat = 2    !< index for vegetation surfaces, used to assess surface-fraction, albedo, etc.                         
655                   
656!
657!-- Land surface parameters according to the respective classes (vegetation_type)
658    INTEGER(iwp) ::  ind_v_rc_min = 0    !< index for r_canopy_min in vegetation_pars
659    INTEGER(iwp) ::  ind_v_rc_lai = 1    !< index for LAI in vegetation_pars
660    INTEGER(iwp) ::  ind_v_c_veg   = 2   !< index for c_veg in vegetation_pars
661    INTEGER(iwp) ::  ind_v_gd  = 3       !< index for g_d in vegetation_pars
662    INTEGER(iwp) ::  ind_v_z0 = 4        !< index for z0 in vegetation_pars
663    INTEGER(iwp) ::  ind_v_z0qh = 5      !< index for z0h / z0q in vegetation_pars
664    INTEGER(iwp) ::  ind_v_lambda_s = 6  !< index for lambda_s_s in vegetation_pars
665    INTEGER(iwp) ::  ind_v_lambda_u = 7  !< index for lambda_s_u in vegetation_pars
666    INTEGER(iwp) ::  ind_v_f_sw_in = 8   !< index for f_sw_in in vegetation_pars
667    INTEGER(iwp) ::  ind_v_c_surf = 9    !< index for c_surface in vegetation_pars
668    INTEGER(iwp) ::  ind_v_at = 10       !< index for albedo_type in vegetation_pars
669    INTEGER(iwp) ::  ind_v_emis = 11     !< index for emissivity in vegetation_pars
670
671    INTEGER(iwp) ::  ind_w_temp     = 0    !< index for temperature in water_pars
672    INTEGER(iwp) ::  ind_w_z0       = 1    !< index for z0 in water_pars
673    INTEGER(iwp) ::  ind_w_z0h      = 2    !< index for z0h in water_pars
674    INTEGER(iwp) ::  ind_w_lambda_s = 3    !< index for lambda_s_s in water_pars
675    INTEGER(iwp) ::  ind_w_lambda_u = 4    !< index for lambda_s_u in water_pars
676    INTEGER(iwp) ::  ind_w_at       = 5    !< index for albedo type in water_pars
677    INTEGER(iwp) ::  ind_w_emis     = 6    !< index for emissivity in water_pars
678
679    INTEGER(iwp) ::  ind_p_z0       = 0    !< index for z0 in pavement_pars
680    INTEGER(iwp) ::  ind_p_z0h      = 1    !< index for z0h in pavement_pars
681    INTEGER(iwp) ::  ind_p_at       = 2    !< index for albedo type in pavement_pars
682    INTEGER(iwp) ::  ind_p_emis     = 3    !< index for emissivity in pavement_pars
683    INTEGER(iwp) ::  ind_p_lambda_h = 0    !< index for lambda_h in pavement_subsurface_pars
684    INTEGER(iwp) ::  ind_p_rho_c    = 1    !< index for rho_c in pavement_pars
685!
686!-- Land surface parameters
687!-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface, albedo_type, emissivity
688    REAL(wp), DIMENSION(0:11,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
689          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
690        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
691        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
692        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
693        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
694        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
695        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
696        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
697        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
698         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
699        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
700        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
701          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
702        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
703        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
704        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
705        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
706        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
707                                                               /), (/ 12, 18 /) )
708
709                                   
710!
711!-- Root distribution for default soil layer configuration (sum = 1)
712!--                                level 1 - level 4 according to zs_ref
713    REAL(wp), DIMENSION(0:3,1:18), PARAMETER :: root_distribution = RESHAPE( (/ &
714                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  1
715                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  2
716                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  3
717                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  4
718                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  5
719                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  6
720                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  7
721                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  8
722                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  9
723                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & ! 10
724                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 11
725                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 12
726                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 13
727                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 14
728                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 15
729                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16
730                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 17
731                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 18
732                                 /), (/ 4, 18 /) )
733
734!
735!-- Soil parameters according to the following porosity classes (soil_type)
736
737!
738!-- Soil parameters  alpha_vg,      l_vg,    n_vg, gamma_w_sat,    m_sat,     m_fc,   m_wilt,    m_res
739    REAL(wp), DIMENSION(0:7,1:6), PARAMETER :: soil_pars = RESHAPE( (/     &
740                      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
741                      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
742                      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
743                      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
744                      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
745                      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
746                                                                     /), (/ 8, 6 /) )
747
748
749!
750!-- TO BE FILLED
751!-- Pavement parameters      z0,       z0h, albedo_type, emissivity 
752    REAL(wp), DIMENSION(0:3,1:15), PARAMETER :: pavement_pars = RESHAPE( (/ &
753                      1.0E-4_wp, 1.0E-5_wp,     18.0_wp,    0.97_wp,  & !  1
754                      1.0E-4_wp, 1.0E-5_wp,     19.0_wp,    0.94_wp,  & !  2
755                      1.0E-4_wp, 1.0E-5_wp,     20.0_wp,    0.98_wp,  & !  3                                 
756                      1.0E-4_wp, 1.0E-5_wp,     21.0_wp,    0.93_wp,  & !  4
757                      1.0E-4_wp, 1.0E-5_wp,     22.0_wp,    0.97_wp,  & !  5
758                      1.0E-4_wp, 1.0E-5_wp,     23.0_wp,    0.97_wp,  & !  6
759                      1.0E-4_wp, 1.0E-5_wp,     24.0_wp,    0.97_wp,  & !  7
760                      1.0E-4_wp, 1.0E-5_wp,     25.0_wp,    0.94_wp,  & !  8
761                      1.0E-4_wp, 1.0E-5_wp,     26.0_wp,    0.98_wp,  & !  9                                 
762                      1.0E-4_wp, 1.0E-5_wp,     27.0_wp,    0.93_wp,  & ! 10
763                      1.0E-4_wp, 1.0E-5_wp,     28.0_wp,    0.97_wp,  & ! 11
764                      1.0E-4_wp, 1.0E-5_wp,     29.0_wp,    0.97_wp,  & ! 12
765                      1.0E-4_wp, 1.0E-5_wp,     30.0_wp,    0.97_wp,  & ! 13
766                      1.0E-4_wp, 1.0E-5_wp,     31.0_wp,    0.94_wp,  & ! 14
767                      1.0E-4_wp, 1.0E-5_wp,     32.0_wp,    0.98_wp   & ! 15
768                      /), (/ 4, 15 /) )                             
769!
770!-- Pavement subsurface parameters part 1: thermal conductivity (W/m/K)
771!--   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
772    REAL(wp), DIMENSION(0:7,1:15), PARAMETER :: pavement_subsurface_pars_1 = RESHAPE( (/ &
773       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & !  1
774       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & !  2
775       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & !  3
776       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
777       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
778       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
779       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
780       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
781       1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp,   1.00_wp, 9999999.9_wp, 9999999.9_wp, & !  9
782       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
783       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
784       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
785       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
786       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
787       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
788       /), (/ 8, 15 /) )
789
790!
791!-- Pavement subsurface parameters part 2: volumetric heat capacity (J/m3/K)
792!--     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
793    REAL(wp), DIMENSION(0:7,1:15), PARAMETER :: pavement_subsurface_pars_2 = RESHAPE( (/ &
794       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
795       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
796       1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 1.94E6_wp, 9999999.9_wp, 9999999.9_wp, & !  3
797       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
798       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
799       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
800       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
801       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
802       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
803       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
804       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
805       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
806       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
807       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
808       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
809                           /), (/ 8, 15 /) )
810 
811!
812!-- TO BE FILLED
813!-- Water parameters                    temperature,     z0,      z0h, albedo_type, emissivity,
814    REAL(wp), DIMENSION(0:6,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
815       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 1
816       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 2
817       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 3
818       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 4
819       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp  & ! 5
820                                                                     /), (/ 7, 5 /) )                                                                   
821                                                                                                                                     
822    SAVE
823
824
825    PRIVATE
826
827   
828!
829!-- Public functions
830    PUBLIC lsm_boundary_condition, lsm_check_data_output,                      &
831           lsm_check_data_output_pr,                                           &
832           lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,& 
833           lsm_data_output_2d, lsm_data_output_3d, lsm_energy_balance,         &
834           lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model,   &
835           lsm_swap_timelevel, lsm_read_restart_data, lsm_write_restart_data
836! !vegetat
837!-- Public parameters, constants and initial values
838    PUBLIC aero_resist_kray, skip_time_do_lsm
839
840!
841!-- Public grid variables
842    PUBLIC nzb_soil, nzs, nzt_soil, zs
843
844!
845!-- Public prognostic variables
846    PUBLIC m_soil_h, t_soil_h
847
848    INTERFACE lsm_boundary_condition
849       MODULE PROCEDURE lsm_boundary_condition
850    END INTERFACE lsm_boundary_condition
851
852    INTERFACE lsm_check_data_output
853       MODULE PROCEDURE lsm_check_data_output
854    END INTERFACE lsm_check_data_output
855   
856    INTERFACE lsm_check_data_output_pr
857       MODULE PROCEDURE lsm_check_data_output_pr
858    END INTERFACE lsm_check_data_output_pr
859   
860    INTERFACE lsm_check_parameters
861       MODULE PROCEDURE lsm_check_parameters
862    END INTERFACE lsm_check_parameters
863   
864    INTERFACE lsm_3d_data_averaging
865       MODULE PROCEDURE lsm_3d_data_averaging
866    END INTERFACE lsm_3d_data_averaging
867
868    INTERFACE lsm_data_output_2d
869       MODULE PROCEDURE lsm_data_output_2d
870    END INTERFACE lsm_data_output_2d
871
872    INTERFACE lsm_data_output_3d
873       MODULE PROCEDURE lsm_data_output_3d
874    END INTERFACE lsm_data_output_3d
875
876    INTERFACE lsm_define_netcdf_grid
877       MODULE PROCEDURE lsm_define_netcdf_grid
878    END INTERFACE lsm_define_netcdf_grid
879
880    INTERFACE lsm_energy_balance
881       MODULE PROCEDURE lsm_energy_balance
882    END INTERFACE lsm_energy_balance
883
884    INTERFACE lsm_header
885       MODULE PROCEDURE lsm_header
886    END INTERFACE lsm_header
887   
888    INTERFACE lsm_init
889       MODULE PROCEDURE lsm_init
890    END INTERFACE lsm_init
891
892    INTERFACE lsm_init_arrays
893       MODULE PROCEDURE lsm_init_arrays
894    END INTERFACE lsm_init_arrays
895   
896    INTERFACE lsm_parin
897       MODULE PROCEDURE lsm_parin
898    END INTERFACE lsm_parin
899   
900    INTERFACE lsm_soil_model
901       MODULE PROCEDURE lsm_soil_model
902    END INTERFACE lsm_soil_model
903
904    INTERFACE lsm_swap_timelevel
905       MODULE PROCEDURE lsm_swap_timelevel
906    END INTERFACE lsm_swap_timelevel
907
908    INTERFACE lsm_read_restart_data
909       MODULE PROCEDURE lsm_read_restart_data
910    END INTERFACE lsm_read_restart_data
911
912    INTERFACE lsm_write_restart_data
913       MODULE PROCEDURE lsm_write_restart_data
914    END INTERFACE lsm_write_restart_data
915
916 CONTAINS
917
918
919!------------------------------------------------------------------------------!
920! Description:
921! ------------
922!> Set internal Neumann boundary condition at outer soil grid points
923!> for temperature and humidity.
924!------------------------------------------------------------------------------!
925 SUBROUTINE lsm_boundary_condition
926 
927    IMPLICIT NONE
928
929    INTEGER(iwp) :: i      !< grid index x-direction
930    INTEGER(iwp) :: ioff   !< offset index x-direction indicating location of soil grid point
931    INTEGER(iwp) :: j      !< grid index y-direction
932    INTEGER(iwp) :: joff   !< offset index x-direction indicating location of soil grid point
933    INTEGER(iwp) :: k      !< grid index z-direction
934    INTEGER(iwp) :: koff   !< offset index x-direction indicating location of soil grid point
935    INTEGER(iwp) :: l      !< running index surface-orientation
936    INTEGER(iwp) :: m      !< running index surface elements
937
938    koff = surf_lsm_h%koff
939    DO  m = 1, surf_lsm_h%ns
940       i = surf_lsm_h%i(m)
941       j = surf_lsm_h%j(m)
942       k = surf_lsm_h%k(m)
943       pt(k+koff,j,i) = pt(k,j,i)
944    ENDDO
945
946    DO  l = 0, 3
947       ioff = surf_lsm_v(l)%ioff
948       joff = surf_lsm_v(l)%joff
949       DO  m = 1, surf_lsm_v(l)%ns
950          i = surf_lsm_v(l)%i(m)
951          j = surf_lsm_v(l)%j(m)
952          k = surf_lsm_v(l)%k(m)
953          pt(k,j+joff,i+ioff) = pt(k,j,i)
954       ENDDO
955    ENDDO
956!
957!-- In case of humidity, set boundary conditions also for q and vpt.
958    IF ( humidity )  THEN
959       koff = surf_lsm_h%koff
960       DO  m = 1, surf_lsm_h%ns
961          i = surf_lsm_h%i(m)
962          j = surf_lsm_h%j(m)
963          k = surf_lsm_h%k(m)
964          q(k+koff,j,i)   = q(k,j,i)
965          vpt(k+koff,j,i) = vpt(k,j,i)
966       ENDDO
967
968       DO  l = 0, 3
969          ioff = surf_lsm_v(l)%ioff
970          joff = surf_lsm_v(l)%joff
971          DO  m = 1, surf_lsm_v(l)%ns
972             i = surf_lsm_v(l)%i(m)
973             j = surf_lsm_v(l)%j(m)
974             k = surf_lsm_v(l)%k(m)
975             q(k,j+joff,i+ioff)   = q(k,j,i)
976             vpt(k,j+joff,i+ioff) = vpt(k,j,i)
977          ENDDO
978       ENDDO
979    ENDIF
980
981 END SUBROUTINE lsm_boundary_condition
982
983!------------------------------------------------------------------------------!
984! Description:
985! ------------
986!> Check data output for land surface model
987!------------------------------------------------------------------------------!
988 SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
989 
990 
991    USE control_parameters,                                                    &
992        ONLY:  data_output, message_string
993
994    IMPLICIT NONE
995
996    CHARACTER (LEN=*) ::  unit  !<
997    CHARACTER (LEN=*) ::  var   !<
998
999    INTEGER(iwp) :: i
1000    INTEGER(iwp) :: ilen   
1001    INTEGER(iwp) :: k
1002
1003    SELECT CASE ( TRIM( var ) )
1004
1005       CASE ( 'm_soil' )
1006          IF (  .NOT.  land_surface )  THEN
1007             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1008                      'res land_surface = .TRUE.'
1009             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1010          ENDIF
1011          unit = 'm3/m3'
1012           
1013       CASE ( 't_soil' )
1014          IF (  .NOT.  land_surface )  THEN
1015             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1016                      'res land_surface = .TRUE.'
1017             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1018          ENDIF
1019          unit = 'K'   
1020             
1021       CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf*', 'm_liq*',         &
1022              'qsws_liq*', 'qsws_soil*', 'qsws_veg*',                          &
1023              'r_a*', 'r_s*' )
1024          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
1025             message_string = 'illegal value for data_output: "' //            &
1026                              TRIM( var ) // '" & only 2d-horizontal ' //      &
1027                              'cross sections are allowed for this value'
1028             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
1029          ENDIF
1030          IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
1031             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1032                              'res land_surface = .TRUE.'
1033             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1034          ENDIF
1035          IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
1036             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1037                              'res land_surface = .TRUE.'
1038             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1039          ENDIF
1040          IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
1041             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1042                              'res land_surface = .TRUE.'
1043             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1044          ENDIF
1045          IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
1046             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1047                              'res land_surface = .TRUE.'
1048             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1049          ENDIF
1050          IF ( TRIM( var ) == 'ghf*'  .AND.  .NOT.  land_surface )  THEN
1051             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1052                              'res land_surface = .TRUE.'
1053             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1054          ENDIF
1055          IF ( TRIM( var ) == 'm_liq*'  .AND.  .NOT.  land_surface )  THEN
1056             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1057                              'res land_surface = .TRUE.'
1058             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1059          ENDIF
1060          IF ( TRIM( var ) == 'qsws_liq*'  .AND.  .NOT. land_surface )         &
1061          THEN
1062             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1063                              'res land_surface = .TRUE.'
1064             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1065          ENDIF
1066          IF ( TRIM( var ) == 'qsws_soil*'  .AND.  .NOT.  land_surface )       &
1067          THEN
1068             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1069                              'res land_surface = .TRUE.'
1070             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1071          ENDIF
1072          IF ( TRIM( var ) == 'qsws_veg*'  .AND.  .NOT. land_surface )         &
1073          THEN
1074             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1075                              'res land_surface = .TRUE.'
1076             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1077          ENDIF
1078          IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface )             &
1079          THEN
1080             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1081                              'res land_surface = .TRUE.'
1082             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1083          ENDIF
1084          IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface )             &
1085          THEN
1086             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
1087                              'res land_surface = .TRUE.'
1088             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
1089          ENDIF
1090
1091          IF ( TRIM( var ) == 'lai*'   )      unit = 'none' 
1092          IF ( TRIM( var ) == 'c_liq*' )      unit = 'none'
1093          IF ( TRIM( var ) == 'c_soil*')      unit = 'none'
1094          IF ( TRIM( var ) == 'c_veg*' )      unit = 'none'
1095          IF ( TRIM( var ) == 'ghf*')         unit = 'W/m2'
1096          IF ( TRIM( var ) == 'm_liq*'     )  unit = 'm'
1097          IF ( TRIM( var ) == 'qsws_liq*'  )  unit = 'W/m2'
1098          IF ( TRIM( var ) == 'qsws_soil*' )  unit = 'W/m2'
1099          IF ( TRIM( var ) == 'qsws_veg*'  )  unit = 'W/m2'
1100          IF ( TRIM( var ) == 'r_a*')         unit = 's/m'     
1101          IF ( TRIM( var ) == 'r_s*')         unit = 's/m' 
1102             
1103       CASE DEFAULT
1104          unit = 'illegal'
1105
1106    END SELECT
1107
1108
1109 END SUBROUTINE lsm_check_data_output
1110
1111
1112
1113!------------------------------------------------------------------------------!
1114! Description:
1115! ------------
1116!> Check data output of profiles for land surface model
1117!------------------------------------------------------------------------------!
1118 SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
1119 
1120    USE control_parameters,                                                    &
1121        ONLY:  data_output_pr, message_string
1122
1123    USE indices
1124
1125    USE profil_parameter
1126
1127    USE statistics
1128
1129    IMPLICIT NONE
1130   
1131    CHARACTER (LEN=*) ::  unit      !<
1132    CHARACTER (LEN=*) ::  variable  !<
1133    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
1134 
1135    INTEGER(iwp) ::  user_pr_index !<
1136    INTEGER(iwp) ::  var_count     !<
1137
1138    SELECT CASE ( TRIM( variable ) )
1139       
1140       CASE ( 't_soil', '#t_soil' )
1141          IF (  .NOT.  land_surface )  THEN
1142             message_string = 'data_output_pr = ' //                           &
1143                              TRIM( data_output_pr(var_count) ) // ' is' //    &
1144                              'not implemented for land_surface = .FALSE.'
1145             CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
1146          ELSE
1147             dopr_index(var_count) = 89
1148             dopr_unit     = 'K'
1149             hom(0:nzs-1,2,89,:)  = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
1150             IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
1151                dopr_initial_index(var_count) = 90
1152                hom(0:nzs-1,2,90,:)   = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
1153                data_output_pr(var_count)     = data_output_pr(var_count)(2:)
1154             ENDIF
1155             unit = dopr_unit
1156          ENDIF
1157
1158       CASE ( 'm_soil', '#m_soil' )
1159          IF (  .NOT.  land_surface )  THEN
1160             message_string = 'data_output_pr = ' //                           &
1161                              TRIM( data_output_pr(var_count) ) // ' is' //    &
1162                              ' not implemented for land_surface = .FALSE.'
1163             CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
1164          ELSE
1165             dopr_index(var_count) = 91
1166             dopr_unit     = 'm3/m3'
1167             hom(0:nzs-1,2,91,:)  = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
1168             IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
1169                dopr_initial_index(var_count) = 92
1170                hom(0:nzs-1,2,92,:)   = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
1171                data_output_pr(var_count)     = data_output_pr(var_count)(2:)
1172             ENDIF
1173             unit = dopr_unit
1174          ENDIF
1175
1176
1177       CASE DEFAULT
1178          unit = 'illegal'
1179
1180    END SELECT
1181
1182
1183 END SUBROUTINE lsm_check_data_output_pr
1184 
1185 
1186!------------------------------------------------------------------------------!
1187! Description:
1188! ------------
1189!> Check parameters routine for land surface model
1190!------------------------------------------------------------------------------!
1191 SUBROUTINE lsm_check_parameters
1192
1193    USE control_parameters,                                                    &
1194        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string,           &
1195               most_method
1196                     
1197   
1198    IMPLICIT NONE
1199
1200    INTEGER(iwp) ::  k        !< running index, z-dimension
1201
1202!
1203!-- Check for a valid setting of surface_type. The default value is 'netcdf'.
1204!-- In that case, the surface types are read from NetCDF file
1205    IF ( TRIM( surface_type ) /= 'vegetation'  .AND.                           &
1206         TRIM( surface_type ) /= 'pavement'    .AND.                           &
1207         TRIM( surface_type ) /= 'water'       .AND.                           &
1208         TRIM( surface_type ) /= 'netcdf' )  THEN 
1209       message_string = 'unknown surface type surface_type = "' //             &
1210                        TRIM( surface_type ) // '"'
1211       CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
1212    ENDIF
1213
1214!
1215!-- Dirichlet boundary conditions are required as the surface fluxes are
1216!-- calculated from the temperature/humidity gradients in the land surface
1217!-- model
1218    IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
1219       message_string = 'lsm requires setting of'//                            &
1220                        'bc_pt_b = "dirichlet" and '//                         &
1221                        'bc_q_b  = "dirichlet"'
1222       CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
1223    ENDIF
1224
1225    IF (  .NOT.  constant_flux_layer )  THEN
1226       message_string = 'lsm requires '//                                      &
1227                        'constant_flux_layer = .T.'
1228       CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
1229    ENDIF
1230
1231    IF ( TRIM( surface_type ) == 'vegetation' )  THEN
1232   
1233       IF ( vegetation_type == 0 )  THEN
1234          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
1235             message_string = 'vegetation_type = 0 (user defined)'//           &
1236                              'requires setting of min_canopy_resistance'//    &
1237                              '/= 9999999.9'
1238             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1239          ENDIF
1240
1241          IF ( leaf_area_index == 9999999.9_wp )  THEN
1242             message_string = 'vegetation_type = 0 (user_defined)'//           &
1243                              'requires setting of leaf_area_index'//          &
1244                              '/= 9999999.9'
1245             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1246          ENDIF
1247
1248          IF ( vegetation_coverage == 9999999.9_wp )  THEN
1249             message_string = 'vegetation_type = 0 (user_defined)'//           &
1250                              'requires setting of vegetation_coverage'//      &
1251                              '/= 9999999.9'
1252                CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1253          ENDIF
1254
1255          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
1256             message_string = 'vegetation_type = 0 (user_defined)'//           &
1257                              'requires setting of'//                          &
1258                              'canopy_resistance_coefficient /= 9999999.9'
1259             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1260          ENDIF
1261
1262          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
1263             message_string = 'vegetation_type = 0 (user_defined)'//           &
1264                              'requires setting of lambda_surface_stable'//    &
1265                              '/= 9999999.9'
1266             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1267          ENDIF
1268
1269          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
1270             message_string = 'vegetation_type = 0 (user_defined)'//           &
1271                              'requires setting of lambda_surface_unstable'//  &
1272                              '/= 9999999.9'
1273             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1274          ENDIF
1275
1276          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
1277             message_string = 'vegetation_type = 0 (user_defined)'//           &
1278                              'requires setting of f_shortwave_incoming'//     &
1279                              '/= 9999999.9'
1280             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1281          ENDIF
1282
1283          IF ( z0_vegetation == 9999999.9_wp )  THEN
1284             message_string = 'vegetation_type = 0 (user_defined)'//           &
1285                              'requires setting of z0_vegetation'//            &
1286                              '/= 9999999.9'
1287             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1288          ENDIF
1289
1290          IF ( z0h_vegetation == 9999999.9_wp )  THEN
1291             message_string = 'vegetation_type = 0 (user_defined)'//           &
1292                              'requires setting of z0h_vegetation'//           &
1293                              '/= 9999999.9'
1294             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1295          ENDIF
1296       ENDIF
1297
1298       IF ( vegetation_type == 1 )  THEN
1299          IF ( vegetation_coverage /= 9999999.9_wp  .AND.  vegetation_coverage &
1300               /= 0.0_wp )  THEN
1301             message_string = 'vegetation_type = 1 (bare soil)'//              &
1302                              ' requires vegetation_coverage = 0'
1303             CALL message( 'check_parameters', 'PA0471', 1, 2, 0, 6, 0 )
1304          ENDIF
1305       ENDIF
1306 
1307    ENDIF
1308   
1309    IF ( TRIM( surface_type ) == 'water' )  THEN
1310
1311       IF ( TRIM( most_method ) == 'lookup' )  THEN   
1312          WRITE( message_string, * ) 'surface_type = ', surface_type,          &
1313                                     ' is not allowed in combination with ',   &
1314                                     'most_method = ', most_method
1315          CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
1316       ENDIF
1317
1318       IF ( water_type == 0 )  THEN 
1319       
1320          IF ( z0_water == 9999999.9_wp )  THEN
1321             message_string = 'water_type = 0 (user_defined)'//                &
1322                              'requires setting of z0_water'//                 &
1323                              '/= 9999999.9'
1324             CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
1325          ENDIF
1326
1327          IF ( z0h_water == 9999999.9_wp )  THEN
1328             message_string = 'water_type = 0 (user_defined)'//                &
1329                              'requires setting of z0h_water'//                &
1330                              '/= 9999999.9'
1331             CALL message( 'check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
1332          ENDIF
1333         
1334          IF ( water_temperature == 9999999.9_wp )  THEN
1335             message_string = 'water_type = 0 (user_defined)'//                &
1336                              'requires setting of water_temperature'//        &
1337                              '/= 9999999.9'
1338             CALL message( 'check_parameters', 'PA0379', 1, 2, 0, 6, 0 )
1339          ENDIF       
1340         
1341       ENDIF
1342       
1343    ENDIF
1344   
1345    IF ( TRIM( surface_type ) == 'pavement' )  THEN
1346
1347       IF ( ANY( dz_soil /= 9999999.9_wp )  .AND.  pavement_type /= 0 )  THEN
1348          message_string = 'non-default setting of dz_soil '//                  &
1349                           'does not allow to use pavement_type /= 0)'
1350             CALL message( 'check_parameters', 'PA0341', 1, 2, 0, 6, 0 )
1351          ENDIF
1352
1353       IF ( pavement_type == 0 )  THEN 
1354       
1355          IF ( z0_pavement == 9999999.9_wp )  THEN
1356             message_string = 'pavement_type = 0 (user_defined)'//             &
1357                              'requires setting of z0_pavement'//              &
1358                              '/= 9999999.9'
1359             CALL message( 'check_parameters', 'PA0352', 1, 2, 0, 6, 0 )
1360          ENDIF
1361
1362          IF ( z0h_pavement == 9999999.9_wp )  THEN
1363             message_string = 'pavement_type = 0 (user_defined)'//             &
1364                              'requires setting of z0h_pavement'//             &
1365                              '/= 9999999.9'
1366             CALL message( 'check_parameters', 'PA0353', 1, 2, 0, 6, 0 )
1367          ENDIF
1368         
1369          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
1370             message_string = 'pavement_type = 0 (user_defined)'//             &
1371                              'requires setting of pavement_heat_conduct'//    &
1372                              '/= 9999999.9'
1373             CALL message( 'check_parameters', 'PA0342', 1, 2, 0, 6, 0 )
1374          ENDIF
1375         
1376           IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
1377             message_string = 'pavement_type = 0 (user_defined)'//             &
1378                              'requires setting of pavement_heat_capacity'//   &
1379                              '/= 9999999.9'
1380             CALL message( 'check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
1381          ENDIF
1382
1383          IF ( pavement_depth_level == 0 )  THEN
1384             message_string = 'pavement_type = 0 (user_defined)'//             &
1385                              'requires setting of pavement_depth_level'//     &
1386                              '/= 0'
1387             CALL message( 'check_parameters', 'PA0474', 1, 2, 0, 6, 0 )
1388          ENDIF
1389
1390       ENDIF
1391   
1392    ENDIF
1393
1394    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
1395       IF ( ANY( water_type_f%var /= water_type_f%fill )  .AND.                &
1396            TRIM( most_method ) == 'lookup' )  THEN   
1397          WRITE( message_string, * ) 'water-surfaces are not allowed in ' //   &
1398                                     'combination with most_method = ',        &
1399                                     TRIM( most_method )
1400          CALL message( 'check_parameters', 'PA0999', 2, 2, 0, 6, 0 )
1401       ENDIF
1402!
1403!--    MS: Some problme here, after calling message everythings stucks at
1404!--        MPI_FINALIZE call.
1405       IF ( ANY( pavement_type_f%var /= pavement_type_f%fill )  .AND.           &
1406            ANY( dz_soil /= 9999999.9_wp ) )  THEN
1407          message_string = 'pavement-surfaces are not allowed in ' //           &
1408                           'combination with a non-default setting of dz_soil'
1409          CALL message( 'check_parameters', 'PA0999', 2, 2, 0, 6, 0 )
1410       ENDIF
1411    ENDIF
1412   
1413!
1414!-- Temporary message as long as NetCDF input is not available
1415    IF ( TRIM( surface_type ) == 'netcdf'  .AND.  .NOT.  input_pids_static )   &
1416    THEN
1417       message_string = 'surface_type = netcdf requires static input file.'
1418       CALL message( 'check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
1419    ENDIF
1420
1421    IF ( soil_type == 0 )  THEN
1422
1423       IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
1424          message_string = 'soil_type = 0 (user_defined)'//                    &
1425                           'requires setting of alpha_vangenuchten'//          &
1426                           '/= 9999999.9'
1427          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1428       ENDIF
1429
1430       IF ( l_vangenuchten == 9999999.9_wp )  THEN
1431          message_string = 'soil_type = 0 (user_defined)'//                    &
1432                           'requires setting of l_vangenuchten'//              &
1433                           '/= 9999999.9'
1434          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1435       ENDIF
1436
1437       IF ( n_vangenuchten == 9999999.9_wp )  THEN
1438          message_string = 'soil_type = 0 (user_defined)'//                    &
1439                           'requires setting of n_vangenuchten'//              &
1440                           '/= 9999999.9'
1441          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1442       ENDIF
1443
1444       IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
1445          message_string = 'soil_type = 0 (user_defined)'//                    &
1446                           'requires setting of hydraulic_conductivity'//      &
1447                           '/= 9999999.9'
1448          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1449       ENDIF
1450
1451       IF ( saturation_moisture == 9999999.9_wp )  THEN
1452          message_string = 'soil_type = 0 (user_defined)'//                    &
1453                           'requires setting of saturation_moisture'//         &
1454                           '/= 9999999.9'
1455          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1456       ENDIF
1457
1458       IF ( field_capacity == 9999999.9_wp )  THEN
1459          message_string = 'soil_type = 0 (user_defined)'//                    &
1460                           'requires setting of field_capacity'//              &
1461                           '/= 9999999.9'
1462          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1463       ENDIF
1464
1465       IF ( wilting_point == 9999999.9_wp )  THEN
1466          message_string = 'soil_type = 0 (user_defined)'//                    &
1467                           'requires setting of wilting_point'//               &
1468                           '/= 9999999.9'
1469          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1470       ENDIF
1471
1472       IF ( residual_moisture == 9999999.9_wp )  THEN
1473          message_string = 'soil_type = 0 (user_defined)'//                    &
1474                           'requires setting of residual_moisture'//           &
1475                           '/= 9999999.9'
1476          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1477       ENDIF
1478
1479    ENDIF
1480
1481!
1482!-- Determine number of soil layers to be used and check whether an appropriate
1483!-- root fraction is prescribed
1484    nzb_soil = 0
1485    nzt_soil = -1
1486    IF ( ALL( dz_soil == 9999999.9_wp ) )  THEN
1487       nzt_soil = 7
1488       dz_soil(nzb_soil:nzt_soil) = dz_soil_default
1489    ELSE
1490       DO k = 0, 19
1491          IF ( dz_soil(k) /= 9999999.9_wp )  THEN
1492             nzt_soil = nzt_soil + 1
1493          ENDIF
1494       ENDDO   
1495    ENDIF
1496    nzs = nzt_soil + 1
1497
1498
1499    IF ( vegetation_type == 0 )  THEN
1500       IF ( SUM( root_fraction(nzb_soil:nzt_soil) ) /= 1.0_wp )  THEN
1501          message_string = 'vegetation_type = 0 (user_defined)'//              &
1502                           'requires setting of root_fraction'//               &
1503                           '/= 9999999.9 and SUM(root_fraction) = 1'
1504          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1505       ENDIF
1506    ENDIF   
1507   
1508   
1509!
1510!-- Check for proper setting of soil moisture, must not be larger than its
1511!-- saturation value.
1512    DO  k = nzb_soil, nzt_soil
1513       IF ( soil_moisture(k) > saturation_moisture )  THEN
1514          message_string = 'soil_moisture must not exceed its saturation' //    &
1515                            ' value'
1516          CALL message( 'check_parameters', 'PA0458', 1, 2, 0, 6, 0 )
1517       ENDIF
1518    ENDDO
1519 
1520!
1521!-- Calculate grid spacings. Temperature and moisture are defined at
1522!-- the center of the soil layers, whereas gradients/fluxes are
1523!-- defined at the edges (_layer)
1524!
1525!-- Allocate global 1D arrays
1526    ALLOCATE ( ddz_soil_center(nzb_soil:nzt_soil) )
1527    ALLOCATE ( ddz_soil(nzb_soil:nzt_soil+1) )
1528    ALLOCATE ( dz_soil_center(nzb_soil:nzt_soil) )
1529    ALLOCATE ( zs(nzb_soil:nzt_soil+1) )
1530
1531
1532    zs(nzb_soil) = 0.5_wp * dz_soil(nzb_soil)
1533    zs_layer(nzb_soil) = dz_soil(nzb_soil)
1534
1535    DO  k = nzb_soil+1, nzt_soil
1536       zs_layer(k) = zs_layer(k-1) + dz_soil(k)
1537       zs(k) = (zs_layer(k) +  zs_layer(k-1)) * 0.5_wp
1538    ENDDO
1539
1540    dz_soil(nzt_soil+1) = zs_layer(nzt_soil) + dz_soil(nzt_soil)
1541    zs(nzt_soil+1) = zs_layer(nzt_soil) + 0.5_wp * dz_soil(nzt_soil)
1542 
1543    DO  k = nzb_soil, nzt_soil-1
1544       dz_soil_center(k) = zs(k+1) - zs(k)
1545       IF ( dz_soil_center(k) == 0.0_wp )  THEN
1546          message_string = 'invalid soil layer configuration found ' //        &
1547                           '(dz_soil_center(k) = 0.0)'
1548          CALL message( 'lsm_read_restart_data', 'PA0140', 1, 2, 0, 6, 0 )
1549       ENDIF 
1550    ENDDO
1551 
1552    dz_soil_center(nzt_soil) = zs_layer(k-1) + dz_soil(k) - zs(nzt_soil)
1553       
1554    ddz_soil_center = 1.0_wp / dz_soil_center
1555    ddz_soil(nzb_soil:nzt_soil) = 1.0_wp / dz_soil(nzb_soil:nzt_soil)
1556
1557
1558
1559 END SUBROUTINE lsm_check_parameters
1560 
1561!------------------------------------------------------------------------------!
1562! Description:
1563! ------------
1564!> Solver for the energy balance at the surface.
1565!------------------------------------------------------------------------------!
1566 SUBROUTINE lsm_energy_balance( horizontal, l )
1567
1568    USE diagnostic_quantities_mod,                                             &
1569        ONLY:  magnus
1570
1571    USE pegrid
1572
1573    IMPLICIT NONE
1574
1575    INTEGER(iwp) ::  i         !< running index
1576    INTEGER(iwp) ::  i_off     !< offset to determine index of surface element, seen from atmospheric grid point, for x
1577    INTEGER(iwp) ::  j         !< running index
1578    INTEGER(iwp) ::  j_off     !< offset to determine index of surface element, seen from atmospheric grid point, for y
1579    INTEGER(iwp) ::  k         !< running index
1580    INTEGER(iwp) ::  k_off     !< offset to determine index of surface element, seen from atmospheric grid point, for z
1581    INTEGER(iwp) ::  ks        !< running index
1582    INTEGER(iwp) ::  l         !< surface-facing index
1583    INTEGER(iwp) ::  m         !< running index concerning wall elements
1584
1585    LOGICAL      ::  horizontal !< Flag indicating horizontal or vertical surfaces
1586
1587    REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
1588                f1,          & !< resistance correction term 1
1589                f2,          & !< resistance correction term 2
1590                f3,          & !< resistance correction term 3
1591                m_min,       & !< minimum soil moisture
1592                e,           & !< water vapour pressure
1593                e_s,         & !< water vapour saturation pressure
1594                e_s_dt,      & !< derivate of e_s with respect to T
1595                tend,        & !< tendency
1596                dq_s_dt,     & !< derivate of q_s with respect to T
1597                coef_1,      & !< coef. for prognostic equation
1598                coef_2,      & !< coef. for prognostic equation
1599                f_qsws,      & !< factor for qsws
1600                f_qsws_veg,  & !< factor for qsws_veg
1601                f_qsws_soil, & !< factor for qsws_soil
1602                f_qsws_liq,  & !< factor for qsws_liq
1603                f_shf,       & !< factor for shf
1604                lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
1605                lambda_surface, & !< Current value of lambda_surface (W/m2/K)
1606                m_liq_max      !< maxmimum value of the liq. water reservoir
1607
1608    TYPE(surf_type_lsm), POINTER ::  surf_t_surface
1609    TYPE(surf_type_lsm), POINTER ::  surf_t_surface_p
1610    TYPE(surf_type_lsm), POINTER ::  surf_tt_surface_m
1611    TYPE(surf_type_lsm), POINTER ::  surf_m_liq
1612    TYPE(surf_type_lsm), POINTER ::  surf_m_liq_p
1613    TYPE(surf_type_lsm), POINTER ::  surf_tm_liq_m
1614
1615    TYPE(surf_type_lsm), POINTER ::  surf_m_soil
1616    TYPE(surf_type_lsm), POINTER ::  surf_t_soil
1617
1618    TYPE(surf_type), POINTER  ::  surf  !< surface-date type variable
1619
1620    IF ( horizontal )  THEN
1621       surf              => surf_lsm_h
1622
1623       surf_t_surface    => t_surface_h
1624       surf_t_surface_p  => t_surface_h_p
1625       surf_tt_surface_m => tt_surface_h_m
1626       surf_m_liq        => m_liq_h
1627       surf_m_liq_p      => m_liq_h_p
1628       surf_tm_liq_m     => tm_liq_h_m
1629       surf_m_soil       => m_soil_h
1630       surf_t_soil       => t_soil_h
1631    ELSE
1632       surf              => surf_lsm_v(l)
1633
1634       surf_t_surface    => t_surface_v(l)
1635       surf_t_surface_p  => t_surface_v_p(l)
1636       surf_tt_surface_m => tt_surface_v_m(l)
1637       surf_m_liq        => m_liq_v(l)
1638       surf_m_liq_p      => m_liq_v_p(l)
1639       surf_tm_liq_m     => tm_liq_v_m(l)
1640       surf_m_soil       => m_soil_v(l)
1641       surf_t_soil       => t_soil_v(l)
1642    ENDIF
1643
1644!
1645!-- Index offset of surface element point with respect to adjoining
1646!-- atmospheric grid point
1647    k_off = surf%koff
1648    j_off = surf%joff
1649    i_off = surf%ioff
1650
1651!
1652!-- Calculate the exner function for the current time step
1653    exn = ( surface_pressure / 1000.0_wp )**0.286_wp
1654
1655    DO  m = 1, surf%ns
1656
1657       i   = surf%i(m)           
1658       j   = surf%j(m)
1659       k   = surf%k(m)
1660
1661!
1662!--    Define heat conductivity between surface and soil depending on surface
1663!--    type. For vegetation, a skin layer parameterization is used. The new
1664!--    parameterization uses a combination of two conductivities: a constant
1665!--    conductivity for the skin layer, and a conductivity according to the
1666!--    uppermost soil layer. For bare soil and pavements, no skin layer is
1667!--    applied. In these cases, the temperature is assumed to be constant
1668!--    between the surface and the first soil layer. The heat conductivity is
1669!--    then derived from the soil/pavement properties.
1670!--    For water surfaces, the conductivity is already set to 1E10.
1671!--    Moreover, the heat capacity is set. For bare soil the heat capacity is
1672!--    the capacity of the uppermost soil layer, for pavement it is that of
1673!--    the material involved.
1674
1675!
1676!--    for vegetation type surfaces, the thermal conductivity of the soil is
1677!--    needed
1678
1679       IF ( surf%vegetation_surface(m) )  THEN
1680
1681          lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(nzb_soil,m)) *      &
1682                         lambda_h_water ** surf_m_soil%var_2d(nzb_soil,m)
1683                         
1684          ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(nzb_soil,m) /   &
1685                                                     surf%m_sat(nzb_soil,m) ) )                   
1686                         
1687          lambda_soil = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry )   &
1688                           * ddz_soil(nzb_soil) * 2.0_wp
1689
1690!
1691!--       When bare soil is set without a thermal conductivity (no skin layer),
1692!--       a heat capacity is that of the soil layer, otherwise it is a
1693!--       combination of the conductivities from the skin and the soil layer
1694          IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN
1695            surf%c_surface(m) = (rho_c_soil * (1.0_wp - surf%m_sat(nzb_soil,m))&
1696                              + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) ) &
1697                              * dz_soil(nzb_soil) * 0.5_wp   
1698            lambda_surface = lambda_soil
1699
1700          ELSE IF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m))&
1701          THEN
1702             lambda_surface = surf%lambda_surface_s(m) * lambda_soil           &
1703                              / ( surf%lambda_surface_s(m) + lambda_soil )
1704          ELSE
1705
1706             lambda_surface = surf%lambda_surface_u(m) * lambda_soil           &
1707                              / ( surf%lambda_surface_u(m) + lambda_soil )
1708          ENDIF
1709       ELSE
1710          lambda_surface = surf%lambda_surface_s(m)
1711       ENDIF
1712
1713!
1714!--    Set heat capacity of the skin/surface. It is ususally zero when a skin
1715!--    layer is used, and non-zero otherwise.
1716       c_surface_tmp = surf%c_surface(m) 
1717
1718!
1719!--    First step: calculate aerodyamic resistance. As pt, us, ts
1720!--    are not available for the prognostic time step, data from the last
1721!--    time step is used here. Note that this formulation is the
1722!--    equivalent to the ECMWF formulation using drag coefficients
1723!        IF ( cloud_physics )  THEN
1724!           pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
1725!           qv1 = q(k,j,i) - ql(k,j,i)
1726!        ELSEIF ( cloud_droplets ) THEN
1727!           pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
1728!           qv1 = q(k,j,i)
1729!        ELSE
1730!           pt1 = pt(k,j,i)
1731!           IF ( humidity )  THEN
1732!              qv1 = q(k,j,i)
1733!           ELSE
1734!              qv1 = 0.0_wp
1735!           ENDIF
1736!        ENDIF
1737!
1738!--    Calculate aerodynamical resistance. For horizontal and vertical
1739!--    surfaces MOST is applied. Moreover, for vertical surfaces, resistance
1740!--    can be obtain via parameterization of Mason (2000) /
1741!--    Krayenhoff and Voogt (2006).
1742!--    To do: detailed investigation which approach gives more reliable results!
1743!--    Please note, in case of very small friction velocity, e.g. in little
1744!--    holes, the resistance can become negative. For this reason, limit r_a
1745!--    to positive values.
1746       IF ( horizontal  .OR.  .NOT. aero_resist_kray )  THEN
1747          surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) /            &
1748                             ( surf%ts(m) * surf%us(m) + 1.0E-20_wp ) )
1749       ELSE
1750          surf%r_a(m) = 1.0_wp / ( 11.8_wp + 4.2_wp *                          &
1751                        SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
1752                                   ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
1753                                   ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2,  &
1754                              0.01_wp ) )                                      &
1755                                 )
1756       ENDIF
1757!
1758!--    Make sure that the resistance does not drop to zero for neutral
1759!--    stratification.
1760       IF ( surf%r_a(m) < 1.0_wp )  surf%r_a(m) = 1.0_wp
1761!
1762!--    Second step: calculate canopy resistance r_canopy
1763!--    f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
1764 
1765!--    f1: correction for incoming shortwave radiation (stomata close at
1766!--    night)
1767       f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) /          &
1768                        (0.81_wp * (0.004_wp * surf%rad_sw_in(m)               &
1769                         + 1.0_wp)) )
1770!
1771!--    f2: correction for soil moisture availability to plants (the
1772!--    integrated soil moisture must thus be considered here)
1773!--    f2 = 0 for very dry soils
1774       m_total = 0.0_wp
1775       DO  ks = nzb_soil, nzt_soil
1776           m_total = m_total + surf%root_fr(ks,m)                              &
1777                     * MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(ks,m) )
1778       ENDDO
1779
1780!
1781!--    The calculation of f2 is based on only one wilting point value for all
1782!--    soil layers. The value at k=nzb_soil is used here as a proxy but might
1783!--    need refinement in the future.
1784       IF ( m_total > surf%m_wilt(nzb_soil,m)  .AND.                           &
1785            m_total < surf%m_fc(nzb_soil,m) )  THEN
1786          f2 = ( m_total - surf%m_wilt(nzb_soil,m) ) /                         &
1787               ( surf%m_fc(nzb_soil,m) - surf%m_wilt(nzb_soil,m) )
1788       ELSEIF ( m_total >= surf%m_fc(nzb_soil,m) )  THEN
1789          f2 = 1.0_wp
1790       ELSE
1791          f2 = 1.0E-20_wp
1792       ENDIF
1793
1794!
1795!--    Calculate water vapour pressure at saturation and convert to hPa
1796       e_s = 0.01_wp * magnus( surf_t_surface%var_1d(m) )
1797
1798!
1799!--    f3: correction for vapour pressure deficit
1800       IF ( surf%g_d(m) /= 0.0_wp )  THEN
1801!
1802!--       Calculate vapour pressure
1803          e  = surf%qv1(m) * surface_pressure / ( surf%qv1(m) + 0.622_wp )
1804          f3 = EXP ( - surf%g_d(m) * (e_s - e) )
1805       ELSE
1806          f3 = 1.0_wp
1807       ENDIF
1808!
1809!--    Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1810!--    this calculation is obsolete, as r_canopy is not used below.
1811!--    To do: check for very dry soil -> r_canopy goes to infinity
1812       surf%r_canopy(m) = surf%r_canopy_min(m) /                               &
1813                              ( surf%lai(m) * f1 * f2 * f3 + 1.0E-20_wp )
1814!
1815!--    Third step: calculate bare soil resistance r_soil.
1816       m_min = surf%c_veg(m) * surf%m_wilt(nzb_soil,m) +                       &
1817                         ( 1.0_wp - surf%c_veg(m) ) * surf%m_res(nzb_soil,m)
1818
1819
1820       f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) /                       &
1821            ( surf%m_fc(nzb_soil,m) - m_min )
1822       f2 = MAX( f2, 1.0E-20_wp )
1823       f2 = MIN( f2, 1.0_wp     )
1824
1825       surf%r_soil(m) = surf%r_soil_min(m) / f2
1826       
1827!
1828!--    Calculate the maximum possible liquid water amount on plants and
1829!--    bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
1830!--    assumed, while paved surfaces might hold up 1 mm of water. The
1831!--    liquid water fraction for paved surfaces is calculated after
1832!--    Noilhan & Planton (1989), while the ECMWF formulation is used for
1833!--    vegetated surfaces and bare soils.
1834       IF ( surf%pavement_surface(m) )  THEN
1835          m_liq_max = m_max_depth * 5.0_wp
1836          surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq%var_1d(m) / m_liq_max)**0.67 )
1837       ELSE
1838          m_liq_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m)              &
1839                      + ( 1.0_wp - surf%c_veg(m) ) )
1840          surf%c_liq(m) = MIN( 1.0_wp, surf_m_liq%var_1d(m) / m_liq_max )
1841       ENDIF
1842!
1843!--    Calculate saturation specific humidity
1844       q_s = 0.622_wp * e_s / ( surface_pressure - e_s )
1845!
1846!--    In case of dewfall, set evapotranspiration to zero
1847!--    All super-saturated water is then removed from the air
1848       IF ( humidity  .AND.  q_s <= surf%qv1(m) )  THEN
1849          surf%r_canopy(m) = 0.0_wp
1850          surf%r_soil(m)   = 0.0_wp
1851       ENDIF
1852
1853!
1854!--    Calculate coefficients for the total evapotranspiration
1855!--    In case of water surface, set vegetation and soil fluxes to zero.
1856!--    For pavements, only evaporation of liquid water is possible.
1857       IF ( surf%water_surface(m) )  THEN
1858          f_qsws_veg  = 0.0_wp
1859          f_qsws_soil = 0.0_wp
1860          f_qsws_liq  = rho_lv / surf%r_a(m)
1861       ELSEIF ( surf%pavement_surface (m) )  THEN
1862          f_qsws_veg  = 0.0_wp
1863          f_qsws_soil = 0.0_wp
1864          f_qsws_liq  = rho_lv * surf%c_liq(m) / surf%r_a(m)
1865       ELSE
1866          f_qsws_veg  = rho_lv * surf%c_veg(m) *                               &
1867                            ( 1.0_wp        - surf%c_liq(m)    ) /             &
1868                            ( surf%r_a(m) + surf%r_canopy(m) )
1869          f_qsws_soil = rho_lv * (1.0_wp    - surf%c_veg(m)    ) /             &
1870                            ( surf%r_a(m) + surf%r_soil(m)   )
1871          f_qsws_liq  = rho_lv * surf%c_veg(m) * surf%c_liq(m)   /             &
1872                              surf%r_a(m)
1873       ENDIF
1874
1875       f_shf  = rho_cp / surf%r_a(m)
1876       f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1877!
1878!--    Calculate derivative of q_s for Taylor series expansion
1879       e_s_dt = e_s * ( 17.62_wp / ( surf_t_surface%var_1d(m) - 29.65_wp) -   &
1880                        17.62_wp*( surf_t_surface%var_1d(m) - 273.15_wp)      &
1881                       / ( surf_t_surface%var_1d(m) - 29.65_wp)**2 )
1882
1883       dq_s_dt = 0.622_wp * e_s_dt / ( surface_pressure - e_s_dt )
1884!
1885!--    Calculate net radiation radiation without longwave outgoing flux because
1886!--    it has a dependency on surface temperature and thus enters the prognostic
1887!--    equations directly
1888       surf%rad_net_l(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)              &
1889                           + surf%rad_lw_in(m)
1890!
1891!--    Calculate new skin temperature
1892       IF ( humidity )  THEN
1893!
1894!--       Numerator of the prognostic equation
1895          coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m)             &
1896                   * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)             &
1897                   + f_shf * surf%pt1(m) + f_qsws * ( surf%qv1(m) - q_s        &
1898                   + dq_s_dt * surf_t_surface%var_1d(m) ) + lambda_surface     &
1899                   * surf_t_soil%var_2d(nzb_soil,m)
1900
1901!
1902!--       Denominator of the prognostic equation
1903          coef_2 = surf%rad_lw_out_change_0(m) + f_qsws * dq_s_dt              &
1904                   + lambda_surface + f_shf / exn
1905       ELSE
1906!
1907!--       Numerator of the prognostic equation
1908          coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m)             &
1909                   * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)             &
1910                   + f_shf * surf%pt1(m)  + lambda_surface                     &
1911                   * surf_t_soil%var_2d(nzb_soil,m)
1912!
1913!--       Denominator of the prognostic equation
1914          coef_2 = surf%rad_lw_out_change_0(m) + lambda_surface + f_shf / exn
1915
1916       ENDIF
1917
1918
1919       tend = 0.0_wp
1920
1921!
1922!--    Implicit solution when the surface layer has no heat capacity,
1923!--    otherwise use RK3 scheme.
1924       surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *&
1925                          surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2&
1926                                             * dt_3d * tsc(2) ) 
1927
1928!
1929!--    Add RK3 term
1930       IF ( c_surface_tmp /= 0.0_wp )  THEN
1931
1932          surf_t_surface_p%var_1d(m) = surf_t_surface_p%var_1d(m) + dt_3d *    &
1933                                       tsc(3) * surf_tt_surface_m%var_1d(m)
1934
1935!
1936!--       Calculate true tendency
1937          tend = ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) -     &
1938                   dt_3d * tsc(3) * surf_tt_surface_m%var_1d(m)) / (dt_3d  * tsc(2))
1939!
1940!--       Calculate t_surface tendencies for the next Runge-Kutta step
1941          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1942             IF ( intermediate_timestep_count == 1 )  THEN
1943                surf_tt_surface_m%var_1d(m) = tend
1944             ELSEIF ( intermediate_timestep_count <                            &
1945                      intermediate_timestep_count_max )  THEN
1946                surf_tt_surface_m%var_1d(m) = -9.5625_wp * tend +              &
1947                                               5.3125_wp * surf_tt_surface_m%var_1d(m)
1948             ENDIF
1949          ENDIF
1950       ENDIF
1951
1952!
1953!--    In case of fast changes in the skin temperature, it is possible to
1954!--    update the radiative fluxes independently from the prescribed
1955!--    radiation call frequency. This effectively prevents oscillations,
1956!--    especially when setting skip_time_do_radiation /= 0. The threshold
1957!--    value of 0.2 used here is just a first guess. This method should be
1958!--    revised in the future as tests have shown that the threshold is
1959!--    often reached, when no oscillations would occur (causes immense
1960!--    computing time for the radiation code).
1961       IF ( ABS( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )       &
1962            > 0.2_wp  .AND. &
1963            unscheduled_radiation_calls )  THEN
1964          force_radiation_call_l = .TRUE.
1965       ENDIF
1966
1967
1968!        IF ( surf_t_surface_p%var_1d(m) < 270.0_wp .OR. surf_t_surface_p%var_1d(m) > 350.0_wp )  THEN
1969!           PRINT*, "myid: ", myid
1970!           PRINT*, "m: ", m
1971!           PRINT*, "i,j,k: ", i, j, k
1972!           PRINT*, "pt_p: ", surf_t_surface_p%var_1d(m)
1973!           PRINT*, "f_shf: ", f_shf
1974!           PRINT*, "f_qsws: ", f_qsws
1975!         ENDIF
1976
1977
1978!        pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn  !is actually no air temperature
1979       surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exn
1980
1981!
1982!--    Calculate fluxes
1983       surf%rad_net_l(m) = surf%rad_net_l(m) +                                 &
1984                            surf%rad_lw_out_change_0(m)                        &
1985                          * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)      &
1986                          - surf%rad_lw_out_change_0(m) * surf_t_surface_p%var_1d(m)
1987
1988       surf%rad_net(m) = surf%rad_net_l(m)
1989       surf%rad_lw_out(m) = surf%rad_lw_out(m) + surf%rad_lw_out_change_0(m) * &
1990                     ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )
1991
1992       surf%ghf(m) = lambda_surface * ( surf_t_surface_p%var_1d(m)             &
1993                                             - surf_t_soil%var_2d(nzb_soil,m) )
1994
1995       surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / cp
1996
1997       IF ( humidity )  THEN
1998          surf%qsws(m)  = - f_qsws * ( surf%qv1(m) - q_s + dq_s_dt             &
1999                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
2000                            surf_t_surface_p%var_1d(m) )
2001
2002          surf%qsws_veg(m)  = - f_qsws_veg  * ( surf%qv1(m) - q_s              &
2003                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
2004                              * surf_t_surface_p%var_1d(m) )
2005
2006          surf%qsws_soil(m) = - f_qsws_soil * ( surf%qv1(m) - q_s              &
2007                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
2008                              * surf_t_surface_p%var_1d(m) )
2009
2010          surf%qsws_liq(m)  = - f_qsws_liq  * ( surf%qv1(m) - q_s              &
2011                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
2012                              * surf_t_surface_p%var_1d(m) )
2013       ENDIF
2014
2015!
2016!--    Calculate the true surface resistance
2017       IF ( .NOT.  humidity )  THEN
2018          surf%r_s(m) = 1.0E10_wp
2019       ELSE
2020          surf%r_s(m) = - rho_lv * ( surf%qv1(m) - q_s + dq_s_dt               &
2021                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
2022                            surf_t_surface_p%var_1d(m) ) /                     &
2023                            (surf%qsws(m) + 1.0E-20)  - surf%r_a(m)
2024       ENDIF
2025
2026!
2027!--    Calculate change in liquid water reservoir due to dew fall or
2028!--    evaporation of liquid water
2029       IF ( humidity )  THEN
2030!
2031!--       If precipitation is activated, add rain water to qsws_liq
2032!--       and qsws_soil according the the vegetation coverage.
2033!--       precipitation_rate is given in mm.
2034          IF ( precipitation )  THEN
2035
2036!
2037!--          Add precipitation to liquid water reservoir, if possible.
2038!--          Otherwise, add the water to soil. In case of
2039!--          pavements, the exceeding water amount is implicitely removed
2040!--          as runoff as qsws_soil is then not used in the soil model
2041             IF ( surf_m_liq%var_1d(m) /= m_liq_max )  THEN
2042                surf%qsws_liq(m) = surf%qsws_liq(m)                            &
2043                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
2044                                 * hyrho(k+k_off)                              &
2045                                 * 0.001_wp * rho_l * l_v
2046             ELSE
2047                surf%qsws_soil(m) = surf%qsws_soil(m)                          &
2048                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
2049                                 * hyrho(k+k_off)                              &
2050                                 * 0.001_wp * rho_l * l_v
2051             ENDIF
2052
2053!--          Add precipitation to bare soil according to the bare soil
2054!--          coverage.
2055             surf%qsws_soil(m) = surf%qsws_soil(m) + ( 1.0_wp                  &
2056                               - surf%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
2057                               * hyrho(k+k_off)                                &
2058                               * 0.001_wp * rho_l * l_v
2059          ENDIF
2060
2061!
2062!--       If the air is saturated, check the reservoir water level
2063          IF ( surf%qsws(m) < 0.0_wp )  THEN
2064!
2065!--          Check if reservoir is full (avoid values > m_liq_max)
2066!--          In that case, qsws_liq goes to qsws_soil. In this
2067!--          case qsws_veg is zero anyway (because c_liq = 1),       
2068!--          so that tend is zero and no further check is needed
2069             IF ( surf_m_liq%var_1d(m) == m_liq_max )  THEN
2070                surf%qsws_soil(m) = surf%qsws_soil(m) + surf%qsws_liq(m)
2071
2072                surf%qsws_liq(m)  = 0.0_wp
2073             ENDIF
2074
2075!
2076!--          In case qsws_veg becomes negative (unphysical behavior),
2077!--          let the water enter the liquid water reservoir as dew on the
2078!--          plant
2079             IF ( surf%qsws_veg(m) < 0.0_wp )  THEN
2080                surf%qsws_liq(m) = surf%qsws_liq(m) + surf%qsws_veg(m)
2081                surf%qsws_veg(m) = 0.0_wp
2082             ENDIF
2083          ENDIF                   
2084 
2085          surf%qsws(m) = surf%qsws(m) / l_v
2086 
2087          tend = - surf%qsws_liq(m) * drho_l_lv
2088          surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d *              &
2089                                        ( tsc(2) * tend +                      &
2090                                          tsc(3) * surf_tm_liq_m%var_1d(m) )
2091!
2092!--       Check if reservoir is overfull -> reduce to maximum
2093!--       (conservation of water is violated here)
2094          surf_m_liq_p%var_1d(m) = MIN( surf_m_liq_p%var_1d(m),m_liq_max )
2095
2096!
2097!--       Check if reservoir is empty (avoid values < 0.0)
2098!--       (conservation of water is violated here)
2099          surf_m_liq_p%var_1d(m) = MAX( surf_m_liq_p%var_1d(m), 0.0_wp )
2100!
2101!--       Calculate m_liq tendencies for the next Runge-Kutta step
2102          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2103             IF ( intermediate_timestep_count == 1 )  THEN
2104                surf_tm_liq_m%var_1d(m) = tend
2105             ELSEIF ( intermediate_timestep_count <                            &
2106                      intermediate_timestep_count_max )  THEN
2107                surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend +                  &
2108                                           5.3125_wp * surf_tm_liq_m%var_1d(m)
2109             ENDIF
2110          ENDIF
2111
2112       ENDIF
2113
2114    ENDDO
2115
2116!
2117!-- Make a logical OR for all processes. Force radiation call if at
2118!-- least one processor reached the threshold change in skin temperature
2119    IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count       &
2120         == intermediate_timestep_count_max-1 )  THEN
2121#if defined( __parallel )
2122       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2123       CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,       &
2124                           1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
2125#else
2126       force_radiation_call = force_radiation_call_l
2127#endif
2128       force_radiation_call_l = .FALSE.
2129    ENDIF
2130
2131!
2132!-- Calculate surface specific humidity
2133    IF ( humidity )  THEN
2134       CALL calc_q_surface
2135    ENDIF
2136
2137!
2138!-- Calculate new roughness lengths (for water surfaces only)
2139    IF ( horizontal  .AND.  .NOT. constant_roughness )  CALL calc_z0_water_surface
2140
2141    CONTAINS
2142!------------------------------------------------------------------------------!
2143! Description:
2144! ------------
2145!> Calculation of specific humidity of the skin layer (surface). It is assumend
2146!> that the skin is always saturated.
2147!------------------------------------------------------------------------------!
2148    SUBROUTINE calc_q_surface
2149
2150       USE diagnostic_quantities_mod
2151
2152       IMPLICIT NONE
2153
2154       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
2155
2156       DO  m = 1, surf%ns
2157
2158          i   = surf%i(m)           
2159          j   = surf%j(m)
2160          k   = surf%k(m)
2161!
2162!--       Calculate water vapour pressure at saturation and convert to hPa
2163          e_s = 0.01_wp * magnus( surf_t_surface_p%var_1d(m) )                   
2164
2165!
2166!--       Calculate specific humidity at saturation
2167          q_s = 0.622_wp * e_s / ( surface_pressure - e_s )
2168
2169          resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) + 1E-5_wp )
2170
2171!
2172!--       Calculate specific humidity at surface
2173          IF ( cloud_physics )  THEN
2174             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
2175                                        ( 1.0_wp - resistance ) *              &
2176                                        ( q(k,j,i) - ql(k,j,i) )
2177          ELSE
2178             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
2179                                        ( 1.0_wp - resistance ) *              &
2180                                          q(k,j,i)
2181          ENDIF
2182!
2183!--       Update virtual potential temperature
2184          vpt(k+k_off,j+j_off,i+i_off) = pt(k+k_off,j+j_off,i+i_off) *         &
2185                     ( 1.0_wp + 0.61_wp * q(k+k_off,j+j_off,i+i_off) )
2186
2187       ENDDO
2188
2189    END SUBROUTINE calc_q_surface
2190
2191
2192
2193 END SUBROUTINE lsm_energy_balance
2194
2195
2196!------------------------------------------------------------------------------!
2197! Description:
2198! ------------
2199!> Header output for land surface model
2200!------------------------------------------------------------------------------!
2201    SUBROUTINE lsm_header ( io )
2202
2203
2204       IMPLICIT NONE
2205
2206       CHARACTER (LEN=86) ::  t_soil_chr          !< String for soil temperature profile
2207       CHARACTER (LEN=86) ::  roots_chr           !< String for root profile
2208       CHARACTER (LEN=86) ::  vertical_index_chr  !< String for the vertical index
2209       CHARACTER (LEN=86) ::  m_soil_chr          !< String for soil moisture
2210       CHARACTER (LEN=86) ::  soil_depth_chr      !< String for soil depth
2211       CHARACTER (LEN=10) ::  coor_chr            !< Temporary string
2212   
2213       INTEGER(iwp) ::  i                         !< Loop index over soil layers
2214 
2215       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
2216 
2217       t_soil_chr = ''
2218       m_soil_chr    = ''
2219       soil_depth_chr  = '' 
2220       roots_chr        = '' 
2221       vertical_index_chr   = ''
2222
2223       i = 1
2224       DO i = nzb_soil, nzt_soil
2225          WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i)
2226          t_soil_chr = TRIM( t_soil_chr ) // ' ' // TRIM( coor_chr )
2227
2228          WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i)
2229          m_soil_chr = TRIM( m_soil_chr ) // ' ' // TRIM( coor_chr )
2230
2231          WRITE (coor_chr,'(F10.2,7X)')  - zs(i)
2232          soil_depth_chr = TRIM( soil_depth_chr ) // ' '  // TRIM( coor_chr )
2233
2234          WRITE (coor_chr,'(F10.2,7X)')  root_fraction(i)
2235          roots_chr = TRIM( roots_chr ) // ' '  // TRIM( coor_chr )
2236
2237          WRITE (coor_chr,'(I10,7X)')  i
2238          vertical_index_chr = TRIM( vertical_index_chr ) // ' '  //           &
2239                               TRIM( coor_chr )
2240       ENDDO
2241
2242!
2243!--    Write land surface model header
2244       WRITE( io,  1 )
2245       IF ( conserve_water_content )  THEN
2246          WRITE( io, 2 )
2247       ELSE
2248          WRITE( io, 3 )
2249       ENDIF
2250
2251       IF ( vegetation_type_f%from_file )  THEN
2252          WRITE( io, 5 )
2253       ELSE
2254          WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ),        &
2255                         TRIM (soil_type_name(soil_type) )
2256       ENDIF
2257       WRITE( io, 6 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
2258                        TRIM( m_soil_chr ), TRIM( roots_chr ),                 &
2259                        TRIM( vertical_index_chr )
2260
22611   FORMAT (//' Land surface model information:'/                              &
2262              ' ------------------------------'/)
22632   FORMAT ('    --> Soil bottom is closed (water content is conserved',       &
2264            ', default)')
22653   FORMAT ('    --> Soil bottom is open (water content is not conserved)')         
22664   FORMAT ('    --> Land surface type  : ',A,/                                &
2267            '    --> Soil porosity type : ',A)
22685   FORMAT ('    --> Land surface type  : read from file' /                    &
2269            '    --> Soil porosity type : read from file' )
22706   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
2271            '       Height:        ',A,'  m'/                                  &
2272            '       Temperature:   ',A,'  K'/                                  &
2273            '       Moisture:      ',A,'  m**3/m**3'/                          &
2274            '       Root fraction: ',A,'  '/                                   &
2275            '       Grid point:    ',A)
2276
2277
2278    END SUBROUTINE lsm_header
2279
2280
2281!------------------------------------------------------------------------------!
2282! Description:
2283! ------------
2284!> Initialization of the land surface model
2285!------------------------------------------------------------------------------!
2286    SUBROUTINE lsm_init
2287   
2288       USE control_parameters,                                                 &
2289           ONLY:  message_string
2290   
2291       IMPLICIT NONE
2292
2293       INTEGER(iwp) ::  i                       !< running index
2294       INTEGER(iwp) ::  i_off                   !< index offset of surface element, seen from atmospheric grid point
2295       INTEGER(iwp) ::  j                       !< running index
2296       INTEGER(iwp) ::  j_off                   !< index offset of surface element, seen from atmospheric grid point
2297       INTEGER(iwp) ::  k                       !< running index
2298       INTEGER(iwp) ::  kn                      !< running index
2299       INTEGER(iwp) ::  ko                      !< running index
2300       INTEGER(iwp) ::  kroot                   !< running index
2301       INTEGER(iwp) ::  kzs                     !< running index
2302       INTEGER(iwp) ::  l                       !< running index surface facing
2303       INTEGER(iwp) ::  m                       !< running index
2304       INTEGER(iwp) ::  st                      !< soil-type index
2305       INTEGER(iwp) ::  n_soil_layers_total     !< temperature variable, stores the total number of soil layers + 4
2306       INTEGER(iwp) ::  n_surf                  !< number of surface types of given surface element
2307
2308       REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
2309
2310!
2311!--    Calculate Exner function
2312       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
2313!
2314!--    If no cloud physics is used, rho_surface has not been calculated before
2315       IF (  .NOT.  cloud_physics )  THEN
2316          CALL calc_mean_profile( pt, 4 )
2317          rho_surface = surface_pressure * 100.0_wp / ( r_d * hom(nzb+1,1,4,0) * exn )
2318       ENDIF
2319
2320!
2321!--    Calculate frequently used parameters
2322       rho_cp    = cp * rho_surface
2323       rd_d_rv   = r_d / r_v
2324       rho_lv    = rho_surface * l_v
2325       drho_l_lv = 1.0_wp / (rho_l * l_v)
2326
2327!
2328!--    Set initial values for prognostic quantities
2329!--    Horizontal surfaces
2330       tt_surface_h_m%var_1d = 0.0_wp
2331       tt_soil_h_m%var_2d    = 0.0_wp
2332       tm_soil_h_m%var_2d    = 0.0_wp
2333       tm_liq_h_m%var_1d  = 0.0_wp
2334       surf_lsm_h%c_liq = 0.0_wp
2335
2336       surf_lsm_h%ghf = 0.0_wp
2337
2338       surf_lsm_h%qsws_liq  = 0.0_wp
2339       surf_lsm_h%qsws_soil = 0.0_wp
2340       surf_lsm_h%qsws_veg  = 0.0_wp
2341
2342       surf_lsm_h%r_a        = 50.0_wp
2343       surf_lsm_h%r_s        = 50.0_wp
2344       surf_lsm_h%r_canopy   = 0.0_wp
2345       surf_lsm_h%r_soil     = 0.0_wp
2346!
2347!--    Do the same for vertical surfaces
2348       DO  l = 0, 3
2349          tt_surface_v_m(l)%var_1d = 0.0_wp
2350          tt_soil_v_m(l)%var_2d    = 0.0_wp
2351          tm_soil_v_m(l)%var_2d    = 0.0_wp
2352          tm_liq_v_m(l)%var_1d  = 0.0_wp
2353          surf_lsm_v(l)%c_liq = 0.0_wp
2354
2355          surf_lsm_v(l)%ghf = 0.0_wp
2356
2357          surf_lsm_v(l)%qsws_liq  = 0.0_wp
2358          surf_lsm_v(l)%qsws_soil = 0.0_wp
2359          surf_lsm_v(l)%qsws_veg  = 0.0_wp
2360
2361          surf_lsm_v(l)%r_a        = 50.0_wp
2362          surf_lsm_v(l)%r_s        = 50.0_wp
2363          surf_lsm_v(l)%r_canopy   = 0.0_wp
2364          surf_lsm_v(l)%r_soil     = 0.0_wp
2365       ENDDO
2366
2367!
2368!--    Set initial values for prognostic soil quantities
2369       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2370          t_soil_h%var_2d = 0.0_wp
2371          m_soil_h%var_2d = 0.0_wp
2372          m_liq_h%var_1d  = 0.0_wp
2373
2374          DO  l = 0, 3
2375             t_soil_v(l)%var_2d = 0.0_wp
2376             m_soil_v(l)%var_2d = 0.0_wp
2377             m_liq_v(l)%var_1d  = 0.0_wp
2378          ENDDO
2379       ENDIF
2380!
2381!--    Allocate 3D soil model arrays
2382!--    First, for horizontal surfaces
2383       ALLOCATE ( surf_lsm_h%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
2384       ALLOCATE ( surf_lsm_h%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2385       ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
2386       ALLOCATE ( surf_lsm_h%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns))
2387       ALLOCATE ( surf_lsm_h%l_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
2388       ALLOCATE ( surf_lsm_h%m_fc(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
2389       ALLOCATE ( surf_lsm_h%m_res(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
2390       ALLOCATE ( surf_lsm_h%m_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
2391       ALLOCATE ( surf_lsm_h%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_h%ns)      )
2392       ALLOCATE ( surf_lsm_h%n_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
2393       ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2394       ALLOCATE ( surf_lsm_h%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2395       ALLOCATE ( surf_lsm_h%root_fr(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
2396   
2397       surf_lsm_h%lambda_h     = 0.0_wp
2398!
2399!--    If required, allocate humidity-related variables for the soil model
2400       IF ( humidity )  THEN
2401          ALLOCATE ( surf_lsm_h%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2402          ALLOCATE ( surf_lsm_h%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  ) 
2403
2404          surf_lsm_h%lambda_w = 0.0_wp 
2405       ENDIF
2406!
2407!--    For vertical surfaces
2408       DO  l = 0, 3
2409          ALLOCATE ( surf_lsm_v(l)%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
2410          ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2411          ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
2412          ALLOCATE ( surf_lsm_v(l)%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns))
2413          ALLOCATE ( surf_lsm_v(l)%l_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
2414          ALLOCATE ( surf_lsm_v(l)%m_fc(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
2415          ALLOCATE ( surf_lsm_v(l)%m_res(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
2416          ALLOCATE ( surf_lsm_v(l)%m_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
2417          ALLOCATE ( surf_lsm_v(l)%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)      )
2418          ALLOCATE ( surf_lsm_v(l)%n_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
2419          ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) ) 
2420          ALLOCATE ( surf_lsm_v(l)%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) ) 
2421          ALLOCATE ( surf_lsm_v(l)%root_fr(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )
2422
2423          surf_lsm_v(l)%lambda_h     = 0.0_wp 
2424         
2425!
2426!--       If required, allocate humidity-related variables for the soil model
2427          IF ( humidity )  THEN
2428             ALLOCATE ( surf_lsm_v(l)%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2429             ALLOCATE ( surf_lsm_v(l)%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  ) 
2430
2431             surf_lsm_v(l)%lambda_w = 0.0_wp 
2432          ENDIF     
2433       ENDDO
2434!
2435!--    Allocate albedo type and emissivity for vegetation, water and pavement
2436!--    fraction.
2437!--    Set default values at each surface element.
2438       ALLOCATE ( surf_lsm_h%albedo_type(0:2,1:surf_lsm_h%ns) )
2439       ALLOCATE ( surf_lsm_h%emissivity(0:2,1:surf_lsm_h%ns) )
2440       surf_lsm_h%albedo_type = 0     
2441       surf_lsm_h%emissivity  = emissivity
2442       DO  l = 0, 3
2443          ALLOCATE ( surf_lsm_v(l)%albedo_type(0:2,1:surf_lsm_v(l)%ns) )
2444          ALLOCATE ( surf_lsm_v(l)%emissivity(0:2,1:surf_lsm_v(l)%ns)  )
2445          surf_lsm_v(l)%albedo_type = 0
2446          surf_lsm_v(l)%emissivity  = emissivity
2447       ENDDO
2448!
2449!--    Allocate arrays for relative surface fraction.
2450!--    0 - vegetation fraction, 2 - water fraction, 1 - pavement fraction
2451       ALLOCATE( surf_lsm_h%frac(0:2,1:surf_lsm_h%ns) )
2452       surf_lsm_h%frac = 0.0_wp
2453       DO  l = 0, 3
2454          ALLOCATE( surf_lsm_v(l)%frac(0:2,1:surf_lsm_v(l)%ns) )
2455          surf_lsm_v(l)%frac = 0.0_wp
2456       ENDDO
2457!
2458!--    For vertical walls only - allocate special flag indicating if any building is on
2459!--    top of any natural surfaces. Used for initialization only.
2460       DO  l = 0, 3
2461          ALLOCATE( surf_lsm_v(l)%building_covered(1:surf_lsm_v(l)%ns) )
2462       ENDDO
2463!
2464!--    Set flag parameter for the prescribed surface type depending on user
2465!--    input. Set surface fraction to 1 for the respective type.
2466       SELECT CASE ( TRIM( surface_type ) )
2467         
2468          CASE ( 'vegetation' )
2469         
2470             surf_lsm_h%vegetation_surface = .TRUE.
2471             surf_lsm_h%frac(ind_veg,:) = 1.0_wp
2472             DO  l = 0, 3
2473                surf_lsm_v(l)%vegetation_surface = .TRUE.
2474                surf_lsm_v(l)%frac(ind_veg,:) = 1.0_wp
2475             ENDDO
2476   
2477          CASE ( 'water' )
2478             
2479             surf_lsm_h%water_surface = .TRUE.
2480             surf_lsm_h%frac(ind_wat,:) = 1.0_wp
2481!
2482!--          Note, vertical water surface does not really make sense.
2483             DO  l = 0, 3 
2484                surf_lsm_v(l)%water_surface   = .TRUE.
2485                surf_lsm_v(l)%frac(ind_wat,:) = 1.0_wp
2486             ENDDO
2487
2488          CASE ( 'pavement' )
2489             
2490             surf_lsm_h%pavement_surface = .TRUE.
2491                surf_lsm_h%frac(ind_pav,:) = 1.0_wp
2492             DO  l = 0, 3
2493                surf_lsm_v(l)%pavement_surface   = .TRUE.
2494                surf_lsm_v(l)%frac(ind_pav,:) = 1.0_wp
2495             ENDDO
2496
2497          CASE ( 'netcdf' )
2498
2499             DO  m = 1, surf_lsm_h%ns
2500                i = surf_lsm_h%i(m)
2501                j = surf_lsm_h%j(m)
2502                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )        &
2503                   surf_lsm_h%vegetation_surface(m) = .TRUE.
2504                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )          &
2505                   surf_lsm_h%pavement_surface(m) = .TRUE.
2506                IF ( water_type_f%var(j,i)      /= water_type_f%fill )             &
2507                   surf_lsm_h%water_surface(m) = .TRUE.
2508             ENDDO
2509             DO  l = 0, 3
2510                DO  m = 1, surf_lsm_v(l)%ns
2511!
2512!--                Only for vertical surfaces. Check if natural walls at reference
2513!--                grid point are covered by any building. This case, problems
2514!--                with initialization will aris if index offsets are used.
2515!--                In order to deal with this, set special flag.
2516                   surf_lsm_v(l)%building_covered(m) = .FALSE.
2517                   IF ( building_type_f%from_file )  THEN
2518                      i = surf_lsm_v(l)%i(m) + surf_lsm_v(l)%ioff
2519                      j = surf_lsm_v(l)%j(m) + surf_lsm_v(l)%joff
2520                      IF ( building_type_f%var(j,i) /= 0 )                     &
2521                         surf_lsm_v(l)%building_covered(m) = .TRUE.
2522                   ENDIF
2523!
2524!--                Normally proceed with setting surface types.
2525                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,     &
2526                                            surf_lsm_v(l)%building_covered(m) )
2527                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,     &
2528                                            surf_lsm_v(l)%building_covered(m) )
2529                   IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
2530                      surf_lsm_v(l)%vegetation_surface(m) = .TRUE.
2531                   IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )   &
2532                      surf_lsm_v(l)%pavement_surface(m) = .TRUE.
2533                   IF ( water_type_f%var(j,i)      /= water_type_f%fill )      &
2534                      surf_lsm_v(l)%water_surface(m) = .TRUE.
2535                ENDDO
2536             ENDDO
2537
2538       END SELECT
2539!
2540!--    In case of netcdf input file, further initialize surface fractions.
2541!--    At the moment only 1 surface is given at a location, so that the fraction
2542!--    is either 0 or 1. This will be revised later.
2543       IF ( input_pids_static )  THEN
2544          DO  m = 1, surf_lsm_h%ns
2545             i = surf_lsm_h%i(m)
2546             j = surf_lsm_h%j(m)
2547!
2548!--          0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction             
2549             surf_lsm_h%frac(ind_veg,m) = surface_fraction_f%frac(ind_veg,j,i)         
2550             surf_lsm_h%frac(ind_pav,m) = surface_fraction_f%frac(ind_pav,j,i)       
2551             surf_lsm_h%frac(ind_wat,m) = surface_fraction_f%frac(ind_wat,j,i)
2552
2553          ENDDO
2554          DO  l = 0, 3
2555             DO  m = 1, surf_lsm_v(l)%ns
2556                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
2557                                                surf_lsm_v(l)%building_covered(m) ) 
2558                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
2559                                                surf_lsm_v(l)%building_covered(m) ) 
2560!
2561!--             0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction       
2562                surf_lsm_v(l)%frac(ind_veg,m) = surface_fraction_f%frac(ind_veg,j,i)         
2563                surf_lsm_v(l)%frac(ind_pav,m) = surface_fraction_f%frac(ind_pav,j,i)       
2564                surf_lsm_v(l)%frac(ind_wat,m) = surface_fraction_f%frac(ind_wat,j,i)
2565
2566             ENDDO
2567          ENDDO
2568       ENDIF
2569!
2570!--    Level 1, initialization of soil parameters.
2571!--    It is possible to overwrite each parameter by setting the respecticy
2572!--    NAMELIST variable to a value /= 9999999.9.
2573       IF ( soil_type /= 0 )  THEN 
2574 
2575          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
2576             alpha_vangenuchten = soil_pars(0,soil_type)
2577          ENDIF
2578
2579          IF ( l_vangenuchten == 9999999.9_wp )  THEN
2580             l_vangenuchten = soil_pars(1,soil_type)
2581          ENDIF
2582
2583          IF ( n_vangenuchten == 9999999.9_wp )  THEN
2584             n_vangenuchten = soil_pars(2,soil_type)           
2585          ENDIF
2586
2587          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
2588             hydraulic_conductivity = soil_pars(3,soil_type)           
2589          ENDIF
2590
2591          IF ( saturation_moisture == 9999999.9_wp )  THEN
2592             saturation_moisture = soil_pars(4,soil_type)           
2593          ENDIF
2594
2595          IF ( field_capacity == 9999999.9_wp )  THEN
2596             field_capacity = soil_pars(5,soil_type)           
2597          ENDIF
2598
2599          IF ( wilting_point == 9999999.9_wp )  THEN
2600             wilting_point = soil_pars(6,soil_type)           
2601          ENDIF
2602
2603          IF ( residual_moisture == 9999999.9_wp )  THEN
2604             residual_moisture = soil_pars(7,soil_type)       
2605          ENDIF
2606
2607       ENDIF
2608!
2609!--    Map values to the respective 2D/3D arrays
2610!--    Horizontal surfaces
2611       surf_lsm_h%alpha_vg      = alpha_vangenuchten
2612       surf_lsm_h%l_vg          = l_vangenuchten
2613       surf_lsm_h%n_vg          = n_vangenuchten
2614       surf_lsm_h%gamma_w_sat   = hydraulic_conductivity
2615       surf_lsm_h%m_sat         = saturation_moisture
2616       surf_lsm_h%m_fc          = field_capacity
2617       surf_lsm_h%m_wilt        = wilting_point
2618       surf_lsm_h%m_res         = residual_moisture
2619       surf_lsm_h%r_soil_min    = min_soil_resistance
2620!
2621!--    Vertical surfaces
2622       DO  l = 0, 3
2623          surf_lsm_v(l)%alpha_vg      = alpha_vangenuchten
2624          surf_lsm_v(l)%l_vg          = l_vangenuchten
2625          surf_lsm_v(l)%n_vg          = n_vangenuchten
2626          surf_lsm_v(l)%gamma_w_sat   = hydraulic_conductivity
2627          surf_lsm_v(l)%m_sat         = saturation_moisture
2628          surf_lsm_v(l)%m_fc          = field_capacity
2629          surf_lsm_v(l)%m_wilt        = wilting_point
2630          surf_lsm_v(l)%m_res         = residual_moisture
2631          surf_lsm_v(l)%r_soil_min    = min_soil_resistance
2632       ENDDO
2633!
2634!--    Level 2, initialization of soil parameters via soil_type read from file.
2635!--    Soil parameters are initialized for each (y,x)-grid point
2636!--    individually using default paramter settings according to the given
2637!--    soil type.
2638       IF ( soil_type_f%from_file )  THEN
2639!
2640!--       Level of detail = 1, i.e. a homogeneous soil distribution along the
2641!--       vertical dimension is assumed.
2642          IF ( soil_type_f%lod == 1 )  THEN
2643!
2644!--          Horizontal surfaces
2645             DO  m = 1, surf_lsm_h%ns
2646                i = surf_lsm_h%i(m)
2647                j = surf_lsm_h%j(m)
2648             
2649                st = soil_type_f%var_2d(j,i)
2650                IF ( st /= soil_type_f%fill )  THEN
2651                   surf_lsm_h%alpha_vg(:,m)    = soil_pars(0,st)
2652                   surf_lsm_h%l_vg(:,m)        = soil_pars(1,st)
2653                   surf_lsm_h%n_vg(:,m)        = soil_pars(2,st)
2654                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars(3,st)
2655                   surf_lsm_h%m_sat(:,m)       = soil_pars(4,st)
2656                   surf_lsm_h%m_fc(:,m)        = soil_pars(5,st)
2657                   surf_lsm_h%m_wilt(:,m)      = soil_pars(6,st)
2658                   surf_lsm_h%m_res(:,m)       = soil_pars(7,st)
2659                ENDIF
2660             ENDDO
2661!
2662!--          Vertical surfaces ( assumes the soil type given at respective (x,y)
2663             DO  l = 0, 3
2664                DO  m = 1, surf_lsm_v(l)%ns
2665                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2666                                                   surf_lsm_v(l)%building_covered(m) ) 
2667                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2668                                                   surf_lsm_v(l)%building_covered(m) ) 
2669
2670                   st = soil_type_f%var_2d(j,i)
2671                   IF ( st /= soil_type_f%fill )  THEN
2672                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars(0,st)
2673                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars(1,st)
2674                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars(2,st)
2675                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars(3,st)
2676                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars(4,st)
2677                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars(5,st)
2678                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars(6,st)
2679                      surf_lsm_v(l)%m_res(:,m)       = soil_pars(7,st)
2680                   ENDIF
2681                ENDDO
2682             ENDDO
2683!
2684!--       Level of detail = 2, i.e. soil type and thus the soil parameters
2685!--       can be heterogeneous along the vertical dimension.
2686          ELSE
2687!
2688!--          Horizontal surfaces
2689             DO  m = 1, surf_lsm_h%ns
2690                i = surf_lsm_h%i(m)
2691                j = surf_lsm_h%j(m)
2692             
2693                DO  k = nzb_soil, nzt_soil
2694                   st = soil_type_f%var_3d(k,j,i)
2695                   IF ( st /= soil_type_f%fill )  THEN
2696                      surf_lsm_h%alpha_vg(k,m)    = soil_pars(0,st)
2697                      surf_lsm_h%l_vg(k,m)        = soil_pars(1,st)
2698                      surf_lsm_h%n_vg(k,m)        = soil_pars(2,st)
2699                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars(3,st)
2700                      surf_lsm_h%m_sat(k,m)       = soil_pars(4,st)
2701                      surf_lsm_h%m_fc(k,m)        = soil_pars(5,st)
2702                      surf_lsm_h%m_wilt(k,m)      = soil_pars(6,st)
2703                      surf_lsm_h%m_res(k,m)       = soil_pars(7,st)
2704                   ENDIF
2705                ENDDO
2706             ENDDO
2707!
2708!--          Vertical surfaces ( assumes the soil type given at respective (x,y)
2709             DO  l = 0, 3
2710                DO  m = 1, surf_lsm_v(l)%ns
2711                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2712                                                   surf_lsm_v(l)%building_covered(m) ) 
2713                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2714                                                   surf_lsm_v(l)%building_covered(m) ) 
2715
2716                   DO  k = nzb_soil, nzt_soil
2717                      st = soil_type_f%var_3d(k,j,i)
2718                      IF ( st /= soil_type_f%fill )  THEN
2719                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars(0,st)
2720                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars(1,st)
2721                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars(2,st)
2722                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars(3,st)
2723                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars(4,st)
2724                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars(5,st)
2725                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars(6,st)
2726                         surf_lsm_v(l)%m_res(k,m)       = soil_pars(7,st)
2727                      ENDIF
2728                   ENDDO
2729                ENDDO
2730             ENDDO
2731          ENDIF
2732       ENDIF
2733!
2734!--    Level 3, initialization of single soil parameters at single z,x,y
2735!--    position via soil_pars read from file.
2736       IF ( soil_pars_f%from_file )  THEN
2737!
2738!--       Level of detail = 1, i.e. a homogeneous vertical distribution of soil
2739!--       parameters is assumed.
2740!--       Horizontal surfaces
2741          IF ( soil_pars_f%lod == 1 )  THEN
2742!
2743!--          Horizontal surfaces
2744             DO  m = 1, surf_lsm_h%ns
2745                i = surf_lsm_h%i(m)
2746                j = surf_lsm_h%j(m)
2747
2748                IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )              &
2749                   surf_lsm_h%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
2750                IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )              &
2751                   surf_lsm_h%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
2752                IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )              &
2753                   surf_lsm_h%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
2754                IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )              &
2755                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
2756                IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )              &
2757                   surf_lsm_h%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
2758                IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )              &
2759                   surf_lsm_h%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
2760                IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )              &
2761                   surf_lsm_h%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
2762                IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )              &
2763                   surf_lsm_h%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)
2764
2765             ENDDO
2766!
2767!--          Vertical surfaces
2768             DO  l = 0, 3
2769                DO  m = 1, surf_lsm_v(l)%ns
2770                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2771                                                   surf_lsm_v(l)%building_covered(m) ) 
2772                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2773                                                   surf_lsm_v(l)%building_covered(m) ) 
2774
2775                   IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )           &
2776                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
2777                   IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )           &
2778                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
2779                   IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )           &
2780                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
2781                   IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )           &
2782                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
2783                   IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )           &
2784                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
2785                   IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )           &
2786                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
2787                   IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )           &
2788                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
2789                   IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )           &
2790                      surf_lsm_v(l)%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)
2791
2792                ENDDO
2793             ENDDO
2794!
2795!--       Level of detail = 2, i.e. soil parameters can be set at each soil
2796!--       layer individually.
2797          ELSE
2798!
2799!--          Horizontal surfaces
2800             DO  m = 1, surf_lsm_h%ns
2801                i = surf_lsm_h%i(m)
2802                j = surf_lsm_h%j(m)
2803
2804                DO  k = nzb_soil, nzt_soil
2805                   IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
2806                      surf_lsm_h%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
2807                   IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
2808                      surf_lsm_h%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
2809                   IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
2810                      surf_lsm_h%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
2811                   IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
2812                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
2813                   IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
2814                      surf_lsm_h%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
2815                   IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
2816                      surf_lsm_h%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
2817                   IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
2818                      surf_lsm_h%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
2819                   IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
2820                      surf_lsm_h%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
2821                ENDDO
2822
2823             ENDDO
2824!
2825!--          Vertical surfaces
2826             DO  l = 0, 3
2827                DO  m = 1, surf_lsm_v(l)%ns
2828                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2829                                                   surf_lsm_v(l)%building_covered(m) ) 
2830                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2831                                                   surf_lsm_v(l)%building_covered(m) ) 
2832
2833                   DO  k = nzb_soil, nzt_soil
2834                      IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
2835                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
2836                      IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
2837                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
2838                      IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
2839                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
2840                      IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
2841                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
2842                      IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
2843                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
2844                      IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
2845                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
2846                      IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
2847                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
2848                      IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
2849                         surf_lsm_v(l)%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
2850                   ENDDO
2851
2852                ENDDO
2853             ENDDO
2854
2855          ENDIF
2856       ENDIF
2857
2858!
2859!--    Level 1, initialization of vegetation parameters. A horizontally
2860!--    homogeneous distribution is assumed here.
2861       IF ( vegetation_type /= 0 )  THEN
2862
2863          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
2864             min_canopy_resistance = vegetation_pars(ind_v_rc_min,vegetation_type)
2865          ENDIF
2866
2867          IF ( leaf_area_index == 9999999.9_wp )  THEN
2868             leaf_area_index = vegetation_pars(ind_v_rc_lai,vegetation_type)         
2869          ENDIF
2870
2871          IF ( vegetation_coverage == 9999999.9_wp )  THEN
2872             vegetation_coverage = vegetation_pars(ind_v_c_veg,vegetation_type)     
2873          ENDIF
2874
2875          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
2876              canopy_resistance_coefficient= vegetation_pars(ind_v_gd,vegetation_type)     
2877          ENDIF
2878
2879          IF ( z0_vegetation == 9999999.9_wp )  THEN
2880             z0_vegetation  = vegetation_pars(ind_v_z0,vegetation_type) 
2881          ENDIF
2882
2883          IF ( z0h_vegetation == 9999999.9_wp )  THEN
2884             z0h_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
2885          ENDIF   
2886         
2887          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
2888             lambda_surface_stable = vegetation_pars(ind_v_lambda_s,vegetation_type) 
2889          ENDIF
2890
2891          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
2892             lambda_surface_unstable = vegetation_pars(ind_v_lambda_u,vegetation_type)           
2893          ENDIF
2894
2895          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
2896             f_shortwave_incoming = vegetation_pars(ind_v_f_sw_in,vegetation_type)       
2897          ENDIF
2898
2899          IF ( c_surface == 9999999.9_wp )  THEN
2900             c_surface = vegetation_pars(ind_v_c_surf,vegetation_type)       
2901          ENDIF
2902
2903          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
2904             albedo_type = INT(vegetation_pars(ind_v_at,vegetation_type))       
2905          ENDIF
2906   
2907          IF ( emissivity == 9999999.9_wp )  THEN
2908             emissivity = vegetation_pars(ind_v_emis,vegetation_type)     
2909          ENDIF
2910
2911       ENDIF
2912!
2913!--    Map values onto horizontal elemements
2914       DO  m = 1, surf_lsm_h%ns
2915          IF ( surf_lsm_h%vegetation_surface(m) )  THEN
2916             surf_lsm_h%r_canopy_min(m)     = min_canopy_resistance
2917             surf_lsm_h%lai(m)              = leaf_area_index
2918             surf_lsm_h%c_veg(m)            = vegetation_coverage
2919             surf_lsm_h%g_d(m)              = canopy_resistance_coefficient
2920             surf_lsm_h%z0(m)               = z0_vegetation
2921             surf_lsm_h%z0h(m)              = z0h_vegetation
2922             surf_lsm_h%z0q(m)              = z0h_vegetation
2923             surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
2924             surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
2925             surf_lsm_h%f_sw_in(m)          = f_shortwave_incoming
2926             surf_lsm_h%c_surface(m)        = c_surface
2927             surf_lsm_h%albedo_type(ind_veg,m) = albedo_type
2928             surf_lsm_h%emissivity(ind_veg,m)  = emissivity
2929          ELSE
2930             surf_lsm_h%lai(m)   = 0.0_wp
2931             surf_lsm_h%c_veg(m) = 0.0_wp
2932             surf_lsm_h%g_d(m)   = 0.0_wp
2933          ENDIF
2934 
2935       ENDDO
2936!
2937!--    Map values onto vertical elements, even though this does not make
2938!--    much sense.
2939       DO  l = 0, 3
2940          DO  m = 1, surf_lsm_v(l)%ns
2941             IF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
2942                surf_lsm_v(l)%r_canopy_min(m)     = min_canopy_resistance
2943                surf_lsm_v(l)%lai(m)              = leaf_area_index
2944                surf_lsm_v(l)%c_veg(m)            = vegetation_coverage
2945                surf_lsm_v(l)%g_d(m)              = canopy_resistance_coefficient
2946                surf_lsm_v(l)%z0(m)               = z0_vegetation
2947                surf_lsm_v(l)%z0h(m)              = z0h_vegetation
2948                surf_lsm_v(l)%z0q(m)              = z0h_vegetation
2949                surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
2950                surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
2951                surf_lsm_v(l)%f_sw_in(m)          = f_shortwave_incoming
2952                surf_lsm_v(l)%c_surface(m)        = c_surface
2953                surf_lsm_v(l)%albedo_type(ind_veg,m) = albedo_type
2954                surf_lsm_v(l)%emissivity(ind_veg,m)  = emissivity
2955             ELSE
2956                surf_lsm_v(l)%lai(m)   = 0.0_wp
2957                surf_lsm_v(l)%c_veg(m) = 0.0_wp
2958                surf_lsm_v(l)%g_d(m)   = 0.0_wp
2959             ENDIF
2960          ENDDO
2961       ENDDO
2962
2963!
2964!--    Level 2, initialization of vegation parameters via vegetation_type read
2965!--    from file. Vegetation parameters are initialized for each (y,x)-grid point
2966!--    individually using default paramter settings according to the given
2967!--    vegetation type.
2968       IF ( vegetation_type_f%from_file )  THEN
2969!
2970!--       Horizontal surfaces
2971          DO  m = 1, surf_lsm_h%ns
2972             i = surf_lsm_h%i(m)
2973             j = surf_lsm_h%j(m)
2974             
2975             st = vegetation_type_f%var(j,i)
2976             IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
2977                surf_lsm_h%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
2978                surf_lsm_h%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
2979                surf_lsm_h%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
2980                surf_lsm_h%g_d(m)              = vegetation_pars(ind_v_gd,st)
2981                surf_lsm_h%z0(m)               = vegetation_pars(ind_v_z0,st)
2982                surf_lsm_h%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
2983                surf_lsm_h%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
2984                surf_lsm_h%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
2985                surf_lsm_h%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
2986                surf_lsm_h%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
2987                surf_lsm_h%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
2988                surf_lsm_h%albedo_type(ind_veg,m) = INT( vegetation_pars(ind_v_at,st) )
2989                surf_lsm_h%emissivity(ind_veg,m)  = vegetation_pars(ind_v_emis,st)
2990             ENDIF
2991          ENDDO
2992!
2993!--       Vertical surfaces
2994          DO  l = 0, 3
2995             DO  m = 1, surf_lsm_v(l)%ns
2996                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2997                                                surf_lsm_v(l)%building_covered(m) ) 
2998                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2999                                                   surf_lsm_v(l)%building_covered(m) ) 
3000             
3001                st = vegetation_type_f%var(j,i)
3002                IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
3003                   surf_lsm_v(l)%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
3004                   surf_lsm_v(l)%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
3005                   surf_lsm_v(l)%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
3006                   surf_lsm_v(l)%g_d(m)              = vegetation_pars(ind_v_gd,st)
3007                   surf_lsm_v(l)%z0(m)               = vegetation_pars(ind_v_z0,st)
3008                   surf_lsm_v(l)%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
3009                   surf_lsm_v(l)%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
3010                   surf_lsm_v(l)%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
3011                   surf_lsm_v(l)%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
3012                   surf_lsm_v(l)%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
3013                   surf_lsm_v(l)%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
3014                   surf_lsm_v(l)%albedo_type(ind_veg,m) = INT( vegetation_pars(ind_v_at,st) )
3015                   surf_lsm_v(l)%emissivity(ind_veg,m)  = vegetation_pars(ind_v_emis,st)
3016                ENDIF
3017             ENDDO
3018          ENDDO
3019       ENDIF
3020!
3021!--    Level 3, initialization of vegation parameters at single (x,y)
3022!--    position via vegetation_pars read from file.
3023       IF ( vegetation_pars_f%from_file )  THEN
3024!
3025!--       Horizontal surfaces
3026          DO  m = 1, surf_lsm_h%ns
3027
3028             i = surf_lsm_h%i(m)
3029             j = surf_lsm_h%j(m)
3030!
3031!--          If surface element is not a vegetation surface and any value in
3032!--          vegetation_pars is given, neglect this information and give an
3033!--          informative message that this value will not be used.   
3034             IF ( .NOT. surf_lsm_h%vegetation_surface(m)  .AND.                &
3035                   ANY( vegetation_pars_f%pars_xy(:,j,i) /=                    &
3036                   vegetation_pars_f%fill ) )  THEN
3037                WRITE( message_string, * )                                     &
3038                                 'surface element at grid point (j,i) = (',    &
3039                                 j, i, ') is not a vegation surface, ',        &
3040                                 'so that information given in ',              &
3041                                 'vegetation_pars at this point is neglected.' 
3042                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
3043             ELSE
3044
3045                IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=            &
3046                     vegetation_pars_f%fill )                                  &
3047                   surf_lsm_h%r_canopy_min(m)  =                               &
3048                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
3049                IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=            &
3050                     vegetation_pars_f%fill )                                  &
3051                   surf_lsm_h%lai(m)           =                               &
3052                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
3053                IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=             &
3054                     vegetation_pars_f%fill )                                  &
3055                   surf_lsm_h%c_veg(m)         =                               &
3056                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
3057                IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=                &
3058                     vegetation_pars_f%fill )                                  &
3059                   surf_lsm_h%g_d(m)           =                               &
3060                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
3061                IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=                &
3062                     vegetation_pars_f%fill )                                  &
3063                   surf_lsm_h%z0(m)            =                               &
3064                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
3065                IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=              &
3066                     vegetation_pars_f%fill )  THEN
3067                   surf_lsm_h%z0h(m)           =                               &
3068                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3069                   surf_lsm_h%z0q(m)           =                               &
3070                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3071                ENDIF
3072                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=          &
3073                     vegetation_pars_f%fill )                                  &
3074                   surf_lsm_h%lambda_surface_s(m) =                            &
3075                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
3076                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=          &
3077                     vegetation_pars_f%fill )                                  &
3078                   surf_lsm_h%lambda_surface_u(m) =                            &
3079                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
3080                IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=           &
3081                     vegetation_pars_f%fill )                                  &
3082                   surf_lsm_h%f_sw_in(m)          =                            &
3083                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
3084                IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=            &
3085                     vegetation_pars_f%fill )                                  &
3086                   surf_lsm_h%c_surface(m)        =                            &
3087                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
3088                IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=                &
3089                     vegetation_pars_f%fill )                                  &
3090                   surf_lsm_h%albedo_type(ind_veg,m) =                         &
3091                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
3092                IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=              &
3093                     vegetation_pars_f%fill )                                  &
3094                   surf_lsm_h%emissivity(ind_veg,m)  =                         &
3095                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
3096             ENDIF
3097          ENDDO
3098!
3099!--       Vertical surfaces
3100          DO  l = 0, 3
3101             DO  m = 1, surf_lsm_v(l)%ns
3102                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3103                                                surf_lsm_v(l)%building_covered(m) ) 
3104                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3105                                                surf_lsm_v(l)%building_covered(m) ) 
3106!
3107!--             If surface element is not a vegetation surface and any value in
3108!--             vegetation_pars is given, neglect this information and give an
3109!--             informative message that this value will not be used.   
3110                IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m)  .AND.          &
3111                      ANY( vegetation_pars_f%pars_xy(:,j,i) /=                 &
3112                      vegetation_pars_f%fill ) )  THEN
3113                   WRITE( message_string, * )                                  &
3114                                 'surface element at grid point (j,i) = (',    &
3115                                 j, i, ') is not a vegation surface, ',        &
3116                                 'so that information given in ',              &
3117                                 'vegetation_pars at this point is neglected.' 
3118                   CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
3119                ELSE
3120
3121                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=         &
3122                        vegetation_pars_f%fill )                               &
3123                      surf_lsm_v(l)%r_canopy_min(m)  =                         &
3124                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
3125                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=         &
3126                        vegetation_pars_f%fill )                               &
3127                      surf_lsm_v(l)%lai(m)           =                         &
3128                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
3129                   IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=          &
3130                        vegetation_pars_f%fill )                               &
3131                      surf_lsm_v(l)%c_veg(m)         =                         &
3132                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
3133                   IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=             &
3134                        vegetation_pars_f%fill )                               &
3135                     surf_lsm_v(l)%g_d(m)            =                         &
3136                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
3137                   IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=             &
3138                        vegetation_pars_f%fill )                               &
3139                      surf_lsm_v(l)%z0(m)            =                         &
3140                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
3141                   IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=           &
3142                        vegetation_pars_f%fill )  THEN
3143                      surf_lsm_v(l)%z0h(m)           =                         &
3144                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3145                      surf_lsm_v(l)%z0q(m)           =                         &
3146                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3147                   ENDIF
3148                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=       &
3149                        vegetation_pars_f%fill )                               &
3150                      surf_lsm_v(l)%lambda_surface_s(m)  =                     &
3151                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
3152                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=       &
3153                        vegetation_pars_f%fill )                               &
3154                      surf_lsm_v(l)%lambda_surface_u(m)  =                     &
3155                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
3156                   IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=        &
3157                        vegetation_pars_f%fill )                               &
3158                      surf_lsm_v(l)%f_sw_in(m)           =                     &
3159                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
3160                   IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=         &
3161                        vegetation_pars_f%fill )                               &
3162                      surf_lsm_v(l)%c_surface(m)         =                     &
3163                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
3164                   IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=             &
3165                        vegetation_pars_f%fill )                               &
3166                      surf_lsm_v(l)%albedo_type(ind_veg,m) =                   &
3167                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
3168                   IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=           &
3169                        vegetation_pars_f%fill )                               &
3170                      surf_lsm_v(l)%emissivity(ind_veg,m)  =                   &
3171                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
3172                ENDIF
3173
3174             ENDDO
3175          ENDDO
3176       ENDIF
3177
3178!
3179!--    Level 1, initialization of water parameters. A horizontally
3180!--    homogeneous distribution is assumed here.
3181       IF ( water_type /= 0 )  THEN
3182
3183          IF ( water_temperature == 9999999.9_wp )  THEN
3184             water_temperature = water_pars(ind_w_temp,water_type)       
3185          ENDIF
3186
3187          IF ( z0_water == 9999999.9_wp )  THEN
3188             z0_water = water_pars(ind_w_z0,water_type)       
3189          ENDIF       
3190
3191          IF ( z0h_water == 9999999.9_wp )  THEN
3192             z0h_water = water_pars(ind_w_z0h,water_type)       
3193          ENDIF 
3194
3195          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
3196             albedo_type = INT(water_pars(ind_w_at,water_type))       
3197          ENDIF
3198   
3199          IF ( emissivity == 9999999.9_wp )  THEN
3200             emissivity = water_pars(ind_w_emis,water_type)       
3201          ENDIF
3202
3203       ENDIF
3204!
3205!--    Map values onto horizontal elemements
3206       DO  m = 1, surf_lsm_h%ns
3207          IF ( surf_lsm_h%water_surface(m) )  THEN
3208             IF ( TRIM( initializing_actions ) /= 'read_restart_data' )        &
3209                t_soil_h%var_2d(:,m)           = water_temperature
3210             surf_lsm_h%z0(m)               = z0_water
3211             surf_lsm_h%z0h(m)              = z0h_water
3212             surf_lsm_h%z0q(m)              = z0h_water
3213             surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
3214             surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp               
3215             surf_lsm_h%c_surface(m)        = 0.0_wp
3216             surf_lsm_h%albedo_type(ind_wat,m) = albedo_type
3217             surf_lsm_h%emissivity(ind_wat,m)  = emissivity
3218          ENDIF
3219       ENDDO
3220!
3221!--    Map values onto vertical elements, even though this does not make
3222!--    much sense.
3223       DO  l = 0, 3
3224          DO  m = 1, surf_lsm_v(l)%ns
3225             IF ( surf_lsm_v(l)%water_surface(m) )  THEN
3226                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
3227                   t_soil_v(l)%var_2d(:,m)           = water_temperature
3228                surf_lsm_v(l)%z0(m)               = z0_water
3229                surf_lsm_v(l)%z0h(m)              = z0h_water
3230                surf_lsm_v(l)%z0q(m)              = z0h_water
3231                surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp
3232                surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp               
3233                surf_lsm_v(l)%c_surface(m)        = 0.0_wp
3234                surf_lsm_v(l)%albedo_type(ind_wat,m) = albedo_type
3235                surf_lsm_v(l)%emissivity(ind_wat,m)  = emissivity
3236             ENDIF
3237          ENDDO
3238       ENDDO
3239!
3240!
3241!--    Level 2, initialization of water parameters via water_type read
3242!--    from file. Water surfaces are initialized for each (y,x)-grid point
3243!--    individually using default paramter settings according to the given
3244!--    water type.
3245!--    Note, parameter 3/4 of water_pars are albedo and emissivity,
3246!--    whereas paramter 3/4 of water_pars_f are heat conductivities!
3247       IF ( water_type_f%from_file )  THEN
3248!
3249!--       Horizontal surfaces
3250          DO  m = 1, surf_lsm_h%ns
3251             i = surf_lsm_h%i(m)
3252             j = surf_lsm_h%j(m)
3253             
3254             st = water_type_f%var(j,i)
3255             IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
3256                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
3257                   t_soil_h%var_2d(:,m) = water_pars(ind_w_temp,st)
3258                surf_lsm_h%z0(m)     = water_pars(ind_w_z0,st)
3259                surf_lsm_h%z0h(m)    = water_pars(ind_w_z0h,st)
3260                surf_lsm_h%z0q(m)    = water_pars(ind_w_z0h,st)
3261                surf_lsm_h%lambda_surface_s(m) = water_pars(ind_w_lambda_s,st)
3262                surf_lsm_h%lambda_surface_u(m) = water_pars(ind_w_lambda_u,st)             
3263                surf_lsm_h%c_surface(m)        = 0.0_wp
3264                surf_lsm_h%albedo_type(ind_wat,m) = INT( water_pars(ind_w_at,st) )
3265                surf_lsm_h%emissivity(ind_wat,m)  = water_pars(ind_w_emis,st)
3266             ENDIF
3267          ENDDO
3268!
3269!--       Vertical surfaces
3270          DO  l = 0, 3
3271             DO  m = 1, surf_lsm_v(l)%ns
3272                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3273                                                surf_lsm_v(l)%building_covered(m) ) 
3274                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3275                                                surf_lsm_v(l)%building_covered(m) ) 
3276             
3277                st = water_type_f%var(j,i)
3278                IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
3279                   IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  &
3280                      t_soil_v(l)%var_2d(:,m) = water_pars(ind_w_temp,st)
3281                   surf_lsm_v(l)%z0(m)     = water_pars(ind_w_z0,st)
3282                   surf_lsm_v(l)%z0h(m)    = water_pars(ind_w_z0h,st)
3283                   surf_lsm_v(l)%z0q(m)    = water_pars(ind_w_z0h,st)
3284                   surf_lsm_v(l)%lambda_surface_s(m) =                         &
3285                                                   water_pars(ind_w_lambda_s,st)
3286                   surf_lsm_v(l)%lambda_surface_u(m) =                         &
3287                                                   water_pars(ind_w_lambda_u,st)           
3288                   surf_lsm_v(l)%c_surface(m)     = 0.0_wp
3289                   surf_lsm_v(l)%albedo_type(ind_wat,m) =                      &
3290                                                  INT( water_pars(ind_w_at,st) )
3291                   surf_lsm_v(l)%emissivity(ind_wat,m)  =                      &
3292                                                  water_pars(ind_w_emis,st)
3293                ENDIF
3294             ENDDO
3295          ENDDO
3296       ENDIF     
3297
3298!
3299!--    Level 3, initialization of water parameters at single (x,y)
3300!--    position via water_pars read from file.
3301       IF ( water_pars_f%from_file )  THEN
3302!
3303!--       Horizontal surfaces
3304          DO  m = 1, surf_lsm_h%ns
3305             i = surf_lsm_h%i(m)
3306             j = surf_lsm_h%j(m)
3307!
3308!--          If surface element is not a water surface and any value in
3309!--          water_pars is given, neglect this information and give an
3310!--          informative message that this value will not be used.   
3311             IF ( .NOT. surf_lsm_h%water_surface(m)  .AND.                     &
3312                   ANY( water_pars_f%pars_xy(:,j,i) /= water_pars_f%fill ) )  THEN
3313                WRITE( message_string, * )                                     &
3314                              'surface element at grid point (j,i) = (',       &
3315                              j, i, ') is not a water surface, ',              &
3316                              'so that information given in ',                 &
3317                              'water_pars at this point is neglected.' 
3318                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
3319             ELSE
3320                IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                   &
3321                     water_pars_f%fill  .AND.                                  &
3322                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
3323                      t_soil_h%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
3324
3325                IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /= water_pars_f%fill ) &
3326                   surf_lsm_h%z0(m)     = water_pars_f%pars_xy(ind_w_z0,j,i)
3327
3328                IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /= water_pars_f%fill )&
3329                THEN
3330                   surf_lsm_h%z0h(m)    = water_pars_f%pars_xy(ind_w_z0h,j,i)
3331                   surf_lsm_h%z0q(m)    = water_pars_f%pars_xy(ind_w_z0h,j,i)
3332                ENDIF
3333
3334                IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /=               &
3335                     water_pars_f%fill )                                       &
3336                   surf_lsm_h%lambda_surface_s(m) =                            &
3337                                        water_pars_f%pars_xy(ind_w_lambda_s,j,i)
3338
3339                IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /=               &
3340                      water_pars_f%fill )                                      &
3341                   surf_lsm_h%lambda_surface_u(m) =                            &
3342                                        water_pars_f%pars_xy(ind_w_lambda_u,j,i)     
3343       
3344                IF ( water_pars_f%pars_xy(ind_w_at,j,i) /=                     &
3345                     water_pars_f%fill )                                       &
3346                   surf_lsm_h%albedo_type(ind_wat,m) =                         &
3347                                       INT( water_pars_f%pars_xy(ind_w_at,j,i) )
3348
3349                IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /=                   &
3350                     water_pars_f%fill )                                       &
3351                   surf_lsm_h%emissivity(ind_wat,m) =                          &
3352                   water_pars_f%pars_xy(ind_w_emis,j,i) 
3353             ENDIF
3354          ENDDO
3355!
3356!--       Vertical surfaces
3357          DO  l = 0, 3
3358             DO  m = 1, surf_lsm_v(l)%ns
3359                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3360                                                surf_lsm_v(l)%building_covered(m) ) 
3361                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3362                                                surf_lsm_v(l)%building_covered(m) ) 
3363!
3364!--             If surface element is not a water surface and any value in
3365!--             water_pars is given, neglect this information and give an
3366!--             informative message that this value will not be used.   
3367                IF ( .NOT. surf_lsm_v(l)%water_surface(m)  .AND.               &
3368                      ANY( water_pars_f%pars_xy(:,j,i) /=                      &
3369                      water_pars_f%fill ) )  THEN
3370                   WRITE( message_string, * )                                  &
3371                              'surface element at grid point (j,i) = (',       &
3372                              j, i, ') is not a water surface, ',              &
3373                              'so that information given in ',                 &
3374                              'water_pars at this point is neglected.' 
3375                   CALL message( 'land_surface_model_mod', 'PA0999',           &
3376                                  0, 0, 0, 6, 0 )
3377                ELSE
3378
3379                   IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                &
3380                     water_pars_f%fill  .AND.                                  &
3381                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
3382                      t_soil_v(l)%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
3383
3384                   IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /=                  &
3385                        water_pars_f%fill )                                    &
3386                      surf_lsm_v(l)%z0(m)   = water_pars_f%pars_xy(ind_w_z0,j,i)
3387
3388                   IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /=                 &
3389                       water_pars_f%fill )  THEN
3390                      surf_lsm_v(l)%z0h(m)  = water_pars_f%pars_xy(ind_w_z0h,j,i)
3391                      surf_lsm_v(l)%z0q(m)  = water_pars_f%pars_xy(ind_w_z0h,j,i)
3392                   ENDIF
3393
3394                   IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /=            &
3395                        water_pars_f%fill )                                    &
3396                      surf_lsm_v(l)%lambda_surface_s(m) =                      &
3397                                      water_pars_f%pars_xy(ind_w_lambda_s,j,i)
3398
3399                   IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /=            &
3400                        water_pars_f%fill )                                    &
3401                      surf_lsm_v(l)%lambda_surface_u(m) =                      &
3402                                      water_pars_f%pars_xy(ind_w_lambda_u,j,i)   
3403 
3404                   IF ( water_pars_f%pars_xy(ind_w_at,j,i) /=                  &
3405                        water_pars_f%fill )                                    &
3406                      surf_lsm_v(l)%albedo_type(ind_wat,m)    =                &
3407                                      INT( water_pars_f%pars_xy(ind_w_at,j,i) )
3408
3409                   IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /=                &
3410                        water_pars_f%fill )                                    &
3411                      surf_lsm_v(l)%emissivity(ind_wat,m)     =                &
3412                                      water_pars_f%pars_xy(ind_w_emis,j,i) 
3413                ENDIF
3414             ENDDO
3415          ENDDO
3416
3417       ENDIF
3418!
3419!--    Initialize pavement-type surfaces, level 1
3420       IF ( pavement_type /= 0 )  THEN 
3421
3422!
3423!--       When a pavement_type is used, overwrite a possible setting of
3424!--       the pavement depth as it is already defined by the pavement type
3425          pavement_depth_level = 0
3426
3427          IF ( z0_pavement == 9999999.9_wp )  THEN
3428             z0_pavement  = pavement_pars(ind_p_z0,pavement_type) 
3429          ENDIF
3430
3431          IF ( z0h_pavement == 9999999.9_wp )  THEN
3432             z0h_pavement = pavement_pars(ind_p_z0h,pavement_type)
3433          ENDIF
3434
3435          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
3436             pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type)
3437          ENDIF
3438
3439          IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
3440             pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)
3441          ENDIF   
3442   
3443          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
3444             albedo_type = INT(pavement_pars(ind_p_at,pavement_type))       
3445          ENDIF
3446   
3447          IF ( emissivity == 9999999.9_wp )  THEN
3448             emissivity = pavement_pars(ind_p_emis,pavement_type)       
3449          ENDIF
3450
3451!
3452!--       If the depth level of the pavement is not set, determine it from
3453!--       lookup table.
3454          IF ( pavement_depth_level == 0 )  THEN
3455             DO  k = nzb_soil, nzt_soil 
3456                IF ( pavement_subsurface_pars_1(k,pavement_type) == 9999999.9_wp &
3457                .OR. pavement_subsurface_pars_2(k,pavement_type) == 9999999.9_wp)&
3458                THEN
3459                   nzt_pavement = k-1
3460                   EXIT
3461                ENDIF
3462             ENDDO
3463          ELSE
3464             nzt_pavement = pavement_depth_level
3465          ENDIF
3466
3467       ENDIF
3468!
3469!--    Level 1 initialization of pavement type surfaces. Horizontally
3470!--    homogeneous characteristics are assumed
3471       surf_lsm_h%nzt_pavement = pavement_depth_level
3472       DO  m = 1, surf_lsm_h%ns
3473          IF ( surf_lsm_h%pavement_surface(m) )  THEN
3474             surf_lsm_h%nzt_pavement(m)        = nzt_pavement
3475             surf_lsm_h%z0(m)                  = z0_pavement
3476             surf_lsm_h%z0h(m)                 = z0h_pavement
3477             surf_lsm_h%z0q(m)                 = z0h_pavement
3478             surf_lsm_h%lambda_surface_s(m)    = pavement_heat_conduct         &
3479                                                  * ddz_soil(nzb_soil)         &
3480                                                  * 2.0_wp   
3481             surf_lsm_h%lambda_surface_u(m)    = pavement_heat_conduct         &
3482                                                  * ddz_soil(nzb_soil)         &
3483                                                  * 2.0_wp           
3484             surf_lsm_h%c_surface(m)           = pavement_heat_capacity        &
3485                                                        * dz_soil(nzb_soil)    &
3486                                                        * 0.25_wp                                   
3487
3488             surf_lsm_h%albedo_type(ind_pav,m) = albedo_type
3489             surf_lsm_h%emissivity(ind_pav,m)  = emissivity     
3490     
3491             IF ( pavement_type /= 0 )  THEN
3492                DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
3493                   surf_lsm_h%lambda_h_def(k,m)    =                           &
3494                                     pavement_subsurface_pars_1(k,pavement_type)                       
3495                   surf_lsm_h%rho_c_total_def(k,m) =                           &
3496                                     pavement_subsurface_pars_2(k,pavement_type) 
3497                ENDDO
3498             ELSE
3499                surf_lsm_v(l)%lambda_h_def(:,m)     = pavement_heat_conduct
3500                surf_lsm_v(l)%rho_c_total_def(:,m)  = pavement_heat_capacity
3501             ENDIF       
3502          ENDIF
3503       ENDDO                               
3504
3505       DO  l = 0, 3
3506          surf_lsm_v(l)%nzt_pavement = pavement_depth_level
3507          DO  m = 1, surf_lsm_v(l)%ns
3508             IF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
3509                surf_lsm_v(l)%nzt_pavement(m)        = nzt_pavement
3510                surf_lsm_v(l)%z0(m)                  = z0_pavement
3511                surf_lsm_v(l)%z0h(m)                 = z0h_pavement
3512                surf_lsm_v(l)%z0q(m)                 = z0h_pavement
3513                surf_lsm_v(l)%lambda_surface_s(m)    = pavement_heat_conduct   &
3514                                                  * ddz_soil(nzb_soil)         &
3515                                                  * 2.0_wp   
3516                surf_lsm_v(l)%lambda_surface_u(m)    = pavement_heat_conduct   &
3517                                                  * ddz_soil(nzb_soil)         &
3518                                                  * 2.0_wp           
3519                surf_lsm_v(l)%c_surface(m)           = pavement_heat_capacity  &
3520                                                        * dz_soil(nzb_soil)    &
3521                                                        * 0.25_wp                                     
3522
3523                surf_lsm_v(l)%albedo_type(ind_pav,m) = albedo_type
3524                surf_lsm_v(l)%emissivity(ind_pav,m)  = emissivity
3525
3526                IF ( pavement_type /= 0 )  THEN
3527                   DO  k = nzb_soil, surf_lsm_v(l)%nzt_pavement(m)
3528                      surf_lsm_v(l)%lambda_h_def(k,m)    =                     &
3529                                     pavement_subsurface_pars_1(k,pavement_type)                       
3530                      surf_lsm_v(l)%rho_c_total_def(k,m) =                     &
3531                                     pavement_subsurface_pars_2(k,pavement_type) 
3532                   ENDDO
3533                ELSE
3534                   surf_lsm_v(l)%lambda_h_def(:,m)     = pavement_heat_conduct
3535                   surf_lsm_v(l)%rho_c_total_def(:,m)  = pavement_heat_capacity
3536                ENDIF     
3537             ENDIF
3538          ENDDO
3539       ENDDO
3540!
3541!--    Level 2 initialization of pavement type surfaces via pavement_type read
3542!--    from file. Pavement surfaces are initialized for each (y,x)-grid point
3543!--    individually.
3544       IF ( pavement_type_f%from_file )  THEN
3545!
3546!--       Horizontal surfaces
3547          DO  m = 1, surf_lsm_h%ns
3548             i = surf_lsm_h%i(m)
3549             j = surf_lsm_h%j(m)
3550             
3551             st = pavement_type_f%var(j,i)
3552             IF ( st /= pavement_type_f%fill  .AND.  st /= 0 )  THEN
3553!
3554!--             Determine deepmost index of pavement layer
3555                DO  k = nzb_soil, nzt_soil 
3556                   IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp       &
3557                   .OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp)      &
3558                   THEN
3559                      surf_lsm_h%nzt_pavement(m) = k-1
3560                      EXIT
3561                   ENDIF
3562                ENDDO
3563
3564                surf_lsm_h%z0(m)                = pavement_pars(ind_p_z0,st)
3565                surf_lsm_h%z0h(m)               = pavement_pars(ind_p_z0h,st)
3566                surf_lsm_h%z0q(m)               = pavement_pars(ind_p_z0h,st)
3567
3568                surf_lsm_h%lambda_surface_s(m)  =                              &
3569                                              pavement_subsurface_pars_1(0,st) &
3570                                                  * ddz_soil(nzb_soil)         &
3571                                                  * 2.0_wp   
3572                surf_lsm_h%lambda_surface_u(m)  =                              &
3573                                              pavement_subsurface_pars_1(0,st) &
3574                                                  * ddz_soil(nzb_soil)         &
3575                                                  * 2.0_wp       
3576                surf_lsm_h%c_surface(m)         =                              &
3577                                               pavement_subsurface_pars_2(0,st)&
3578                                                        * dz_soil(nzb_soil)    &
3579                                                        * 0.25_wp                               
3580                surf_lsm_h%albedo_type(ind_pav,m) = INT( pavement_pars(ind_p_at,st) )
3581                surf_lsm_h%emissivity(ind_pav,m)  = pavement_pars(ind_p_emis,st) 
3582
3583                DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
3584                   surf_lsm_h%lambda_h_def(k,m)    =                           &
3585                                     pavement_subsurface_pars_1(k,pavement_type)                       
3586                   surf_lsm_h%rho_c_total_def(k,m) =                           &
3587                                     pavement_subsurface_pars_2(k,pavement_type) 
3588                ENDDO   
3589             ENDIF
3590          ENDDO
3591!
3592!--       Vertical surfaces
3593          DO  l = 0, 3
3594             DO  m = 1, surf_lsm_v(l)%ns
3595                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3596                                                surf_lsm_v(l)%building_covered(m) ) 
3597                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3598                                                surf_lsm_v(l)%building_covered(m) ) 
3599             
3600                st = pavement_type_f%var(j,i)
3601                IF ( st /= pavement_type_f%fill  .AND.  st /= 0 )  THEN
3602!
3603!--                Determine deepmost index of pavement layer
3604                   DO  k = nzb_soil, nzt_soil 
3605                      IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp    &
3606                      .OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp)   &
3607                      THEN
3608                         surf_lsm_v(l)%nzt_pavement(m) = k-1
3609                         EXIT
3610                      ENDIF
3611                   ENDDO
3612
3613                   surf_lsm_v(l)%z0(m)  = pavement_pars(ind_p_z0,st)
3614                   surf_lsm_v(l)%z0h(m) = pavement_pars(ind_p_z0h,st)
3615                   surf_lsm_v(l)%z0q(m) = pavement_pars(ind_p_z0h,st)
3616
3617                   surf_lsm_v(l)%lambda_surface_s(m)  =                        &
3618                                              pavement_subsurface_pars_1(0,st) &
3619                                                  * ddz_soil(nzb_soil)         & 
3620                                                  * 2.0_wp   
3621                   surf_lsm_v(l)%lambda_surface_u(m)  =                        &
3622                                              pavement_subsurface_pars_1(0,st) &
3623                                                  * ddz_soil(nzb_soil)         &
3624                                                  * 2.0_wp     
3625
3626                   surf_lsm_v(l)%c_surface(m)    =                             &
3627                                           pavement_subsurface_pars_2(0,st)    &
3628                                                        * dz_soil(nzb_soil)    &
3629                                                        * 0.25_wp                                   
3630                   surf_lsm_v(l)%albedo_type(ind_pav,m) =                      &
3631                                              INT( pavement_pars(ind_p_at,st) )
3632                   surf_lsm_v(l)%emissivity(ind_pav,m)  =                      &
3633                                              pavement_pars(ind_p_emis,st)   
3634
3635                   DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
3636                      surf_lsm_v(l)%lambda_h_def(k,m)    =                    &
3637                                    pavement_subsurface_pars_1(k,pavement_type)                       
3638                      surf_lsm_v(l)%rho_c_total_def(k,m) =                    &
3639                                    pavement_subsurface_pars_2(k,pavement_type) 
3640                   ENDDO   
3641                ENDIF
3642             ENDDO
3643          ENDDO
3644       ENDIF
3645!
3646!--    Level 3, initialization of pavement parameters at single (x,y)
3647!--    position via pavement_pars read from file.
3648       IF ( pavement_pars_f%from_file )  THEN
3649!
3650!--       Horizontal surfaces
3651          DO  m = 1, surf_lsm_h%ns
3652             i = surf_lsm_h%i(m)
3653             j = surf_lsm_h%j(m)
3654!
3655!--          If surface element is not a pavement surface and any value in
3656!--          pavement_pars is given, neglect this information and give an
3657!--          informative message that this value will not be used.   
3658             IF ( .NOT. surf_lsm_h%pavement_surface(m)  .AND.                  &
3659                   ANY( pavement_pars_f%pars_xy(:,j,i) /=                      &
3660                   pavement_pars_f%fill ) )  THEN
3661                WRITE( message_string, * )                                     &
3662                              'surface element at grid point (j,i) = (',       &
3663                              j, i, ') is not a pavement surface, ',           &
3664                              'so that information given in ',                 &
3665                              'pavement_pars at this point is neglected.' 
3666                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
3667             ELSE
3668                IF ( pavement_pars_f%pars_xy(ind_p_z0,j,i) /=                  &
3669                     pavement_pars_f%fill )                                    &
3670                   surf_lsm_h%z0(m)  = pavement_pars_f%pars_xy(ind_p_z0,j,i)
3671                IF ( pavement_pars_f%pars_xy(ind_p_z0h,j,i) /=                 &
3672                     pavement_pars_f%fill )  THEN
3673                   surf_lsm_h%z0h(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
3674                   surf_lsm_h%z0q(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
3675                ENDIF
3676                IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i) &
3677                     /= pavement_subsurface_pars_f%fill )  THEN
3678                   surf_lsm_h%lambda_surface_s(m)  =                           &
3679                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
3680                                                  * ddz_soil(nzb_soil)         &
3681                                                  * 2.0_wp
3682                   surf_lsm_h%lambda_surface_u(m)  =                           &
3683                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
3684                                                  * ddz_soil(nzb_soil)         &
3685                                                  * 2.0_wp   
3686                ENDIF
3687                IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) /= &
3688                     pavement_subsurface_pars_f%fill )  THEN
3689                   surf_lsm_h%c_surface(m)     =                               &
3690                      pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i)   &
3691                                                  * dz_soil(nzb_soil)          &
3692                                                  * 0.25_wp                                   
3693                ENDIF
3694                IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /=                  &
3695                     pavement_pars_f%fill )                                    &
3696                   surf_lsm_h%albedo_type(ind_pav,m) =                         &
3697                                              INT( pavement_pars(ind_p_at,st) )
3698                IF ( pavement_pars_f%pars_xy(ind_p_emis,j,i) /=                &
3699                     pavement_pars_f%fill )                                    &
3700                   surf_lsm_h%emissivity(ind_pav,m)  =                         &
3701                                              pavement_pars(ind_p_emis,st) 
3702             ENDIF
3703
3704          ENDDO
3705!
3706!--       Vertical surfaces
3707          DO  l = 0, 3
3708             DO  m = 1, surf_lsm_v(l)%ns
3709                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3710                                                surf_lsm_v(l)%building_covered(m) ) 
3711                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3712                                                surf_lsm_v(l)%building_covered(m) ) 
3713!
3714!--             If surface element is not a pavement surface and any value in
3715!--             pavement_pars is given, neglect this information and give an
3716!--             informative message that this value will not be used.   
3717                IF ( .NOT. surf_lsm_v(l)%pavement_surface(m)  .AND.            &
3718                      ANY( pavement_pars_f%pars_xy(:,j,i) /=                   &
3719                      pavement_pars_f%fill ) )  THEN
3720                   WRITE( message_string, * )                                  &
3721                                 'surface element at grid point (j,i) = (',    &
3722                                 j, i, ') is not a pavement surface, ',        &
3723                                 'so that information given in ',              &
3724                                 'pavement_pars at this point is neglected.' 
3725                   CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
3726                ELSE
3727
3728                   IF ( pavement_pars_f%pars_xy(ind_p_z0,j,i) /=               &
3729                        pavement_pars_f%fill )                                 &
3730                      surf_lsm_v(l)%z0(m) = pavement_pars_f%pars_xy(ind_p_z0,j,i)
3731                   IF ( pavement_pars_f%pars_xy(ind_p_z0h,j,i) /=              &
3732                        pavement_pars_f%fill )  THEN
3733                      surf_lsm_v(l)%z0h(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
3734                      surf_lsm_v(l)%z0q(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
3735                   ENDIF
3736                   IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
3737                        /= pavement_subsurface_pars_f%fill )  THEN
3738                      surf_lsm_v(l)%lambda_surface_s(m) =                      &
3739                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
3740                                                  * ddz_soil(nzb_soil)         &
3741                                                  * 2.0_wp
3742                      surf_lsm_v(l)%lambda_surface_u(m) =                      &
3743                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
3744                                                  * ddz_soil(nzb_soil)         &
3745                                                  * 2.0_wp   
3746                   ENDIF
3747                   IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) &
3748                        /= pavement_subsurface_pars_f%fill )  THEN
3749                      surf_lsm_v(l)%c_surface(m)    =                          &
3750                         pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i)&
3751                                                  * dz_soil(nzb_soil)          &
3752                                                  * 0.25_wp                                 
3753                   ENDIF
3754                   IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /=               &
3755                        pavement_pars_f%fill )                                 &
3756                      surf_lsm_v(l)%albedo_type(ind_pav,m) =                   &
3757                                            INT( pavement_pars(ind_p_at,st) )
3758
3759                   IF ( pavement_pars_f%pars_xy(ind_p_emis,j,i) /=             &
3760                        pavement_pars_f%fill )                                 &
3761                      surf_lsm_v(l)%emissivity(ind_pav,m)  =                   &
3762                                            pavement_pars(ind_p_emis,st) 
3763                ENDIF
3764             ENDDO
3765          ENDDO
3766       ENDIF
3767!
3768!--    Moreover, for grid points which are flagged with pavement-type 0 or whre
3769!--    pavement_subsurface_pars_f is provided, soil heat conductivity and
3770!--    capacity are initialized with parameters given in       
3771!--    pavement_subsurface_pars read from file.
3772       IF ( pavement_subsurface_pars_f%from_file )  THEN
3773!
3774!--       Set pavement depth to nzt_soil. Please note, this is just a
3775!--       workaround at the moment.
3776          DO  m = 1, surf_lsm_h%ns
3777             IF ( surf_lsm_h%pavement_surface(m) )  THEN
3778
3779                i = surf_lsm_h%i(m)
3780                j = surf_lsm_h%j(m)
3781
3782                surf_lsm_h%nzt_pavement(m) = nzt_soil
3783
3784                DO  k = nzb_soil, nzt_soil 
3785                   surf_lsm_h%lambda_h_def(k,m) =                              &
3786                       pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,k,j,i)
3787                   surf_lsm_h%rho_c_total_def(k,m) =                           &
3788                       pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,k,j,i)
3789                ENDDO
3790
3791             ENDIF
3792          ENDDO
3793          DO  l = 0, 3
3794             DO  m = 1, surf_lsm_v(l)%ns
3795                IF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
3796
3797                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
3798                                                surf_lsm_v(l)%building_covered(m) ) 
3799                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
3800                                                surf_lsm_v(l)%building_covered(m) ) 
3801
3802                   surf_lsm_v(l)%nzt_pavement(m) = nzt_soil
3803
3804                   DO  k = nzb_soil, nzt_soil 
3805                      surf_lsm_v(l)%lambda_h_def(k,m) =                        &
3806                       pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,k,j,i)
3807                      surf_lsm_v(l)%rho_c_total_def(k,m) =                     &
3808                       pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,k,j,i)
3809                   ENDDO
3810
3811                ENDIF
3812             ENDDO
3813          ENDDO
3814       ENDIF
3815!
3816!--    Initial run actions
3817       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
3818!
3819!--       First, initialize soil temperature and moisture.
3820!--       According to the initialization for surface and soil parameters,
3821!--       initialize soil moisture and temperature via a level approach. This
3822!--       is to assure that all surface elements are initialized, even if
3823!--       data provided from input file contains fill values at some locations.
3824!--       Level 1, initialization via profiles given in parameter file
3825          DO  m = 1, surf_lsm_h%ns
3826             IF ( surf_lsm_h%vegetation_surface(m)  .OR.                       &
3827                  surf_lsm_h%pavement_surface(m) )  THEN
3828                DO  k = nzb_soil, nzt_soil 
3829                   t_soil_h%var_2d(k,m) = soil_temperature(k)
3830                   m_soil_h%var_2d(k,m) = soil_moisture(k)
3831                ENDDO
3832                t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
3833             ENDIF
3834          ENDDO
3835          DO  l = 0, 3
3836             DO  m = 1, surf_lsm_v(l)%ns
3837                IF ( surf_lsm_v(l)%vegetation_surface(m)  .OR.                 &
3838                     surf_lsm_v(l)%pavement_surface(m) )  THEN
3839                   DO  k = nzb_soil, nzt_soil 
3840                      t_soil_v(l)%var_2d(k,m) = soil_temperature(k)
3841                      m_soil_v(l)%var_2d(k,m) = soil_moisture(k)
3842                   ENDDO
3843                   t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
3844                ENDIF
3845             ENDDO
3846          ENDDO
3847!
3848!--       Level 2, if soil moisture and/or temperature  are
3849!--       provided from file, interpolate / extrapolate the provided data
3850!--       onto the respective soil layers. Please note, both, zs as well as
3851!--       init_3d%z_soil indicate a depth with positive values, so that no
3852!--       distinction between atmosphere is required concerning interpolation.
3853!--       Start with soil moisture
3854          IF ( init_3d%from_file_msoil )  THEN
3855
3856             IF ( init_3d%lod_msoil == 1 )  THEN
3857                DO  m = 1, surf_lsm_h%ns
3858
3859                   CALL netcdf_data_input_interpolate(                         &
3860                                       m_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
3861                                       init_3d%msoil_init(:),                  &
3862                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
3863                                       nzb_soil, nzt_soil,                     &
3864                                       nzb_soil, init_3d%nzs-1 )
3865                ENDDO
3866                DO  l = 0, 3
3867                   DO  m = 1, surf_lsm_v(l)%ns
3868
3869                      CALL netcdf_data_input_interpolate(                      &
3870                                       m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
3871                                       init_3d%msoil_init(:),                  &
3872                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
3873                                       nzb_soil, nzt_soil,                     &
3874                                       nzb_soil, init_3d%nzs-1 )
3875                   ENDDO
3876                ENDDO
3877             ELSE
3878
3879                DO  m = 1, surf_lsm_h%ns
3880                   i = surf_lsm_h%i(m)
3881                   j = surf_lsm_h%j(m)
3882
3883                   IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )           &
3884                      CALL netcdf_data_input_interpolate(                      &
3885                                       m_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
3886                                       init_3d%msoil(:,j,i),                   &
3887                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
3888                                       nzb_soil, nzt_soil,                     &
3889                                       nzb_soil, init_3d%nzs-1 )
3890                ENDDO
3891                DO  l = 0, 3
3892                   DO  m = 1, surf_lsm_v(l)%ns
3893                      i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,   &
3894                                             surf_lsm_v(l)%building_covered(m) ) 
3895                      j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,   &
3896                                             surf_lsm_v(l)%building_covered(m) ) 
3897
3898                      IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )        &
3899                         CALL netcdf_data_input_interpolate(                   &
3900                                       m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
3901                                       init_3d%msoil(:,j,i),                   &
3902                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
3903                                       nzb_soil, nzt_soil,                     &
3904                                       nzb_soil, init_3d%nzs-1 )
3905                   ENDDO
3906                ENDDO
3907             ENDIF
3908
3909          ENDIF
3910!
3911!--       Soil temperature
3912          IF ( init_3d%from_file_tsoil )  THEN
3913
3914             IF ( init_3d%lod_tsoil == 1 )  THEN ! change to 1 if provided correctly by INIFOR
3915                DO  m = 1, surf_lsm_h%ns
3916
3917                   CALL netcdf_data_input_interpolate(                         &
3918                                       t_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
3919                                       init_3d%tsoil_init(:),                  &
3920                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
3921                                       nzb_soil, nzt_soil,                     &
3922                                       nzb_soil, init_3d%nzs-1 )
3923                   t_soil_h%var_2d(nzt_soil+1,m) = t_soil_h%var_2d(nzt_soil,m)
3924                ENDDO
3925                DO  l = 0, 3
3926                   DO  m = 1, surf_lsm_v(l)%ns
3927
3928                      CALL netcdf_data_input_interpolate(                      &
3929                                       t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
3930                                       init_3d%tsoil_init(:),                  &
3931                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
3932                                       nzb_soil, nzt_soil,                     &
3933                                       nzb_soil, init_3d%nzs-1 )
3934                      t_soil_v(l)%var_2d(nzt_soil+1,m) =                       &
3935                                                 t_soil_v(l)%var_2d(nzt_soil,m)
3936                   ENDDO
3937                ENDDO
3938             ELSE
3939
3940                DO  m = 1, surf_lsm_h%ns
3941                   i = surf_lsm_h%i(m)
3942                   j = surf_lsm_h%j(m)
3943
3944                   IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )           &
3945                      CALL netcdf_data_input_interpolate(                      &
3946                                       t_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
3947                                       init_3d%tsoil(:,j,i),                   &
3948                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
3949                                       nzb_soil, nzt_soil,                     &
3950                                       nzb_soil, init_3d%nzs-1 )
3951                   t_soil_h%var_2d(nzt_soil+1,m) = t_soil_h%var_2d(nzt_soil,m)
3952                ENDDO
3953                DO  l = 0, 3
3954                   DO  m = 1, surf_lsm_v(l)%ns
3955                      i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,   &
3956                                                surf_lsm_v(l)%building_covered(m) ) 
3957                      j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,   &
3958                                                surf_lsm_v(l)%building_covered(m) ) 
3959
3960                      IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )        &
3961                         CALL netcdf_data_input_interpolate(                   &
3962                                       t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
3963                                       init_3d%tsoil(:,j,i),                   &
3964                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
3965                                       nzb_soil, nzt_soil,                     &
3966                                       nzb_soil, init_3d%nzs-1 )
3967                      t_soil_v(l)%var_2d(nzt_soil+1,m) =                       &
3968                                                 t_soil_v(l)%var_2d(nzt_soil,m)
3969                   ENDDO
3970                ENDDO
3971             ENDIF
3972          ENDIF
3973!
3974!--       Further initialization
3975          DO  m = 1, surf_lsm_h%ns
3976
3977             i   = surf_lsm_h%i(m)           
3978             j   = surf_lsm_h%j(m)
3979             k   = surf_lsm_h%k(m)
3980!
3981!--          Calculate surface temperature. In case of bare soil, the surface
3982!--          temperature must be reset to the soil temperature in the first soil
3983!--          layer
3984             IF ( surf_lsm_h%lambda_surface_s(m) == 0.0_wp )  THEN
3985                t_surface_h%var_1d(:)    = t_soil_h%var_2d(nzb_soil,m)
3986                surf_lsm_h%pt_surface(m) = t_soil_h%var_2d(nzb_soil,m) / exn
3987             ELSE
3988                t_surface_h%var_1d(m)    = pt(k-1,j,i) * exn
3989                surf_lsm_h%pt_surface(m) = pt(k-1,j,i) 
3990             ENDIF
3991             
3992             IF ( cloud_physics  .OR. cloud_droplets ) THEN
3993                surf_lsm_h%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
3994             ELSE
3995                surf_lsm_h%pt1(m) = pt(k,j,i)
3996             ENDIF
3997
3998
3999!
4000!--          Assure that r_a cannot be zero at model start
4001             IF ( surf_lsm_h%pt1(m) == surf_lsm_h%pt_surface(m) )              &
4002                surf_lsm_h%pt1(m) = surf_lsm_h%pt1(m) + 1.0E-20_wp
4003
4004             surf_lsm_h%us(m)   = 0.1_wp
4005             surf_lsm_h%ts(m)   = ( surf_lsm_h%pt1(m) - surf_lsm_h%pt_surface(m) )&
4006                                  / surf_lsm_h%r_a(m)
4007             surf_lsm_h%shf(m)  = - surf_lsm_h%us(m) * surf_lsm_h%ts(m)        &
4008                                  * rho_surface
4009         ENDDO
4010!
4011!--      Vertical surfaces
4012         DO  l = 0, 3
4013             DO  m = 1, surf_lsm_v(l)%ns
4014                i   = surf_lsm_v(l)%i(m)           
4015                j   = surf_lsm_v(l)%j(m)
4016                k   = surf_lsm_v(l)%k(m)         
4017!
4018!--             Calculate surface temperature. In case of bare soil, the surface
4019!--             temperature must be reset to the soil temperature in the first soil
4020!--             layer
4021                IF ( surf_lsm_v(l)%lambda_surface_s(m) == 0.0_wp )  THEN
4022                   t_surface_v(l)%var_1d(m)      = t_soil_v(l)%var_2d(nzb_soil,m)
4023                   surf_lsm_v(l)%pt_surface(m)   = t_soil_v(l)%var_2d(nzb_soil,m) / exn
4024                ELSE
4025                   j_off = surf_lsm_v(l)%joff
4026                   i_off = surf_lsm_v(l)%ioff
4027
4028                   t_surface_v(l)%var_1d(m)      = pt(k,j+j_off,i+i_off) * exn
4029                   surf_lsm_v(l)%pt_surface(m)   = pt(k,j+j_off,i+i_off)           
4030                ENDIF
4031
4032
4033                IF ( cloud_physics  .OR. cloud_droplets ) THEN
4034                   surf_lsm_v(l)%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
4035                ELSE
4036                   surf_lsm_v(l)%pt1(m) = pt(k,j,i)
4037                ENDIF
4038
4039!
4040!--             Assure that r_a cannot be zero at model start
4041                IF ( surf_lsm_v(l)%pt1(m) == surf_lsm_v(l)%pt_surface(m) )     &
4042                     surf_lsm_v(l)%pt1(m) = surf_lsm_v(l)%pt1(m) + 1.0E-20_wp
4043!
4044!--             Set artifical values for ts and us so that r_a has its initial value
4045!--             for the first time step. Only for interior core domain, not for ghost points
4046                surf_lsm_v(l)%us(m)   = 0.1_wp
4047                surf_lsm_v(l)%ts(m)   = ( surf_lsm_v(l)%pt1(m) - surf_lsm_v(l)%pt_surface(m) ) /&
4048                                          surf_lsm_v(l)%r_a(m)
4049                surf_lsm_v(l)%shf(m)  = - surf_lsm_v(l)%us(m) *                &
4050                                          surf_lsm_v(l)%ts(m) * rho_surface
4051
4052             ENDDO
4053          ENDDO
4054
4055!
4056!--    Actions for restart runs
4057       ELSE
4058!
4059!--       Horizontal surfaces
4060          DO  m = 1, surf_lsm_h%ns
4061             i   = surf_lsm_h%i(m)           
4062             j   = surf_lsm_h%j(m)
4063             k   = surf_lsm_h%k(m)         
4064             t_surface_h%var_1d(m) = pt(k-1,j,i) * exn
4065          ENDDO
4066!
4067!--       Vertical surfaces
4068          DO  l = 0, 3
4069!
4070!--          Set index offset of surface element, seen from atmospheric grid point
4071             j_off = surf_lsm_v(l)%joff
4072             i_off = surf_lsm_v(l)%ioff
4073
4074             DO  m = 1, surf_lsm_v(l)%ns
4075                i   = surf_lsm_v(l)%i(m)           
4076                j   = surf_lsm_v(l)%j(m)
4077                k   = surf_lsm_v(l)%k(m)         
4078                t_surface_v(l)%var_1d(m) = pt(k,j+j_off,i+i_off) * exn
4079             ENDDO
4080          ENDDO
4081
4082       ENDIF
4083!
4084!--    Level 1 initialization of root distribution - provided by the user via
4085!--    via namelist.
4086       DO  m = 1, surf_lsm_h%ns
4087          DO  k = nzb_soil, nzt_soil
4088             surf_lsm_h%root_fr(k,m) = root_fraction(k)
4089          ENDDO
4090       ENDDO
4091
4092       DO  l = 0, 3
4093          DO  m = 1, surf_lsm_v(l)%ns
4094             DO  k = nzb_soil, nzt_soil
4095                surf_lsm_v(l)%root_fr(k,m) = root_fraction(k)
4096             ENDDO
4097          ENDDO
4098       ENDDO
4099
4100!
4101!--    Level 2 initialization of root distribution.
4102!--    When no root distribution is given by the user, use look-up table to prescribe
4103!--    the root fraction in the individual soil layers.
4104       IF ( ALL( root_fraction == 9999999.9_wp ) )  THEN
4105!
4106!--       First, calculate the index bounds for integration
4107          n_soil_layers_total = nzt_soil - nzb_soil + 6
4108          ALLOCATE ( bound(0:n_soil_layers_total) )
4109          ALLOCATE ( bound_root_fr(0:n_soil_layers_total) )
4110
4111          kn = 0
4112          ko = 0
4113          bound(0) = 0.0_wp
4114          DO k = 1, n_soil_layers_total-1
4115             IF ( zs_layer(kn) <= zs_ref(ko) )  THEN
4116                bound(k) = zs_layer(kn)
4117                bound_root_fr(k) = ko
4118                kn = kn + 1
4119                IF ( kn > nzt_soil+1 )  THEN
4120                   kn = nzt_soil
4121                ENDIF
4122             ELSE
4123                bound(k) = zs_ref(ko)
4124                bound_root_fr(k) = ko
4125                ko = ko + 1
4126                IF ( ko > 3 )  THEN
4127                   ko = 3
4128                ENDIF
4129             ENDIF
4130
4131          ENDDO
4132
4133!
4134!--       Integrate over all soil layers based on the four-layer root fraction
4135          kzs = 1
4136          root_fraction = 0.0_wp
4137          DO k = 0, n_soil_layers_total-2
4138             kroot = bound_root_fr(k+1)
4139             root_fraction(kzs-1) = root_fraction(kzs-1)                       &
4140                                + root_distribution(kroot,vegetation_type)     &
4141                                / dz_soil_ref(kroot) * ( bound(k+1) - bound(k) )
4142
4143             IF ( bound(k+1) == zs_layer(kzs-1) )  THEN
4144                kzs = kzs+1
4145             ENDIF
4146          ENDDO
4147
4148
4149!
4150!--       Normalize so that the sum of all fractions equals one
4151          root_fraction = root_fraction / SUM(root_fraction)
4152
4153          DEALLOCATE ( bound )
4154          DEALLOCATE ( bound_root_fr )
4155
4156!
4157!--       Map calculated root fractions
4158          DO  m = 1, surf_lsm_h%ns
4159             DO  k = nzb_soil, nzt_soil
4160                IF ( surf_lsm_h%pavement_surface(m)  .AND.                     &
4161                     k <= surf_lsm_h%nzt_pavement(m) )  THEN
4162                   surf_lsm_h%root_fr(k,m) = 0.0_wp
4163                ELSE
4164                   surf_lsm_h%root_fr(k,m) = root_fraction(k)
4165                ENDIF
4166
4167             ENDDO
4168!
4169!--          Normalize so that the sum = 1. Only relevant when the root         
4170!--          distribution was set to zero due to pavement at some layers.
4171             IF ( SUM( surf_lsm_h%root_fr(:,m) ) > 0.0_wp )  THEN
4172                DO k = nzb_soil, nzt_soil
4173                   surf_lsm_h%root_fr(k,m) = surf_lsm_h%root_fr(k,m)           &
4174                   / SUM( surf_lsm_h%root_fr(:,m) )
4175                ENDDO
4176             ENDIF
4177          ENDDO
4178          DO  l = 0, 3
4179             DO  m = 1, surf_lsm_v(l)%ns
4180                DO  k = nzb_soil, nzt_soil
4181                   IF ( surf_lsm_v(l)%pavement_surface(m)  .AND.               &
4182                        k <= surf_lsm_h%nzt_pavement(m) )  THEN
4183                      surf_lsm_v(l)%root_fr(k,m) = 0.0_wp
4184                   ELSE
4185                      surf_lsm_v(l)%root_fr(k,m) = root_fraction(k)
4186                   ENDIF
4187                ENDDO
4188!
4189!--             Normalize so that the sum = 1. Only relevant when the root     
4190!--             distribution was set to zero due to pavement at some layers.
4191                IF ( SUM( surf_lsm_v(l)%root_fr(:,m) ) > 0.0_wp )  THEN
4192                   DO  k = nzb_soil, nzt_soil 
4193                      surf_lsm_v(l)%root_fr(k,m) = surf_lsm_v(l)%root_fr(k,m)  &
4194                      / SUM( surf_lsm_v(l)%root_fr(:,m) )
4195                   ENDDO
4196                ENDIF
4197             ENDDO
4198           ENDDO
4199       ENDIF
4200!
4201!--    Level 3 initialization of root distribution.
4202!--    Take value from file
4203       IF ( root_area_density_lsm_f%from_file )  THEN
4204          DO  m = 1, surf_lsm_h%ns
4205             IF ( surf_lsm_h%vegetation_surface(m) )  THEN
4206                i = surf_lsm_h%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,            &
4207                                             surf_lsm_v(l)%building_covered(m) ) 
4208                j = surf_lsm_h%j(m) + MERGE( 0, surf_lsm_v(l)%joff,            &
4209                                             surf_lsm_v(l)%building_covered(m) ) 
4210                DO  k = nzb_soil, nzt_soil
4211                   surf_lsm_h%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i) 
4212                ENDDO
4213
4214             ENDIF
4215          ENDDO
4216
4217          DO  l = 0, 3
4218             DO  m = 1, surf_lsm_v(l)%ns
4219                IF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
4220                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
4221                                                   surf_lsm_v(l)%building_covered(m) ) 
4222                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
4223                                                   surf_lsm_v(l)%building_covered(m) ) 
4224
4225                   DO  k = nzb_soil, nzt_soil
4226                      surf_lsm_v(l)%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i) 
4227                   ENDDO
4228
4229                ENDIF
4230             ENDDO
4231          ENDDO
4232
4233       ENDIF
4234 
4235!
4236!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
4237       CALL user_init_land_surface
4238
4239
4240!
4241!--    Calculate new roughness lengths (for water surfaces only, i.e. only
4242!-     horizontal surfaces)
4243       CALL calc_z0_water_surface
4244
4245       t_soil_h_p    = t_soil_h
4246       m_soil_h_p    = m_soil_h
4247       m_liq_h_p  = m_liq_h
4248       t_surface_h_p = t_surface_h
4249
4250       t_soil_v_p    = t_soil_v
4251       m_soil_v_p    = m_soil_v
4252       m_liq_v_p  = m_liq_v
4253       t_surface_v_p = t_surface_v
4254
4255
4256
4257!--    Store initial profiles of t_soil and m_soil (assuming they are
4258!--    horizontally homogeneous on this PE)
4259!--    DEACTIVATED FOR NOW - leads to error when number of locations with
4260!--    soil model is zero on a PE.
4261!        hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
4262!                                                 2, statistic_regions+1 )
4263!        hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
4264!                                                 2, statistic_regions+1 )
4265
4266!
4267!--    Finally, make some consistency checks.
4268!--    Ceck for eck for illegal combination of LAI and vegetation coverage.
4269       IF ( ANY( .NOT. surf_lsm_h%pavement_surface  .AND.                      &
4270                 surf_lsm_h%lai == 0.0_wp  .AND.  surf_lsm_h%c_veg == 1.0_wp ) &
4271          )  THEN
4272          message_string = 'For non-pavement surfaces the combination ' //     &
4273                           ' lai = 0.0 and c_veg = 1.0 is not allowed.'
4274          CALL message( 'lsm_read_restart_data', 'PA0999', 2, 2, 0, 6, 0 )
4275       ENDIF
4276
4277       DO  l = 0, 3
4278          IF ( ANY( .NOT. surf_lsm_v(l)%pavement_surface  .AND.                &
4279                    surf_lsm_v(l)%lai == 0.0_wp  .AND.                         &
4280                    surf_lsm_v(l)%c_veg == 1.0_wp ) )  THEN
4281             message_string = 'For non-pavement surfaces the combination ' //  &
4282                              ' lai = 0.0 and c_veg = 1.0 is not allowed.'
4283             CALL message( 'lsm_read_restart_data', 'PA0999', 2, 2, 0, 6, 0 )
4284          ENDIF
4285       ENDDO
4286
4287
4288    END SUBROUTINE lsm_init
4289
4290
4291!------------------------------------------------------------------------------!
4292! Description:
4293! ------------
4294!> Allocate land surface model arrays and define pointers
4295!------------------------------------------------------------------------------!
4296    SUBROUTINE lsm_init_arrays
4297   
4298
4299       IMPLICIT NONE
4300
4301       INTEGER(iwp) ::  l !< index indicating facing of surface array
4302   
4303       ALLOCATE ( root_extr(nzb_soil:nzt_soil) )
4304       root_extr = 0.0_wp 
4305       
4306!
4307!--    Allocate surface and soil temperature / humidity. Please note,
4308!--    these arrays are allocated according to surface-data structure,
4309!--    even if they do not belong to the data type due to the
4310!--    pointer arithmetric (TARGET attribute is not allowed in a data-type).
4311#if defined( __nopointer )
4312!
4313!--    Horizontal surfaces
4314       ALLOCATE ( m_liq_h_p%var_1d(1:surf_lsm_h%ns)                      )
4315       ALLOCATE ( t_surface_h%var_1d(1:surf_lsm_h%ns)                    )
4316       ALLOCATE ( t_surface_h_p%var_1d(1:surf_lsm_h%ns)                  )
4317       ALLOCATE ( m_soil_h_p%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
4318       ALLOCATE ( t_soil_h_p%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
4319
4320!
4321!--    Vertical surfaces
4322       DO  l = 0, 3
4323          ALLOCATE ( m_liq_v(l)%var_1d(1:surf_lsm_v(l)%ns)                        )
4324          ALLOCATE ( m_liq_v_p(l)%var_1d(1:surf_lsm_v(l)%ns)                      )
4325          ALLOCATE ( t_surface_v(l)%var_1d(1:surf_lsm_v(l)%ns)                    )
4326          ALLOCATE ( t_surface_v_p(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
4327          ALLOCATE ( m_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )
4328          ALLOCATE ( m_soil_v_p(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)   )
4329          ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns)   )
4330          ALLOCATE ( t_soil_v_p(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
4331       ENDDO
4332!
4333!--    Allocate soil temperature and moisture. As these variables might be
4334!--    already allocated in case of restarts, check this.
4335       IF ( .NOT. ALLOCATED( m_liq_h%var_1d ) )                                &
4336          ALLOCATE ( m_liq_h%var_1d(1:surf_lsm_h%ns) )
4337       IF ( .NOT. ALLOCATED( m_soil_h%var_2d ) )                               &
4338          ALLOCATE ( m_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
4339       IF ( .NOT. ALLOCATED( t_soil_h%var_2d ) )                               &
4340          ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
4341
4342       DO  l = 0, 3
4343          IF ( .NOT. ALLOCATED( m_liq_v(l)%var_1d ) )                          &
4344             ALLOCATE ( m_liq_v(l)%var_1d(1:surf_lsm_v(l)%ns) )
4345          IF ( .NOT. ALLOCATED( m_soil_v(l)%var_2d ) )                         &
4346             ALLOCATE ( m_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
4347          IF ( .NOT. ALLOCATED( t_soil_v(l)%var_2d ) )                         &
4348             ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
4349       ENDDO
4350#else
4351!
4352!--    Horizontal surfaces
4353       ALLOCATE ( m_liq_h_1%var_1d(1:surf_lsm_h%ns)                      )
4354       ALLOCATE ( m_liq_h_2%var_1d(1:surf_lsm_h%ns)                      )
4355       ALLOCATE ( t_surface_h_1%var_1d(1:surf_lsm_h%ns)                  )
4356       ALLOCATE ( t_surface_h_2%var_1d(1:surf_lsm_h%ns)                  )
4357       ALLOCATE ( m_soil_h_1%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
4358       ALLOCATE ( m_soil_h_2%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
4359       ALLOCATE ( t_soil_h_1%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
4360       ALLOCATE ( t_soil_h_2%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
4361!
4362!--    Vertical surfaces
4363       DO  l = 0, 3
4364          ALLOCATE ( m_liq_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                      )
4365          ALLOCATE ( m_liq_v_2(l)%var_1d(1:surf_lsm_v(l)%ns)                      )
4366          ALLOCATE ( t_surface_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
4367          ALLOCATE ( t_surface_v_2(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
4368          ALLOCATE ( m_soil_v_1(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)   )
4369          ALLOCATE ( m_soil_v_2(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)   )
4370          ALLOCATE ( t_soil_v_1(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
4371          ALLOCATE ( t_soil_v_2(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
4372       ENDDO
4373#endif
4374!
4375!--    Allocate array for heat flux in W/m2, required for radiation?
4376!--    Consider to remove this array
4377       ALLOCATE( surf_lsm_h%surfhf(1:surf_lsm_h%ns) )
4378       DO  l = 0, 3
4379          ALLOCATE( surf_lsm_v(l)%surfhf(1:surf_lsm_v(l)%ns) )
4380       ENDDO
4381
4382
4383!
4384!--    Allocate intermediate timestep arrays
4385!--    Horizontal surfaces
4386       ALLOCATE ( tm_liq_h_m%var_1d(1:surf_lsm_h%ns)                     )
4387       ALLOCATE ( tt_surface_h_m%var_1d(1:surf_lsm_h%ns)                 )
4388       ALLOCATE ( tm_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )
4389       ALLOCATE ( tt_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  ) 
4390!
4391!--    Horizontal surfaces
4392       DO  l = 0, 3
4393          ALLOCATE ( tm_liq_v_m(l)%var_1d(1:surf_lsm_v(l)%ns)                     )
4394          ALLOCATE ( tt_surface_v_m(l)%var_1d(1:surf_lsm_v(l)%ns)                 )
4395          ALLOCATE ( tm_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
4396          ALLOCATE ( tt_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
4397       ENDDO
4398
4399!
4400!--    Allocate 2D vegetation model arrays
4401!--    Horizontal surfaces
4402       ALLOCATE ( surf_lsm_h%building_surface(1:surf_lsm_h%ns)    )
4403       ALLOCATE ( surf_lsm_h%c_liq(1:surf_lsm_h%ns)               )
4404       ALLOCATE ( surf_lsm_h%c_surface(1:surf_lsm_h%ns)           )
4405       ALLOCATE ( surf_lsm_h%c_veg(1:surf_lsm_h%ns)               )
4406       ALLOCATE ( surf_lsm_h%f_sw_in(1:surf_lsm_h%ns)             )
4407       ALLOCATE ( surf_lsm_h%ghf(1:surf_lsm_h%ns)                 )
4408       ALLOCATE ( surf_lsm_h%g_d(1:surf_lsm_h%ns)                 )
4409       ALLOCATE ( surf_lsm_h%lai(1:surf_lsm_h%ns)                 )
4410       ALLOCATE ( surf_lsm_h%lambda_surface_u(1:surf_lsm_h%ns)    )
4411       ALLOCATE ( surf_lsm_h%lambda_surface_s(1:surf_lsm_h%ns)    )
4412       ALLOCATE ( surf_lsm_h%nzt_pavement(1:surf_lsm_h%ns)        )
4413       ALLOCATE ( surf_lsm_h%pavement_surface(1:surf_lsm_h%ns)    )
4414       ALLOCATE ( surf_lsm_h%qsws_soil(1:surf_lsm_h%ns)           )
4415       ALLOCATE ( surf_lsm_h%qsws_liq(1:surf_lsm_h%ns)            )
4416       ALLOCATE ( surf_lsm_h%qsws_veg(1:surf_lsm_h%ns)            )
4417       ALLOCATE ( surf_lsm_h%rad_net_l(1:surf_lsm_h%ns)           ) 
4418       ALLOCATE ( surf_lsm_h%r_a(1:surf_lsm_h%ns)                 )
4419       ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns)            )
4420       ALLOCATE ( surf_lsm_h%r_soil(1:surf_lsm_h%ns)              )
4421       ALLOCATE ( surf_lsm_h%r_soil_min(1:surf_lsm_h%ns)          )
4422       ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns)                 )
4423       ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns)        )
4424       ALLOCATE ( surf_lsm_h%vegetation_surface(1:surf_lsm_h%ns)  )
4425       ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns)       )
4426
4427       surf_lsm_h%water_surface        = .FALSE.
4428       surf_lsm_h%pavement_surface     = .FALSE.
4429       surf_lsm_h%vegetation_surface   = .FALSE. 
4430!
4431!--    Vertical surfaces
4432       DO  l = 0, 3
4433          ALLOCATE ( surf_lsm_v(l)%building_surface(1:surf_lsm_v(l)%ns)    )
4434          ALLOCATE ( surf_lsm_v(l)%c_liq(1:surf_lsm_v(l)%ns)               )
4435          ALLOCATE ( surf_lsm_v(l)%c_surface(1:surf_lsm_v(l)%ns)           )
4436          ALLOCATE ( surf_lsm_v(l)%c_veg(1:surf_lsm_v(l)%ns)               )
4437          ALLOCATE ( surf_lsm_v(l)%f_sw_in(1:surf_lsm_v(l)%ns)             )
4438          ALLOCATE ( surf_lsm_v(l)%ghf(1:surf_lsm_v(l)%ns)                 )
4439          ALLOCATE ( surf_lsm_v(l)%g_d(1:surf_lsm_v(l)%ns)                 )
4440          ALLOCATE ( surf_lsm_v(l)%lai(1:surf_lsm_v(l)%ns)                 )
4441          ALLOCATE ( surf_lsm_v(l)%lambda_surface_u(1:surf_lsm_v(l)%ns)    )
4442          ALLOCATE ( surf_lsm_v(l)%lambda_surface_s(1:surf_lsm_v(l)%ns)    )
4443          ALLOCATE ( surf_lsm_v(l)%nzt_pavement(1:surf_lsm_v(l)%ns)        )
4444          ALLOCATE ( surf_lsm_v(l)%pavement_surface(1:surf_lsm_v(l)%ns)    )
4445          ALLOCATE ( surf_lsm_v(l)%qsws_soil(1:surf_lsm_v(l)%ns)           )
4446          ALLOCATE ( surf_lsm_v(l)%qsws_liq(1:surf_lsm_v(l)%ns)            )
4447          ALLOCATE ( surf_lsm_v(l)%qsws_veg(1:surf_lsm_v(l)%ns)            )
4448          ALLOCATE ( surf_lsm_v(l)%rad_net_l(1:surf_lsm_v(l)%ns)           )
4449          ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns)                 )
4450          ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns)            )
4451          ALLOCATE ( surf_lsm_v(l)%r_soil(1:surf_lsm_v(l)%ns)              )
4452          ALLOCATE ( surf_lsm_v(l)%r_soil_min(1:surf_lsm_v(l)%ns)          )
4453          ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns)                 )
4454          ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns)        )
4455          ALLOCATE ( surf_lsm_v(l)%vegetation_surface(1:surf_lsm_v(l)%ns)  )
4456          ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns)       )
4457
4458          surf_lsm_v(l)%water_surface     = .FALSE.
4459          surf_lsm_v(l)%pavement_surface  = .FALSE.
4460          surf_lsm_v(l)%vegetation_surface   = .FALSE. 
4461         
4462       ENDDO
4463
4464   
4465#if ! defined( __nopointer )
4466!
4467!--    Initial assignment of the pointers
4468!--    Horizontal surfaces
4469       t_soil_h    => t_soil_h_1;    t_soil_h_p    => t_soil_h_2
4470       t_surface_h => t_surface_h_1; t_surface_h_p => t_surface_h_2
4471       m_soil_h    => m_soil_h_1;    m_soil_h_p    => m_soil_h_2
4472       m_liq_h     => m_liq_h_1;     m_liq_h_p     => m_liq_h_2
4473!
4474!--    Vertical surfaces
4475       t_soil_v    => t_soil_v_1;    t_soil_v_p    => t_soil_v_2
4476       t_surface_v => t_surface_v_1; t_surface_v_p => t_surface_v_2
4477       m_soil_v    => m_soil_v_1;    m_soil_v_p    => m_soil_v_2
4478       m_liq_v     => m_liq_v_1;     m_liq_v_p     => m_liq_v_2
4479
4480#endif
4481
4482
4483    END SUBROUTINE lsm_init_arrays
4484
4485
4486!------------------------------------------------------------------------------!
4487! Description:
4488! ------------
4489!> Parin for &lsmpar for land surface model
4490!------------------------------------------------------------------------------!
4491    SUBROUTINE lsm_parin
4492
4493
4494       IMPLICIT NONE
4495
4496       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
4497       
4498       NAMELIST /lsm_par/         alpha_vangenuchten, c_surface,               &
4499                                  canopy_resistance_coefficient,               &
4500                                  constant_roughness,                          &
4501                                  conserve_water_content,                      &
4502                                  dz_soil,                                     &
4503                                  f_shortwave_incoming, field_capacity,        & 
4504                                  aero_resist_kray, hydraulic_conductivity,    &
4505                                  lambda_surface_stable,                       &
4506                                  lambda_surface_unstable, leaf_area_index,    &
4507                                  l_vangenuchten, min_canopy_resistance,       &
4508                                  min_soil_resistance, n_vangenuchten,         &
4509                                  pavement_depth_level,                        &
4510                                  pavement_heat_capacity,                      &
4511                                  pavement_heat_conduct, pavement_type,        &
4512                                  residual_moisture, root_fraction,            &
4513                                  saturation_moisture, skip_time_do_lsm,       &
4514                                  soil_moisture, soil_temperature,             &
4515                                  soil_type,                                   &
4516                                  surface_type,                                &
4517                                  vegetation_coverage, vegetation_type,        &
4518                                  water_temperature, water_type,               &
4519                                  wilting_point, z0_vegetation,                &
4520                                  z0h_vegetation, z0q_vegetation, z0_water,    &
4521                                  z0h_water, z0q_water, z0_pavement,           &
4522                                  z0h_pavement, z0q_pavement
4523       
4524       line = ' '
4525       
4526!
4527!--    Try to find land surface model package
4528       REWIND ( 11 )
4529       line = ' '
4530       DO   WHILE ( INDEX( line, '&lsm_par' ) == 0 )
4531          READ ( 11, '(A)', END=10 )  line
4532       ENDDO
4533       BACKSPACE ( 11 )
4534
4535!
4536!--    Read user-defined namelist
4537       READ ( 11, lsm_par )
4538
4539!
4540!--    Set flag that indicates that the land surface model is switched on
4541       land_surface = .TRUE.
4542
4543!
4544!--    Activate spinup
4545       IF ( spinup_time > 0.0_wp )  THEN
4546          coupling_start_time = spinup_time
4547          end_time = end_time + spinup_time
4548          IF ( spinup_pt_mean == 9999999.9_wp )  THEN
4549             spinup_pt_mean = pt_surface
4550          ENDIF
4551          spinup = .TRUE.
4552       ENDIF
4553
4554
4555 10    CONTINUE
4556       
4557
4558    END SUBROUTINE lsm_parin
4559
4560
4561!------------------------------------------------------------------------------!
4562! Description:
4563! ------------
4564!> Soil model as part of the land surface model. The model predicts soil
4565!> temperature and water content.
4566!------------------------------------------------------------------------------!
4567    SUBROUTINE lsm_soil_model( horizontal, l, calc_soil_moisture )
4568
4569
4570       IMPLICIT NONE
4571
4572       INTEGER(iwp) ::  k       !< running index
4573       INTEGER(iwp) ::  l       !< surface-data type index indication facing
4574       INTEGER(iwp) ::  m       !< running index
4575
4576       LOGICAL, INTENT(IN) ::  calc_soil_moisture !< flag indicating whether soil moisture shall be calculated or not.
4577
4578       LOGICAL      ::  horizontal !< flag indication horizontal wall, required to set pointer accordingly
4579
4580       REAL(wp)     ::  h_vg !< Van Genuchten coef. h
4581
4582       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
4583                                                 lambda_temp, & !< temp. lambda
4584                                                 tend           !< tendency
4585
4586       TYPE(surf_type_lsm), POINTER ::  surf_m_soil
4587       TYPE(surf_type_lsm), POINTER ::  surf_m_soil_p
4588       TYPE(surf_type_lsm), POINTER ::  surf_t_soil
4589       TYPE(surf_type_lsm), POINTER ::  surf_t_soil_p
4590       TYPE(surf_type_lsm), POINTER ::  surf_tm_soil_m
4591       TYPE(surf_type_lsm), POINTER ::  surf_tt_soil_m
4592
4593       TYPE(surf_type), POINTER  ::  surf  !< surface-date type variable
4594
4595       IF ( horizontal )  THEN
4596          surf           => surf_lsm_h
4597
4598          surf_m_soil    => m_soil_h
4599          surf_m_soil_p  => m_soil_h_p
4600          surf_t_soil    => t_soil_h
4601          surf_t_soil_p  => t_soil_h_p
4602          surf_tm_soil_m => tm_soil_h_m
4603          surf_tt_soil_m => tt_soil_h_m
4604       ELSE
4605          surf           => surf_lsm_v(l)
4606
4607          surf_m_soil    => m_soil_v(l)
4608          surf_m_soil_p  => m_soil_v_p(l)
4609          surf_t_soil    => t_soil_v(l)
4610          surf_t_soil_p  => t_soil_v_p(l)
4611          surf_tm_soil_m => tm_soil_v_m(l)
4612          surf_tt_soil_m => tt_soil_v_m(l)
4613       ENDIF
4614
4615       DO  m = 1, surf%ns
4616
4617          IF (  .NOT.  surf%water_surface(m) )  THEN
4618             DO  k = nzb_soil, nzt_soil
4619
4620                IF ( surf%pavement_surface(m)  .AND.                           &
4621                     k <= surf%nzt_pavement(m) )  THEN
4622                   
4623                   surf%rho_c_total(k,m) = surf%rho_c_total_def(k,m)
4624                   lambda_temp(k)        = surf%lambda_h_def(k,m) 
4625
4626                ELSE           
4627!
4628!--                Calculate volumetric heat capacity of the soil, taking
4629!--                into account water content
4630                   surf%rho_c_total(k,m) = (rho_c_soil *                       &
4631                                               ( 1.0_wp - surf%m_sat(k,m) )    &
4632                                               + rho_c_water * surf_m_soil%var_2d(k,m) )
4633
4634!
4635!--                Calculate soil heat conductivity at the center of the soil
4636!--                layers
4637                   lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(k,m)) *    &
4638                                  lambda_h_water ** surf_m_soil%var_2d(k,m)
4639
4640                   ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(k,m) / &
4641                                                     surf%m_sat(k,m) ) )
4642
4643                   lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +       &
4644                                    lambda_h_dry
4645                ENDIF
4646             ENDDO
4647
4648!
4649!--          Calculate soil heat conductivity (lambda_h) at the _layer level
4650!--          using linear interpolation. For pavement surface, the
4651!--          true pavement depth is considered
4652             DO  k = nzb_soil, nzt_soil-1
4653                   surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )  &
4654                                        * 0.5_wp
4655             ENDDO
4656             surf%lambda_h(nzt_soil,m) = lambda_temp(nzt_soil)
4657
4658!
4659!--          Prognostic equation for soil temperature t_soil
4660             tend(:) = 0.0_wp
4661
4662             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *            &
4663                    ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) &
4664                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil_center(nzb_soil) &
4665                      + surf%ghf(m) ) * ddz_soil(nzb_soil)
4666
4667             DO  k = nzb_soil+1, nzt_soil
4668                tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) )                   &
4669                          * (   surf%lambda_h(k,m)                             &
4670                     * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
4671                     * ddz_soil_center(k)                                      &
4672                     - surf%lambda_h