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

Last change on this file since 2707 was 2707, checked in by suehring, 4 years ago

changes documented

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