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

Last change on this file since 4194 was 4194, checked in by suehring, 3 years ago

In case z0 exceeds the surface-layer height at water surfaces, apply more strict limitation of z0, in order to avoid instabilities when friction velocity is almost zero.

  • Property svn:keywords set to Id
File size: 327.8 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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: land_surface_model_mod.f90 4194 2019-08-28 08:09:44Z suehring $
27! Apply more strict limitation of z0 over water surfaces in case it exceeds the
28! surface-layer height, in order to avoid instabilities.
29!
30! 4188 2019-08-26 14:15:47Z suehring
31! Minor adjustment in error numbers, typos corrected
32!
33! 4187 2019-08-26 12:43:15Z suehring
34! Adjust message call in case of local checks
35!
36! 4182 2019-08-22 15:20:23Z scharf
37! Corrected "Former revisions" section
38!
39! 4118 2019-07-25 16:11:45Z suehring
40! Initialization of soil temperature and moisture via dynamic input file only
41! for vegetation and pavement surfaces.
42!
43! 4110 2019-07-22 17:05:21Z suehring
44! Relax checks for non-consistent initialization in case static or dynamic
45! input is provided. For example, soil_temperature or deep_soil_temperature
46! is not mandatory any more if dynamic input is available. Also, improper
47! settings of x_type in namelist are only checked if no static file is
48! available.
49!
50! 4109 2019-07-22 17:00:34Z suehring
51! Further revision of last commit in order to avoid any side effects when
52! albedo type is not set in namelist and default albedo type changes.
53!
54! 4024 2019-06-12 14:06:46Z suehring
55! Bugfix in albedo initialization, caused crashes in rrtmg calls
56!
57! 3987 2019-05-22 09:52:13Z kanani
58! Introduce alternative switch for debug output during timestepping
59!
60! 3964 2019-05-09 09:48:32Z suehring
61! In a nested child domain, distinguish between soil moisture and temperature
62! initialization from parent via dynamic input file. Further, initialize soil
63! moisture/temperature from dynamic input file only when initialization via
64! 'inifor' is desired.
65!
66! 3943 2019-05-02 09:50:41Z maronga
67! Removed extra blank character
68!
69! 3941 2019-04-30 09:48:33Z suehring
70! Check that at least one surface type is set at surface element.
71!
72! 3933 2019-04-25 12:33:20Z kanani
73! Remove unused subroutine and allocation of pt_2m, this is done in surface_mod
74! now (surfaces%pt_2m)
75!
76!
77! Changes related to global restructuring of location messages and introduction
78! of additional debug messages
79!
80! 3881 2019-04-10 09:31:22Z suehring
81! Bugfix in level 3 initialization of pavement albedo type and pavement
82! emissivity
83!
84! 3868 2019-04-08 11:52:36Z suehring
85! More strict limitation of roughness length when it is in the order of the
86! vertical grid spacing
87!
88! 3856 2019-04-03 11:06:59Z suehring
89! Bugfix in lsm_init in case no surface-fractions are provided
90!
91! 3847 2019-04-01 14:51:44Z suehring
92! Adjust message-call for checks that are especially carried out locally.
93!
94! 3832 2019-03-28 13:16:58Z raasch
95! instrumented with openmp directives
96!
97! 3786 2019-03-06 16:58:03Z raasch
98! further unused variables removed
99!
100! 3767 2019-02-27 08:18:02Z raasch
101! unused variable for file index removed from rrd-subroutines parameter list
102!
103! 3715 2019-02-04 17:34:55Z suehring
104! Revise check for saturation moisture
105!
106! 3710 2019-01-30 18:11:19Z suehring
107! Check if soil-, water-, pavement- and vegetation types are set within a valid
108! range.
109!
110! 3692 2019-01-23 14:45:49Z suehring
111! Revise check for soil moisture higher than its saturation value
112!
113! 3685 2019-01-21 01:02:11Z knoop
114! Some interface calls moved to module_interface + cleanup
115!
116! 3677 2019-01-17 09:07:06Z moh.hefny
117! Removed most_method
118!
119! 3655 2019-01-07 16:51:22Z knoop
120! nopointer option removed
121!
122! 1496 2014-12-02 17:25:50Z maronga
123! Initial revision
124!
125!
126! Description:
127! ------------
128!> Land surface model, consisting of a solver for the energy balance at the
129!> surface and a multi layer soil scheme. The scheme is similar to the TESSEL
130!> scheme implemented in the ECMWF IFS model, with modifications according to
131!> H-TESSEL. The implementation is based on the formulation implemented in the
132!> DALES and UCLA-LES models.
133!>
134!> @todo Extensive verification energy-balance solver for vertical surfaces,
135!>       e.g. parametrization of r_a
136!> @todo Revise single land-surface processes for vertical surfaces, e.g.
137!>       treatment of humidity, etc.
138!> @todo Consider partial absorption of the net shortwave radiation by the
139!>       skin layer.
140!> @todo Improve surface water parameterization
141!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
142!>       nzt_soil=3)).
143!> @todo Implement surface runoff model (required when performing long-term LES
144!>       with considerable precipitation.
145!> @todo Revise calculation of f2 when wilting point is non-constant in the
146!>       soil
147!> @todo Allow for zero soil moisture (currently, it is set to wilting point)
148!> @note No time step criterion is required as long as the soil layers do not
149!>       become too thin.
150!> @todo Attention, pavement_subpars_1/2 are hardcoded to 8 levels, in case
151!>       more levels are used this may cause an potential bug
152!> @todo Routine calc_q_surface required?
153!> @todo Allow for precipitation water to enter pavements that are semi-pervious
154!------------------------------------------------------------------------------!
155 MODULE land_surface_model_mod
156 
157    USE arrays_3d,                                                             &
158        ONLY:  hyp, pt, prr, q, q_p, ql, vpt, u, v, w, hyrho, exner, d_exner
159
160    USE basic_constants_and_equations_mod,                                     &
161        ONLY:  c_p, g, lv_d_cp, l_v, kappa, magnus, rho_l, r_d, r_v, rd_d_rv
162
163    USE calc_mean_profile_mod,                                                 &
164        ONLY:  calc_mean_profile
165
166    USE control_parameters,                                                    &
167        ONLY:  cloud_droplets, coupling_start_time,                            &
168               debug_output, debug_output_timestep, debug_string,              &
169               dt_3d,                                                          &
170               end_time, humidity, intermediate_timestep_count,                &
171               initializing_actions, intermediate_timestep_count_max,          &
172               land_surface, max_masks, pt_surface,                            &
173               rho_surface, spinup, spinup_pt_mean, spinup_time,               &
174               surface_pressure, timestep_scheme, tsc,                         &
175               time_since_reference_point
176
177    USE indices,                                                               &
178        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb
179
180    USE bulk_cloud_model_mod,                                                  &
181        ONLY: bulk_cloud_model, precipitation
182
183    USE netcdf_data_input_mod,                                                 &
184        ONLY :  building_type_f, init_3d,                                      &
185                input_pids_dynamic,                                            &
186                input_pids_static,                                             &
187                netcdf_data_input_interpolate, netcdf_data_input_init_lsm,     &
188                pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,  &
189                root_area_density_lsm_f, soil_pars_f, soil_type_f,             &
190                surface_fraction_f, vegetation_pars_f, vegetation_type_f,      &
191                water_pars_f, water_type_f
192
193    USE kinds
194
195    USE pegrid
196
197    USE radiation_model_mod,                                                   &
198        ONLY:  albedo, albedo_type, emissivity, force_radiation_call,          &
199               radiation, radiation_scheme, unscheduled_radiation_calls
200       
201    USE statistics,                                                            &
202        ONLY:  hom, statistic_regions
203
204    USE surface_mod,                                                           &
205        ONLY :  ind_pav_green, ind_veg_wall, ind_wat_win,                      &
206                surf_lsm_h, surf_lsm_v, surf_type, surface_restore_elements
207
208    IMPLICIT NONE
209
210    TYPE surf_type_lsm
211       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  var_1d !< 1D prognostic variable
212       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2d !< 2D prognostic variable
213    END TYPE surf_type_lsm
214
215!
216!-- LSM model constants
217
218    REAL(wp), PARAMETER  ::                    &
219              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
220              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil (W/m/K) 
221              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix (W/m/K)
222              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water (W/m/K)
223              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
224              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil (J/m3/K)
225              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water (J/m3/K)
226              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
227
228
229    REAL(wp), DIMENSION(0:7), PARAMETER  :: dz_soil_default =                  & ! default soil layer configuration
230                                            (/ 0.01_wp, 0.02_wp, 0.04_wp,      &
231                                               0.06_wp, 0.14_wp, 0.26_wp,      &
232                                               0.54_wp, 1.86_wp/)
233
234    REAL(wp), DIMENSION(0:3), PARAMETER  :: dz_soil_ref =                      & ! reference four layer soil configuration used for estimating the root fractions
235                                            (/ 0.07_wp, 0.21_wp, 0.72_wp,      &
236                                               1.89_wp /)
237
238    REAL(wp), DIMENSION(0:3), PARAMETER  :: zs_ref =                           & ! reference four layer soil configuration used for estimating the root fractions
239                                            (/ 0.07_wp, 0.28_wp, 1.0_wp,       &
240                                               2.89_wp /)
241
242
243!
244!-- LSM variables
245    CHARACTER(10) :: surface_type = 'netcdf'      !< general classification. Allowed are:
246                                                  !< 'vegetation', 'pavement', ('building'),
247                                                  !< 'water', and 'netcdf'
248
249
250
251    INTEGER(iwp) :: nzb_soil = 0,             & !< bottom of the soil model (Earth's surface)
252                    nzt_soil = 7,             & !< top of the soil model
253                    nzt_pavement = 0,         & !< top of the pavement within the soil
254                    nzs = 8,                  & !< number of soil layers
255                    pavement_depth_level = 0, & !< default NAMELIST nzt_pavement
256                    pavement_type = 1,        & !< default NAMELIST pavement_type                 
257                    soil_type = 3,            & !< default NAMELIST soil_type
258                    vegetation_type = 2,      & !< default NAMELIST vegetation_type
259                    water_type = 1              !< default NAMELISt water_type
260                   
261   
262       
263    LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
264               constant_roughness = .FALSE.,     & !< use fixed/dynamic roughness lengths for water surfaces
265               force_radiation_call_l = .FALSE., & !< flag to force calling of radiation routine
266               aero_resist_kray = .TRUE.           !< flag to control parametrization of aerodynamic resistance at vertical surface elements
267
268!   value 9999999.9_wp -> generic available or user-defined value must be set
269!   otherwise -> no generic variable and user setting is optional
270    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
271                canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
272                c_surface = 9999999.9_wp,               & !< Surface (skin) heat capacity (J/m2/K)
273                deep_soil_temperature =  9999999.9_wp,  & !< Deep soil temperature (bottom boundary condition)
274                drho_l_lv,                              & !< (rho_l * l_v)**-1
275                field_capacity = 9999999.9_wp,          & !< NAMELIST m_fc
276                f_shortwave_incoming = 9999999.9_wp,    & !< NAMELIST f_sw_in
277                hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
278                ke = 0.0_wp,                            & !< Kersten number
279                lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil (W/m/K)
280                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s (W/m2/K)
281                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u (W/m2/K)
282                leaf_area_index = 9999999.9_wp,         & !< NAMELIST lai
283                l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
284                min_canopy_resistance = 9999999.9_wp,   & !< NAMELIST r_canopy_min
285                min_soil_resistance = 50.0_wp,          & !< NAMELIST r_soil_min
286                m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
287                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
288                pavement_heat_capacity = 9999999.9_wp,  & !< volumetric heat capacity of pavement (e.g. roads) (J/m3/K)
289                pavement_heat_conduct  = 9999999.9_wp,  & !< heat conductivity for pavements (e.g. roads) (W/m/K)
290                q_s = 0.0_wp,                           & !< saturation water vapor mixing ratio
291                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
292                rho_cp,                                 & !< rho_surface * cp
293                rho_lv,                                 & !< rho_ocean * l_v
294                saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
295                skip_time_do_lsm = 0.0_wp,              & !< LSM is not called before this time
296                vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
297                water_temperature = 9999999.9_wp,       & !< water temperature
298                wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
299                z0_vegetation  = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
300                z0h_vegetation = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
301                z0q_vegetation = 9999999.9_wp,          & !< NAMELIST z0q (lsm_par)
302                z0_pavement    = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
303                z0h_pavement   = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
304                z0q_pavement   = 9999999.9_wp,          & !< NAMELIST z0q (lsm_par)
305                z0_water       = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
306                z0h_water      = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
307                z0q_water      = 9999999.9_wp             !< NAMELIST z0q (lsm_par) 
308               
309               
310    REAL(wp), DIMENSION(:), ALLOCATABLE  :: ddz_soil_center, & !< 1/dz_soil_center
311                                            ddz_soil,        & !< 1/dz_soil
312                                            dz_soil_center,  & !< soil grid spacing (center-center)
313                                            zs,              & !< depth of the temperature/moisute levels
314                                            root_extr          !< root extraction
315
316
317                                           
318    REAL(wp), DIMENSION(0:20)  ::  root_fraction = 9999999.9_wp,     & !< (NAMELIST) distribution of root surface area to the individual soil layers
319                                   soil_moisture = 0.0_wp,           & !< NAMELIST soil moisture content (m3/m3)
320                                   soil_temperature = 9999999.9_wp,  & !< NAMELIST soil temperature (K) +1
321                                   dz_soil  = 9999999.9_wp,          & !< (NAMELIST) soil layer depths (spacing)
322                                   zs_layer = 9999999.9_wp             !< soil layer depths (edge)
323                                 
324    TYPE(surf_type_lsm), POINTER ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
325                                     t_soil_h_p,  & !< Prog. soil temperature (K), horizontal surface elements
326                                     m_soil_h,    & !< Soil moisture (m3/m3), horizontal surface elements
327                                     m_soil_h_p     !< Prog. soil moisture (m3/m3), horizontal surface elements
328
329    TYPE(surf_type_lsm), TARGET  ::  t_soil_h_1,  & !<
330                                     t_soil_h_2,  & !<
331                                     m_soil_h_1,  & !<
332                                     m_soil_h_2     !<
333
334    TYPE(surf_type_lsm), DIMENSION(:), POINTER :: &
335                                     t_soil_v,    & !< Soil temperature (K), vertical surface elements
336                                     t_soil_v_p,  & !< Prog. soil temperature (K), vertical surface elements
337                                     m_soil_v,    & !< Soil moisture (m3/m3), vertical surface elements
338                                     m_soil_v_p     !< Prog. soil moisture (m3/m3), vertical surface elements   
339
340    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::&
341                                     t_soil_v_1,  & !<
342                                     t_soil_v_2,  & !<
343                                     m_soil_v_1,  & !<
344                                     m_soil_v_2     !<
345
346    TYPE(surf_type_lsm), POINTER  ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements
347                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements
348                                      m_liq_h,        & !< liquid water reservoir (m), horizontal surface elements
349                                      m_liq_h_p         !< progn. liquid water reservoir (m), horizontal surface elements
350
351    TYPE(surf_type_lsm), TARGET   ::  t_surface_h_1,  & !<
352                                      t_surface_h_2,  & !<
353                                      m_liq_h_1,      & !<
354                                      m_liq_h_2         !<
355
356    TYPE(surf_type_lsm), DIMENSION(:), POINTER  ::    &
357                                      t_surface_v,    & !< surface temperature (K), vertical surface elements
358                                      t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements
359                                      m_liq_v,        & !< liquid water reservoir (m), vertical surface elements
360                                      m_liq_v_p         !< progn. liquid water reservoir (m), vertical surface elements
361
362    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
363                                      t_surface_v_1,  & !<
364                                      t_surface_v_2,  & !<
365                                      m_liq_v_1,      & !<
366                                      m_liq_v_2         !<
367
368    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av
369
370    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  t_soil_av, & !< Average of t_soil
371                                                        m_soil_av    !< Average of m_soil
372
373    TYPE(surf_type_lsm), TARGET ::  tm_liq_h_m      !< liquid water reservoir tendency (m), horizontal surface elements
374    TYPE(surf_type_lsm), TARGET ::  tt_surface_h_m  !< surface temperature tendency (K), horizontal surface elements
375    TYPE(surf_type_lsm), TARGET ::  tt_soil_h_m     !< t_soil storage array, horizontal surface elements
376    TYPE(surf_type_lsm), TARGET ::  tm_soil_h_m     !< m_soil storage array, horizontal surface elements
377
378    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_liq_v_m      !< liquid water reservoir tendency (m), vertical surface elements
379    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_surface_v_m  !< surface temperature tendency (K), vertical surface elements
380    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_soil_v_m     !< t_soil storage array, vertical surface elements
381    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_soil_v_m     !< m_soil storage array, vertical surface elements
382
383!
384!-- Energy balance variables               
385    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
386              c_liq_av,         & !< average of c_liq
387              c_soil_av,        & !< average of c_soil
388              c_veg_av,         & !< average of c_veg
389              lai_av,           & !< average of lai       
390              qsws_liq_av,      & !< average of qsws_liq
391              qsws_soil_av,     & !< average of qsws_soil
392              qsws_veg_av,      & !< average of qsws_veg
393              r_s_av              !< average of r_s
394 
395!
396!-- Predefined Land surface classes (vegetation_type)
397    CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ &
398                                   'user defined              ',           & !  0
399                                   'bare soil                 ',           & !  1                           
400                                   'crops, mixed farming      ',           & !  2
401                                   'short grass               ',           & !  3
402                                   'evergreen needleleaf trees',           & !  4
403                                   'deciduous needleleaf trees',           & !  5
404                                   'evergreen broadleaf trees ',           & !  6
405                                   'deciduous broadleaf trees ',           & !  7
406                                   'tall grass                ',           & !  8
407                                   'desert                    ',           & !  9
408                                   'tundra                    ',           & ! 10
409                                   'irrigated crops           ',           & ! 11
410                                   'semidesert                ',           & ! 12
411                                   'ice caps and glaciers     ',           & ! 13
412                                   'bogs and marshes          ',           & ! 14
413                                   'evergreen shrubs          ',           & ! 15
414                                   'deciduous shrubs          ',           & ! 16
415                                   'mixed forest/woodland     ',           & ! 17
416                                   'interrupted forest        '            & ! 18
417                                                                 /)
418
419!
420!-- Soil model classes (soil_type)
421    CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ &
422                                   'user defined',                  & ! 0
423                                   'coarse      ',                  & ! 1
424                                   'medium      ',                  & ! 2
425                                   'medium-fine ',                  & ! 3
426                                   'fine        ',                  & ! 4
427                                   'very fine   ',                  & ! 5
428                                   'organic     '                   & ! 6
429                                                                 /)
430
431!
432!-- Pavement classes
433    CHARACTER(29), DIMENSION(0:15), PARAMETER :: pavement_type_name = (/ &
434                                   'user defined                 ', & ! 0
435                                   'asphalt/concrete mix         ', & ! 1
436                                   'asphalt (asphalt concrete)   ', & ! 2
437                                   'concrete (Portland concrete) ', & ! 3
438                                   'sett                         ', & ! 4
439                                   'paving stones                ', & ! 5
440                                   'cobblestone                  ', & ! 6
441                                   'metal                        ', & ! 7
442                                   'wood                         ', & ! 8
443                                   'gravel                       ', & ! 9
444                                   'fine gravel                  ', & ! 10
445                                   'pebblestone                  ', & ! 11
446                                   'woodchips                    ', & ! 12
447                                   'tartan (sports)              ', & ! 13
448                                   'artifical turf (sports)      ', & ! 14
449                                   'clay (sports)                '  & ! 15
450                                                                 /)                                                             
451                                                                 
452!
453!-- Water classes
454    CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ &
455                                   'user defined',                   & ! 0
456                                   'lake        ',                   & ! 1
457                                   'river       ',                   & ! 2
458                                   'ocean       ',                   & ! 3
459                                   'pond        ',                   & ! 4
460                                   'fountain    '                    & ! 5
461                                                                  /)                                                                                 
462                   
463!
464!-- Land surface parameters according to the respective classes (vegetation_type)
465    INTEGER(iwp) ::  ind_v_rc_min = 0    !< index for r_canopy_min in vegetation_pars
466    INTEGER(iwp) ::  ind_v_rc_lai = 1    !< index for LAI in vegetation_pars
467    INTEGER(iwp) ::  ind_v_c_veg   = 2   !< index for c_veg in vegetation_pars
468    INTEGER(iwp) ::  ind_v_gd  = 3       !< index for g_d in vegetation_pars
469    INTEGER(iwp) ::  ind_v_z0 = 4        !< index for z0 in vegetation_pars
470    INTEGER(iwp) ::  ind_v_z0qh = 5      !< index for z0h / z0q in vegetation_pars
471    INTEGER(iwp) ::  ind_v_lambda_s = 6  !< index for lambda_s_s in vegetation_pars
472    INTEGER(iwp) ::  ind_v_lambda_u = 7  !< index for lambda_s_u in vegetation_pars
473    INTEGER(iwp) ::  ind_v_f_sw_in = 8   !< index for f_sw_in in vegetation_pars
474    INTEGER(iwp) ::  ind_v_c_surf = 9    !< index for c_surface in vegetation_pars
475    INTEGER(iwp) ::  ind_v_at = 10       !< index for albedo_type in vegetation_pars
476    INTEGER(iwp) ::  ind_v_emis = 11     !< index for emissivity in vegetation_pars
477
478    INTEGER(iwp) ::  ind_w_temp     = 0    !< index for temperature in water_pars
479    INTEGER(iwp) ::  ind_w_z0       = 1    !< index for z0 in water_pars
480    INTEGER(iwp) ::  ind_w_z0h      = 2    !< index for z0h in water_pars
481    INTEGER(iwp) ::  ind_w_lambda_s = 3    !< index for lambda_s_s in water_pars
482    INTEGER(iwp) ::  ind_w_lambda_u = 4    !< index for lambda_s_u in water_pars
483    INTEGER(iwp) ::  ind_w_at       = 5    !< index for albedo type in water_pars
484    INTEGER(iwp) ::  ind_w_emis     = 6    !< index for emissivity in water_pars
485
486    INTEGER(iwp) ::  ind_p_z0       = 0    !< index for z0 in pavement_pars
487    INTEGER(iwp) ::  ind_p_z0h      = 1    !< index for z0h in pavement_pars
488    INTEGER(iwp) ::  ind_p_at       = 2    !< index for albedo type in pavement_pars
489    INTEGER(iwp) ::  ind_p_emis     = 3    !< index for emissivity in pavement_pars
490    INTEGER(iwp) ::  ind_p_lambda_h = 0    !< index for lambda_h in pavement_subsurface_pars
491    INTEGER(iwp) ::  ind_p_rho_c    = 1    !< index for rho_c in pavement_pars
492!
493!-- Land surface parameters
494!-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface, albedo_type, emissivity
495    REAL(wp), DIMENSION(0:11,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
496          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
497        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
498        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
499        500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  5.0_wp, 0.97_wp, & !  4
500        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
501        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
502        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
503        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
504        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
505         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
506        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
507        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
508          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
509        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
510        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
511        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
512        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
513        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
514                                                               /), (/ 12, 18 /) )
515
516                                   
517!
518!-- Root distribution for default soil layer configuration (sum = 1)
519!--                                level 1 - level 4 according to zs_ref
520    REAL(wp), DIMENSION(0:3,1:18), PARAMETER :: root_distribution = RESHAPE( (/ &
521                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  1
522                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  2
523                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  3
524                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  4
525                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  5
526                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  6
527                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  7
528                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  8
529                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  9
530                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & ! 10
531                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 11
532                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 12
533                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 13
534                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 14
535                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 15
536                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16
537                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 17
538                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 18
539                                 /), (/ 4, 18 /) )
540
541!
542!-- Soil parameters according to the following porosity classes (soil_type)
543
544!
545!-- Soil parameters  alpha_vg,      l_vg,    n_vg, gamma_w_sat,    m_sat,     m_fc,   m_wilt,    m_res
546    REAL(wp), DIMENSION(0:7,1:6), PARAMETER :: soil_pars = RESHAPE( (/     &
547                      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
548                      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
549                      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
550                      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
551                      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
552                      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
553                                                                     /), (/ 8, 6 /) )
554
555
556!
557!-- TO BE FILLED
558!-- Pavement parameters      z0,       z0h, albedo_type, emissivity 
559    REAL(wp), DIMENSION(0:3,1:15), PARAMETER :: pavement_pars = RESHAPE( (/ &
560                      5.0E-2_wp, 5.0E-4_wp,     18.0_wp,    0.97_wp,  & !  1
561                      5.0E-2_wp, 5.0E-4_wp,     19.0_wp,    0.94_wp,  & !  2
562                      1.0E-2_wp, 1.0E-4_wp,     20.0_wp,    0.98_wp,  & !  3                                 
563                      1.0E-2_wp, 1.0E-4_wp,     21.0_wp,    0.93_wp,  & !  4
564                      1.0E-2_wp, 1.0E-4_wp,     22.0_wp,    0.97_wp,  & !  5
565                      1.0E-2_wp, 1.0E-4_wp,     23.0_wp,    0.97_wp,  & !  6
566                      1.0E-2_wp, 1.0E-4_wp,     24.0_wp,    0.97_wp,  & !  7
567                      1.0E-2_wp, 1.0E-4_wp,     25.0_wp,    0.94_wp,  & !  8
568                      1.0E-2_wp, 1.0E-4_wp,     26.0_wp,    0.98_wp,  & !  9                                 
569                      1.0E-2_wp, 1.0E-4_wp,     27.0_wp,    0.93_wp,  & ! 10
570                      1.0E-2_wp, 1.0E-4_wp,     28.0_wp,    0.97_wp,  & ! 11
571                      1.0E-2_wp, 1.0E-4_wp,     29.0_wp,    0.97_wp,  & ! 12
572                      1.0E-2_wp, 1.0E-4_wp,     30.0_wp,    0.97_wp,  & ! 13
573                      1.0E-2_wp, 1.0E-4_wp,     31.0_wp,    0.94_wp,  & ! 14
574                      1.0E-2_wp, 1.0E-4_wp,     32.0_wp,    0.98_wp   & ! 15
575                      /), (/ 4, 15 /) )                             
576!
577!-- Pavement subsurface parameters part 1: thermal conductivity (W/m/K)
578!--   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
579    REAL(wp), DIMENSION(0:7,1:15), PARAMETER :: pavement_subsurface_pars_1 = RESHAPE( (/ &
580       0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp, 9999999.9_wp, 9999999.9_wp, & !  1
581       0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp,   0.75_wp, 9999999.9_wp, 9999999.9_wp, & !  2
582       0.89_wp,   0.89_wp,   0.89_wp,   0.89_wp,   0.89_wp,   0.89_wp, 9999999.9_wp, 9999999.9_wp, & !  3
583       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
584       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
585       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
586       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
587       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
588       0.70_wp,   0.70_wp,   0.70_wp,   0.70_wp,   0.70_wp,   0.70_wp, 9999999.9_wp, 9999999.9_wp, & !  9
589       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
590       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
591       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
592       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
593       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
594       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
595       /), (/ 8, 15 /) )
596
597!
598!-- Pavement subsurface parameters part 2: volumetric heat capacity (J/m3/K)
599!--     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
600    REAL(wp), DIMENSION(0:7,1:15), PARAMETER :: pavement_subsurface_pars_2 = RESHAPE( (/ &
601       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
602       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
603       1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 1.76E6_wp, 9999999.9_wp, 9999999.9_wp, & !  3
604       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
605       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
606       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
607       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
608       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
609       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
610       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
611       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
612       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
613       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
614       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
615       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
616                           /), (/ 8, 15 /) )
617 
618!
619!-- TO BE FILLED
620!-- Water parameters                    temperature,     z0,      z0h, albedo_type, emissivity,
621    REAL(wp), DIMENSION(0:6,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
622       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 1
623       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 2
624       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 3
625       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 4
626       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp  & ! 5
627                                                                     /), (/ 7, 5 /) )                                                                   
628                                                                                                                                     
629    SAVE
630
631
632    PRIVATE
633
634   
635!
636!-- Public functions
637    PUBLIC lsm_boundary_condition, lsm_check_data_output,                      &
638           lsm_check_data_output_pr,                                           &
639           lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,& 
640           lsm_data_output_2d, lsm_data_output_3d, lsm_energy_balance,         &
641           lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model,   &
642           lsm_swap_timelevel, lsm_rrd_local, lsm_wrd_local
643! !vegetat
644!-- Public parameters, constants and initial values
645    PUBLIC aero_resist_kray, skip_time_do_lsm
646
647!
648!-- Public grid variables
649    PUBLIC nzb_soil, nzs, nzt_soil, zs
650
651!
652!-- Public prognostic variables
653    PUBLIC m_soil_h, t_soil_h
654
655    INTERFACE lsm_boundary_condition
656       MODULE PROCEDURE lsm_boundary_condition
657    END INTERFACE lsm_boundary_condition
658
659    INTERFACE lsm_check_data_output
660       MODULE PROCEDURE lsm_check_data_output
661    END INTERFACE lsm_check_data_output
662   
663    INTERFACE lsm_check_data_output_pr
664       MODULE PROCEDURE lsm_check_data_output_pr
665    END INTERFACE lsm_check_data_output_pr
666   
667    INTERFACE lsm_check_parameters
668       MODULE PROCEDURE lsm_check_parameters
669    END INTERFACE lsm_check_parameters
670   
671    INTERFACE lsm_3d_data_averaging
672       MODULE PROCEDURE lsm_3d_data_averaging
673    END INTERFACE lsm_3d_data_averaging
674
675    INTERFACE lsm_data_output_2d
676       MODULE PROCEDURE lsm_data_output_2d
677    END INTERFACE lsm_data_output_2d
678
679    INTERFACE lsm_data_output_3d
680       MODULE PROCEDURE lsm_data_output_3d
681    END INTERFACE lsm_data_output_3d
682
683    INTERFACE lsm_define_netcdf_grid
684       MODULE PROCEDURE lsm_define_netcdf_grid
685    END INTERFACE lsm_define_netcdf_grid
686
687    INTERFACE lsm_energy_balance
688       MODULE PROCEDURE lsm_energy_balance
689    END INTERFACE lsm_energy_balance
690
691    INTERFACE lsm_header
692       MODULE PROCEDURE lsm_header
693    END INTERFACE lsm_header
694   
695    INTERFACE lsm_init
696       MODULE PROCEDURE lsm_init
697    END INTERFACE lsm_init
698
699    INTERFACE lsm_init_arrays
700       MODULE PROCEDURE lsm_init_arrays
701    END INTERFACE lsm_init_arrays
702   
703    INTERFACE lsm_parin
704       MODULE PROCEDURE lsm_parin
705    END INTERFACE lsm_parin
706   
707    INTERFACE lsm_soil_model
708       MODULE PROCEDURE lsm_soil_model
709    END INTERFACE lsm_soil_model
710
711    INTERFACE lsm_swap_timelevel
712       MODULE PROCEDURE lsm_swap_timelevel
713    END INTERFACE lsm_swap_timelevel
714
715    INTERFACE lsm_rrd_local
716       MODULE PROCEDURE lsm_rrd_local
717    END INTERFACE lsm_rrd_local
718
719    INTERFACE lsm_wrd_local
720       MODULE PROCEDURE lsm_wrd_local
721    END INTERFACE lsm_wrd_local
722
723 CONTAINS
724
725
726!------------------------------------------------------------------------------!
727! Description:
728! ------------
729!> Set internal Neumann boundary condition at outer soil grid points
730!> for temperature and humidity.
731!------------------------------------------------------------------------------!
732 SUBROUTINE lsm_boundary_condition
733 
734    IMPLICIT NONE
735
736    INTEGER(iwp) :: i      !< grid index x-direction
737    INTEGER(iwp) :: ioff   !< offset index x-direction indicating location of soil grid point
738    INTEGER(iwp) :: j      !< grid index y-direction
739    INTEGER(iwp) :: joff   !< offset index x-direction indicating location of soil grid point
740    INTEGER(iwp) :: k      !< grid index z-direction
741    INTEGER(iwp) :: koff   !< offset index x-direction indicating location of soil grid point
742    INTEGER(iwp) :: l      !< running index surface-orientation
743    INTEGER(iwp) :: m      !< running index surface elements
744
745    koff = surf_lsm_h%koff
746    DO  m = 1, surf_lsm_h%ns
747       i = surf_lsm_h%i(m)
748       j = surf_lsm_h%j(m)
749       k = surf_lsm_h%k(m)
750       pt(k+koff,j,i) = pt(k,j,i)
751    ENDDO
752
753    DO  l = 0, 3
754       ioff = surf_lsm_v(l)%ioff
755       joff = surf_lsm_v(l)%joff
756       DO  m = 1, surf_lsm_v(l)%ns
757          i = surf_lsm_v(l)%i(m)
758          j = surf_lsm_v(l)%j(m)
759          k = surf_lsm_v(l)%k(m)
760          pt(k,j+joff,i+ioff) = pt(k,j,i)
761       ENDDO
762    ENDDO
763!
764!-- In case of humidity, set boundary conditions also for q and vpt.
765    IF ( humidity )  THEN
766       koff = surf_lsm_h%koff
767       DO  m = 1, surf_lsm_h%ns
768          i = surf_lsm_h%i(m)
769          j = surf_lsm_h%j(m)
770          k = surf_lsm_h%k(m)
771          q(k+koff,j,i)   = q(k,j,i)
772          vpt(k+koff,j,i) = vpt(k,j,i)
773       ENDDO
774
775       DO  l = 0, 3
776          ioff = surf_lsm_v(l)%ioff
777          joff = surf_lsm_v(l)%joff
778          DO  m = 1, surf_lsm_v(l)%ns
779             i = surf_lsm_v(l)%i(m)
780             j = surf_lsm_v(l)%j(m)
781             k = surf_lsm_v(l)%k(m)
782             q(k,j+joff,i+ioff)   = q(k,j,i)
783             vpt(k,j+joff,i+ioff) = vpt(k,j,i)
784          ENDDO
785       ENDDO
786    ENDIF
787
788 END SUBROUTINE lsm_boundary_condition
789
790!------------------------------------------------------------------------------!
791! Description:
792! ------------
793!> Check data output for land surface model
794!------------------------------------------------------------------------------!
795 SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
796 
797 
798    USE control_parameters,                                                    &
799        ONLY:  data_output, message_string
800
801    IMPLICIT NONE
802
803    CHARACTER (LEN=*) ::  unit  !<
804    CHARACTER (LEN=*) ::  var   !<
805
806    INTEGER(iwp) :: i
807    INTEGER(iwp) :: ilen   
808    INTEGER(iwp) :: k
809
810    SELECT CASE ( TRIM( var ) )
811
812       CASE ( 'm_soil' )
813          IF (  .NOT.  land_surface )  THEN
814             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
815                      'res land_surface = .TRUE.'
816             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
817          ENDIF
818          unit = 'm3/m3'
819           
820       CASE ( 't_soil' )
821          IF (  .NOT.  land_surface )  THEN
822             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
823                      'res land_surface = .TRUE.'
824             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
825          ENDIF
826          unit = 'K'   
827             
828       CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'm_liq*',                 &
829              'qsws_liq*', 'qsws_soil*', 'qsws_veg*', 'r_s*' )
830          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
831             message_string = 'illegal value for data_output: "' //            &
832                              TRIM( var ) // '" & only 2d-horizontal ' //      &
833                              'cross sections are allowed for this value'
834             CALL message( 'lsm_check_data_output', 'PA0111', 1, 2, 0, 6, 0 )
835          ENDIF
836          IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
837             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
838                              'res land_surface = .TRUE.'
839             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
840          ENDIF
841          IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
842             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
843                              'res land_surface = .TRUE.'
844             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
845          ENDIF
846          IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
847             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
848                              'res land_surface = .TRUE.'
849             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
850          ENDIF
851          IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
852             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
853                              'res land_surface = .TRUE.'
854             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
855          ENDIF
856          IF ( TRIM( var ) == 'm_liq*'  .AND.  .NOT.  land_surface )  THEN
857             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
858                              'res land_surface = .TRUE.'
859             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
860          ENDIF
861          IF ( TRIM( var ) == 'qsws_liq*'  .AND.  .NOT. land_surface )         &
862          THEN
863             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
864                              'res land_surface = .TRUE.'
865             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
866          ENDIF
867          IF ( TRIM( var ) == 'qsws_soil*'  .AND.  .NOT.  land_surface )       &
868          THEN
869             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
870                              'res land_surface = .TRUE.'
871             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
872          ENDIF
873          IF ( TRIM( var ) == 'qsws_veg*'  .AND.  .NOT. land_surface )         &
874          THEN
875             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
876                              'res land_surface = .TRUE.'
877             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
878          ENDIF
879          IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface )             &
880          THEN
881             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
882                              'res land_surface = .TRUE.'
883             CALL message( 'lsm_check_data_output', 'PA0404', 1, 2, 0, 6, 0 )
884          ENDIF
885
886          IF ( TRIM( var ) == 'lai*'   )      unit = 'none' 
887          IF ( TRIM( var ) == 'c_liq*' )      unit = 'none'
888          IF ( TRIM( var ) == 'c_soil*')      unit = 'none'
889          IF ( TRIM( var ) == 'c_veg*' )      unit = 'none'
890          IF ( TRIM( var ) == 'm_liq*'     )  unit = 'm'
891          IF ( TRIM( var ) == 'qsws_liq*'  )  unit = 'W/m2'
892          IF ( TRIM( var ) == 'qsws_soil*' )  unit = 'W/m2'
893          IF ( TRIM( var ) == 'qsws_veg*'  )  unit = 'W/m2'
894          IF ( TRIM( var ) == 'r_s*')         unit = 's/m' 
895             
896       CASE DEFAULT
897          unit = 'illegal'
898
899    END SELECT
900
901
902 END SUBROUTINE lsm_check_data_output
903
904
905
906!------------------------------------------------------------------------------!
907! Description:
908! ------------
909!> Check data output of profiles for land surface model
910!------------------------------------------------------------------------------!
911 SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
912 
913    USE control_parameters,                                                    &
914        ONLY:  data_output_pr, message_string
915
916    USE indices
917
918    USE profil_parameter
919
920    USE statistics
921
922    IMPLICIT NONE
923   
924    CHARACTER (LEN=*) ::  unit      !<
925    CHARACTER (LEN=*) ::  variable  !<
926    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
927 
928    INTEGER(iwp) ::  var_count     !<
929
930    SELECT CASE ( TRIM( variable ) )
931       
932       CASE ( 't_soil', '#t_soil' )
933          IF (  .NOT.  land_surface )  THEN
934             message_string = 'data_output_pr = ' //                           &
935                              TRIM( data_output_pr(var_count) ) // ' is' //    &
936                              'not implemented for land_surface = .FALSE.'
937             CALL message( 'lsm_check_data_output_pr', 'PA0402', 1, 2, 0, 6, 0 )
938          ELSE
939             dopr_index(var_count) = 89
940             dopr_unit     = 'K'
941             hom(0:nzs-1,2,89,:)  = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
942             IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
943                dopr_initial_index(var_count) = 90
944                hom(0:nzs-1,2,90,:)   = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
945                data_output_pr(var_count)     = data_output_pr(var_count)(2:)
946             ENDIF
947             unit = dopr_unit
948          ENDIF
949
950       CASE ( 'm_soil', '#m_soil' )
951          IF (  .NOT.  land_surface )  THEN
952             message_string = 'data_output_pr = ' //                           &
953                              TRIM( data_output_pr(var_count) ) // ' is' //    &
954                              ' not implemented for land_surface = .FALSE.'
955             CALL message( 'lsm_check_data_output_pr', 'PA0402', 1, 2, 0, 6, 0 )
956          ELSE
957             dopr_index(var_count) = 91
958             dopr_unit     = 'm3/m3'
959             hom(0:nzs-1,2,91,:)  = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
960             IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
961                dopr_initial_index(var_count) = 92
962                hom(0:nzs-1,2,92,:)   = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
963                data_output_pr(var_count)     = data_output_pr(var_count)(2:)
964             ENDIF
965             unit = dopr_unit
966          ENDIF
967
968
969       CASE DEFAULT
970          unit = 'illegal'
971
972    END SELECT
973
974
975 END SUBROUTINE lsm_check_data_output_pr
976 
977 
978!------------------------------------------------------------------------------!
979! Description:
980! ------------
981!> Check parameters routine for land surface model
982!------------------------------------------------------------------------------!
983 SUBROUTINE lsm_check_parameters
984
985    USE control_parameters,                                                    &
986        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string
987                     
988   
989    IMPLICIT NONE
990
991    INTEGER(iwp) ::  i        !< running index, x-dimension
992    INTEGER(iwp) ::  j        !< running index, y-dimension
993    INTEGER(iwp) ::  k        !< running index, z-dimension
994
995!
996!-- Check for a valid setting of surface_type. The default value is 'netcdf'.
997!-- In that case, the surface types are read from NetCDF file
998    IF ( TRIM( surface_type ) /= 'vegetation'  .AND.                           &
999         TRIM( surface_type ) /= 'pavement'    .AND.                           &
1000         TRIM( surface_type ) /= 'water'       .AND.                           &
1001         TRIM( surface_type ) /= 'netcdf' )  THEN 
1002       message_string = 'unknown surface type: surface_type = "' //            &
1003                        TRIM( surface_type ) // '"'
1004       CALL message( 'lsm_check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
1005    ENDIF
1006
1007!
1008!-- Dirichlet boundary conditions are required as the surface fluxes are
1009!-- calculated from the temperature/humidity gradients in the land surface
1010!-- model
1011    IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
1012       message_string = 'lsm requires setting of'//                            &
1013                        'bc_pt_b = "dirichlet" and '//                         &
1014                        'bc_q_b  = "dirichlet"'
1015       CALL message( 'lsm_check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
1016    ENDIF
1017
1018    IF (  .NOT.  constant_flux_layer )  THEN
1019       message_string = 'lsm requires '//                                      &
1020                        'constant_flux_layer = .T.'
1021       CALL message( 'lsm_check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
1022    ENDIF
1023   
1024    IF (  .NOT.  radiation )  THEN
1025       message_string = 'lsm requires '//                                      &
1026                        'the radiation model to be switched on'
1027       CALL message( 'lsm_check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
1028    ENDIF
1029!
1030!-- Check if soil types are set within a valid range.
1031    IF ( TRIM( surface_type ) == 'vegetation'  .OR.                            &
1032         TRIM( surface_type ) == 'pavement' )  THEN
1033       IF ( soil_type < LBOUND( soil_pars, 2 )  .AND.                          &
1034            soil_type > UBOUND( soil_pars, 2 ) )  THEN
1035          WRITE( message_string, * ) 'soil_type = ', soil_type, ' is out ' //  &
1036                                     'of the valid range'
1037          CALL message( 'lsm_check_parameters', 'PA0452', 2, 2, 0, 6, 0 )
1038       ENDIF
1039    ENDIF
1040    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
1041       IF ( soil_type_f%from_file )  THEN
1042          DO  i = nxl, nxr
1043             DO  j = nys, nyn
1044                IF ( soil_type_f%var_2d(j,i) /= soil_type_f%fill  .AND.        &
1045                     ( soil_type_f%var_2d(j,i) < LBOUND( soil_pars, 2 )  .OR.  &
1046                       soil_type_f%var_2d(j,i) > UBOUND( soil_pars, 2 ) ) )  THEN
1047                   WRITE( message_string, * ) 'soil_type = is out  of ' //     &
1048                                        'the valid range at (j,i) = ', j, i
1049                   CALL message( 'lsm_check_parameters', 'PA0452',             &
1050                                  2, 2, myid, 6, 0 )
1051                ENDIF
1052             ENDDO
1053          ENDDO
1054       ENDIF
1055    ENDIF
1056!
1057!-- Check if vegetation types are set within a valid range.   
1058    IF ( TRIM( surface_type ) == 'vegetation' )  THEN
1059       IF ( vegetation_type < LBOUND( vegetation_pars, 2 )  .AND.              &
1060            vegetation_type > UBOUND( vegetation_pars, 2 ) )  THEN
1061          WRITE( message_string, * ) 'vegetation_type = ', vegetation_type,    &
1062                                     ' is out of the valid range'
1063          CALL message( 'lsm_check_parameters', 'PA0526', 2, 2, 0, 6, 0 )
1064       ENDIF
1065    ENDIF
1066    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
1067       IF ( vegetation_type_f%from_file )  THEN
1068          DO  i = nxl, nxr
1069             DO  j = nys, nyn
1070                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill  .AND.&
1071              ( vegetation_type_f%var(j,i) < LBOUND( vegetation_pars, 2 )  .OR.&
1072                vegetation_type_f%var(j,i) > UBOUND( vegetation_pars, 2 ) ) )  &
1073                THEN
1074                   WRITE( message_string, * ) 'vegetation_type = is out of ' //&
1075                                        'the valid range at (j,i) = ', j, i
1076                   CALL message( 'lsm_check_parameters', 'PA0526',             &
1077                                  2, 2, myid, 6, 0 ) 
1078                ENDIF
1079             ENDDO
1080          ENDDO
1081       ENDIF
1082    ENDIF
1083!
1084!-- Check if pavement types are set within a valid range.   
1085    IF ( TRIM( surface_type ) == 'pavement' )  THEN
1086       IF ( pavement_type < LBOUND( pavement_pars, 2 )  .AND.                  &
1087            pavement_type > UBOUND( pavement_pars, 2 ) )  THEN
1088          WRITE( message_string, * ) 'pavement_type = ', pavement_type,        &
1089                                     ' is out of the valid range'
1090          CALL message( 'lsm_check_parameters', 'PA0527', 2, 2, 0, 6, 0 )
1091       ENDIF
1092    ENDIF
1093    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
1094       IF ( pavement_type_f%from_file )  THEN
1095          DO  i = nxl, nxr
1096             DO  j = nys, nyn
1097                IF ( pavement_type_f%var(j,i) /= pavement_type_f%fill  .AND.   &
1098              ( pavement_type_f%var(j,i) < LBOUND( pavement_pars, 2 )  .OR.    &
1099                pavement_type_f%var(j,i) > UBOUND( pavement_pars, 2 ) ) )  THEN
1100                   WRITE( message_string, * ) 'pavement_type = is out of ' //  &
1101                                        'the valid range at (j,i) = ', j, i
1102                   CALL message( 'lsm_check_parameters', 'PA0527',             &
1103                                  2, 2, myid, 6, 0 ) 
1104                ENDIF
1105             ENDDO
1106          ENDDO
1107       ENDIF
1108    ENDIF
1109!
1110!-- Check if water types are set within a valid range.   
1111    IF ( TRIM( surface_type ) == 'water' )  THEN
1112       IF ( water_type < LBOUND( water_pars, 2 )  .AND.                        &
1113            water_type > UBOUND( water_pars, 2 ) )  THEN
1114          WRITE( message_string, * ) 'water_type = ', water_type,              &
1115                                     ' is out of the valid range'
1116          CALL message( 'lsm_check_parameters', 'PA0528', 2, 2, 0, 6, 0 )
1117       ENDIF
1118    ENDIF
1119    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
1120       IF ( water_type_f%from_file )  THEN
1121          DO  i = nxl, nxr
1122             DO  j = nys, nyn
1123                IF ( water_type_f%var(j,i) /= water_type_f%fill  .AND.         &
1124              ( water_type_f%var(j,i) < LBOUND( water_pars, 2 )  .OR.          &
1125                water_type_f%var(j,i) > UBOUND( water_pars, 2 ) ) )  THEN
1126                   WRITE( message_string, * ) 'water_type = is out  of ' //    &
1127                                        'the valid range at (j,i) = ', j, i
1128                   CALL message( 'lsm_check_parameters', 'PA0528',             &
1129                                 2, 2, myid, 6, 0 ) 
1130                ENDIF
1131             ENDDO
1132          ENDDO
1133       ENDIF
1134    ENDIF
1135!
1136!-- Check further settings for consistency.
1137    IF ( TRIM( surface_type ) == 'vegetation' )  THEN
1138   
1139       IF ( vegetation_type == 0 )  THEN
1140          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
1141             message_string = 'vegetation_type = 0 (user defined)'//           &
1142                              'requires setting of min_canopy_resistance'//    &
1143                              '/= 9999999.9'
1144             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1145          ENDIF
1146
1147          IF ( leaf_area_index == 9999999.9_wp )  THEN
1148             message_string = 'vegetation_type = 0 (user_defined)'//           &
1149                              'requires setting of leaf_area_index'//          &
1150                              '/= 9999999.9'
1151             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1152          ENDIF
1153
1154          IF ( vegetation_coverage == 9999999.9_wp )  THEN
1155             message_string = 'vegetation_type = 0 (user_defined)'//           &
1156                              'requires setting of vegetation_coverage'//      &
1157                              '/= 9999999.9'
1158             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1159          ENDIF
1160
1161          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
1162             message_string = 'vegetation_type = 0 (user_defined)'//           &
1163                              'requires setting of'//                          &
1164                              'canopy_resistance_coefficient /= 9999999.9'
1165             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1166          ENDIF
1167
1168          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
1169             message_string = 'vegetation_type = 0 (user_defined)'//           &
1170                              'requires setting of lambda_surface_stable'//    &
1171                              '/= 9999999.9'
1172             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1173          ENDIF
1174
1175          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
1176             message_string = 'vegetation_type = 0 (user_defined)'//           &
1177                              'requires setting of lambda_surface_unstable'//  &
1178                              '/= 9999999.9'
1179             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1180          ENDIF
1181
1182          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
1183             message_string = 'vegetation_type = 0 (user_defined)'//           &
1184                              'requires setting of f_shortwave_incoming'//     &
1185                              '/= 9999999.9'
1186             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1187          ENDIF
1188
1189          IF ( z0_vegetation == 9999999.9_wp )  THEN
1190             message_string = 'vegetation_type = 0 (user_defined)'//           &
1191                              'requires setting of z0_vegetation'//            &
1192                              '/= 9999999.9'
1193             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1194          ENDIF
1195
1196          IF ( z0h_vegetation == 9999999.9_wp )  THEN
1197             message_string = 'vegetation_type = 0 (user_defined)'//           &
1198                              'requires setting of z0h_vegetation'//           &
1199                              '/= 9999999.9'
1200             CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1201          ENDIF
1202       ENDIF
1203
1204       IF ( vegetation_type == 1 )  THEN
1205          IF ( vegetation_coverage /= 9999999.9_wp  .AND.  vegetation_coverage &
1206               /= 0.0_wp )  THEN
1207             message_string = 'vegetation_type = 1 (bare soil)'//              &
1208                              ' requires vegetation_coverage = 0'
1209             CALL message( 'lsm_check_parameters', 'PA0294', 1, 2, 0, 6, 0 )
1210          ENDIF
1211       ENDIF
1212 
1213    ENDIF
1214   
1215    IF ( TRIM( surface_type ) == 'water' )  THEN
1216
1217       IF ( water_type == 0 )  THEN 
1218       
1219          IF ( z0_water == 9999999.9_wp )  THEN
1220             message_string = 'water_type = 0 (user_defined)'//                &
1221                              'requires setting of z0_water'//                 &
1222                              '/= 9999999.9'
1223             CALL message( 'lsm_check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
1224          ENDIF
1225
1226          IF ( z0h_water == 9999999.9_wp )  THEN
1227             message_string = 'water_type = 0 (user_defined)'//                &
1228                              'requires setting of z0h_water'//                &
1229                              '/= 9999999.9'
1230             CALL message( 'lsm_check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
1231          ENDIF
1232         
1233          IF ( water_temperature == 9999999.9_wp )  THEN
1234             message_string = 'water_type = 0 (user_defined)'//                &
1235                              'requires setting of water_temperature'//        &
1236                              '/= 9999999.9'
1237             CALL message( 'lsm_check_parameters', 'PA0379', 1, 2, 0, 6, 0 )
1238          ENDIF       
1239         
1240       ENDIF
1241       
1242    ENDIF
1243   
1244    IF ( TRIM( surface_type ) == 'pavement' )  THEN
1245
1246       IF ( ANY( dz_soil /= 9999999.9_wp )  .AND.  pavement_type /= 0 )  THEN
1247          message_string = 'non-default setting of dz_soil '//                  &
1248                           'does not allow to use pavement_type /= 0)'
1249          CALL message( 'lsm_check_parameters', 'PA0341', 1, 2, 0, 6, 0 )
1250       ENDIF
1251
1252       IF ( pavement_type == 0 )  THEN 
1253       
1254          IF ( z0_pavement == 9999999.9_wp )  THEN
1255             message_string = 'pavement_type = 0 (user_defined)'//             &
1256                              'requires setting of z0_pavement'//              &
1257                              '/= 9999999.9'
1258             CALL message( 'lsm_check_parameters', 'PA0352', 1, 2, 0, 6, 0 )
1259          ENDIF
1260         
1261          IF ( z0h_pavement == 9999999.9_wp )  THEN
1262             message_string = 'pavement_type = 0 (user_defined)'//             &
1263                              'requires setting of z0h_pavement'//             &
1264                              '/= 9999999.9'
1265             CALL message( 'lsm_check_parameters', 'PA0353', 1, 2, 0, 6, 0 )
1266          ENDIF
1267         
1268          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
1269             message_string = 'pavement_type = 0 (user_defined)'//             &
1270                              'requires setting of pavement_heat_conduct'//    &
1271                              '/= 9999999.9'
1272             CALL message( 'lsm_check_parameters', 'PA0342', 1, 2, 0, 6, 0 )
1273          ENDIF
1274         
1275           IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
1276             message_string = 'pavement_type = 0 (user_defined)'//             &
1277                              'requires setting of pavement_heat_capacity'//   &
1278                              '/= 9999999.9'
1279             CALL message( 'lsm_check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
1280          ENDIF
1281
1282          IF ( pavement_depth_level == 0 )  THEN
1283             message_string = 'pavement_type = 0 (user_defined)'//             &
1284                              'requires setting of pavement_depth_level'//     &
1285                              '/= 0'
1286             CALL message( 'lsm_check_parameters', 'PA0474', 1, 2, 0, 6, 0 )
1287          ENDIF
1288
1289       ENDIF
1290   
1291    ENDIF
1292
1293    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
1294!
1295!--    MS: Some problme here, after calling message everythings stucks at
1296!--        MPI_FINALIZE call.
1297       IF ( ANY( pavement_type_f%var /= pavement_type_f%fill )  .AND.           &
1298            ANY( dz_soil /= 9999999.9_wp ) )  THEN
1299          message_string = 'pavement-surfaces are not allowed in ' //           &
1300                           'combination with a non-default setting of dz_soil'
1301          CALL message( 'lsm_check_parameters', 'PA0316', 2, 2, 0, 6, 0 )
1302       ENDIF
1303    ENDIF
1304   
1305!
1306!-- Temporary message as long as NetCDF input is not available
1307    IF ( TRIM( surface_type ) == 'netcdf'  .AND.  .NOT. input_pids_static )   &
1308    THEN
1309       message_string = 'surface_type = netcdf requires static input file.'
1310       CALL message( 'lsm_check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
1311    ENDIF
1312
1313    IF ( soil_type == 0  .AND.  .NOT. input_pids_static )  THEN
1314
1315       IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
1316          message_string = 'soil_type = 0 (user_defined)'//                    &
1317                           'requires setting of alpha_vangenuchten'//          &
1318                           '/= 9999999.9'
1319          CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1320       ENDIF
1321
1322       IF ( l_vangenuchten == 9999999.9_wp )  THEN
1323          message_string = 'soil_type = 0 (user_defined)'//                    &
1324                           'requires setting of l_vangenuchten'//              &
1325                           '/= 9999999.9'
1326          CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1327       ENDIF
1328
1329       IF ( n_vangenuchten == 9999999.9_wp )  THEN
1330          message_string = 'soil_type = 0 (user_defined)'//                    &
1331                           'requires setting of n_vangenuchten'//              &
1332                           '/= 9999999.9'
1333          CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1334       ENDIF
1335
1336       IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
1337          message_string = 'soil_type = 0 (user_defined)'//                    &
1338                           'requires setting of hydraulic_conductivity'//      &
1339                           '/= 9999999.9'
1340          CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1341       ENDIF
1342
1343       IF ( saturation_moisture == 9999999.9_wp )  THEN
1344          message_string = 'soil_type = 0 (user_defined)'//                    &
1345                           'requires setting of saturation_moisture'//         &
1346                           '/= 9999999.9'
1347          CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1348       ENDIF
1349
1350       IF ( field_capacity == 9999999.9_wp )  THEN
1351          message_string = 'soil_type = 0 (user_defined)'//                    &
1352                           'requires setting of field_capacity'//              &
1353                           '/= 9999999.9'
1354          CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1355       ENDIF
1356
1357       IF ( wilting_point == 9999999.9_wp )  THEN
1358          message_string = 'soil_type = 0 (user_defined)'//                    &
1359                           'requires setting of wilting_point'//               &
1360                           '/= 9999999.9'
1361          CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1362       ENDIF
1363
1364       IF ( residual_moisture == 9999999.9_wp )  THEN
1365          message_string = 'soil_type = 0 (user_defined)'//                    &
1366                           'requires setting of residual_moisture'//           &
1367                           '/= 9999999.9'
1368          CALL message( 'lsm_check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
1369       ENDIF
1370
1371    ENDIF
1372
1373
1374!!! these checks are not needed for water surfaces??
1375
1376!
1377!-- Determine number of soil layers to be used and check whether an appropriate
1378!-- root fraction is prescribed
1379    nzb_soil = 0
1380    nzt_soil = -1
1381    IF ( ALL( dz_soil == 9999999.9_wp ) )  THEN
1382       nzt_soil = 7
1383       dz_soil(nzb_soil:nzt_soil) = dz_soil_default
1384    ELSE
1385       DO k = 0, 19
1386          IF ( dz_soil(k) /= 9999999.9_wp )  THEN
1387             nzt_soil = nzt_soil + 1
1388          ENDIF
1389       ENDDO   
1390    ENDIF
1391    nzs = nzt_soil + 1
1392
1393!
1394!-- Check whether valid soil temperatures are prescribed
1395    IF ( .NOT. input_pids_dynamic )  THEN
1396       IF ( COUNT( soil_temperature /= 9999999.9_wp ) /= nzs )  THEN
1397          WRITE( message_string, * )                                           &
1398                                  'number of soil layers (', nzs, ') does not',&
1399                                  ' match to the number of layers specified',  &
1400                                  ' in soil_temperature (', COUNT(             &
1401                                  soil_temperature /= 9999999.9_wp ), ')'
1402             CALL message( 'lsm_check_parameters', 'PA0471', 1, 2, 0, 6, 0 )
1403       ENDIF
1404
1405       IF ( deep_soil_temperature == 9999999.9_wp ) THEN
1406             message_string = 'deep_soil_temperature is not set but must be'// &
1407                              '/= 9999999.9'
1408             CALL message( 'lsm_check_parameters', 'PA0472', 1, 2, 0, 6, 0 )
1409       ENDIF
1410    ENDIF
1411
1412!
1413!-- Check whether the sum of all root fractions equals one
1414    IF ( vegetation_type == 0 )  THEN
1415       IF ( SUM( root_fraction(nzb_soil:nzt_soil) ) /= 1.0_wp )  THEN
1416          message_string = 'vegetation_type = 0 (user_defined)'//              &
1417                           'requires setting of root_fraction'//               &
1418                           '/= 9999999.9 and SUM(root_fraction) = 1'
1419          CALL message( 'lsm_check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
1420       ENDIF
1421    ENDIF     
1422!
1423!-- Calculate grid spacings. Temperature and moisture are defined at
1424!-- the center of the soil layers, whereas gradients/fluxes are
1425!-- defined at the edges (_layer)
1426!
1427!-- Allocate global 1D arrays
1428    ALLOCATE ( ddz_soil_center(nzb_soil:nzt_soil) )
1429    ALLOCATE ( ddz_soil(nzb_soil:nzt_soil+1) )
1430    ALLOCATE ( dz_soil_center(nzb_soil:nzt_soil) )
1431    ALLOCATE ( zs(nzb_soil:nzt_soil+1) )
1432
1433
1434    zs(nzb_soil) = 0.5_wp * dz_soil(nzb_soil)
1435    zs_layer(nzb_soil) = dz_soil(nzb_soil)
1436
1437    DO  k = nzb_soil+1, nzt_soil
1438       zs_layer(k) = zs_layer(k-1) + dz_soil(k)
1439       zs(k) = (zs_layer(k) +  zs_layer(k-1)) * 0.5_wp
1440    ENDDO
1441
1442    dz_soil(nzt_soil+1) = zs_layer(nzt_soil) + dz_soil(nzt_soil)
1443    zs(nzt_soil+1) = zs_layer(nzt_soil) + 0.5_wp * dz_soil(nzt_soil)
1444 
1445    DO  k = nzb_soil, nzt_soil-1
1446       dz_soil_center(k) = zs(k+1) - zs(k)
1447       IF ( dz_soil_center(k) <= 0.0_wp )  THEN
1448          message_string = 'invalid soil layer configuration found ' //        &
1449                           '(dz_soil_center(k) <= 0.0)'
1450          CALL message( 'lsm_rrd_local', 'PA0140', 1, 2, 0, 6, 0 )
1451       ENDIF 
1452    ENDDO
1453 
1454    dz_soil_center(nzt_soil) = zs_layer(k-1) + dz_soil(k) - zs(nzt_soil)
1455       
1456    ddz_soil_center = 1.0_wp / dz_soil_center
1457    ddz_soil(nzb_soil:nzt_soil) = 1.0_wp / dz_soil(nzb_soil:nzt_soil)
1458
1459
1460
1461 END SUBROUTINE lsm_check_parameters
1462 
1463!------------------------------------------------------------------------------!
1464! Description:
1465! ------------
1466!> Solver for the energy balance at the surface.
1467!------------------------------------------------------------------------------!
1468 SUBROUTINE lsm_energy_balance( horizontal, l )
1469
1470    USE pegrid
1471    USE radiation_model_mod,  ONLY:  rad_lw_out
1472
1473    IMPLICIT NONE
1474
1475    INTEGER(iwp) ::  i         !< running index
1476    INTEGER(iwp) ::  i_off     !< offset to determine index of surface element, seen from atmospheric grid point, for x
1477    INTEGER(iwp) ::  j         !< running index
1478    INTEGER(iwp) ::  j_off     !< offset to determine index of surface element, seen from atmospheric grid point, for y
1479    INTEGER(iwp) ::  k         !< running index
1480    INTEGER(iwp) ::  k_off     !< offset to determine index of surface element, seen from atmospheric grid point, for z
1481    INTEGER(iwp) ::  ks        !< running index
1482    INTEGER(iwp) ::  l         !< surface-facing index
1483    INTEGER(iwp) ::  m         !< running index concerning wall elements
1484
1485    LOGICAL      ::  horizontal !< Flag indicating horizontal or vertical surfaces
1486
1487    REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
1488                f1,          & !< resistance correction term 1
1489                f2,          & !< resistance correction term 2
1490                f3,          & !< resistance correction term 3
1491                m_min,       & !< minimum soil moisture
1492                e,           & !< water vapour pressure
1493                e_s,         & !< water vapour saturation pressure
1494                e_s_dt,      & !< derivate of e_s with respect to T
1495                tend,        & !< tendency
1496                dq_s_dt,     & !< derivate of q_s with respect to T
1497                coef_1,      & !< coef. for prognostic equation
1498                coef_2,      & !< coef. for prognostic equation
1499                f_qsws,      & !< factor for qsws
1500                f_qsws_veg,  & !< factor for qsws_veg
1501                f_qsws_soil, & !< factor for qsws_soil
1502                f_qsws_liq,  & !< factor for qsws_liq
1503                f_shf,       & !< factor for shf
1504                lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
1505                lambda_surface, & !< Current value of lambda_surface (W/m2/K)
1506                m_liq_max      !< maxmimum value of the liq. water reservoir
1507
1508    TYPE(surf_type_lsm), POINTER ::  surf_t_surface
1509    TYPE(surf_type_lsm), POINTER ::  surf_t_surface_p
1510    TYPE(surf_type_lsm), POINTER ::  surf_tt_surface_m
1511    TYPE(surf_type_lsm), POINTER ::  surf_m_liq
1512    TYPE(surf_type_lsm), POINTER ::  surf_m_liq_p
1513    TYPE(surf_type_lsm), POINTER ::  surf_tm_liq_m
1514
1515    TYPE(surf_type_lsm), POINTER ::  surf_m_soil
1516    TYPE(surf_type_lsm), POINTER ::  surf_t_soil
1517
1518    TYPE(surf_type), POINTER  ::  surf  !< surface-date type variable
1519
1520
1521    IF ( debug_output_timestep )  THEN
1522       WRITE( debug_string, * ) 'lsm_energy_balance', horizontal, l
1523       CALL debug_message( debug_string, 'start' )
1524    ENDIF
1525
1526    IF ( horizontal )  THEN
1527       surf              => surf_lsm_h
1528
1529       surf_t_surface    => t_surface_h
1530       surf_t_surface_p  => t_surface_h_p
1531       surf_tt_surface_m => tt_surface_h_m
1532       surf_m_liq        => m_liq_h
1533       surf_m_liq_p      => m_liq_h_p
1534       surf_tm_liq_m     => tm_liq_h_m
1535       surf_m_soil       => m_soil_h
1536       surf_t_soil       => t_soil_h
1537    ELSE
1538       surf              => surf_lsm_v(l)
1539
1540       surf_t_surface    => t_surface_v(l)
1541       surf_t_surface_p  => t_surface_v_p(l)
1542       surf_tt_surface_m => tt_surface_v_m(l)
1543       surf_m_liq        => m_liq_v(l)
1544       surf_m_liq_p      => m_liq_v_p(l)
1545       surf_tm_liq_m     => tm_liq_v_m(l)
1546       surf_m_soil       => m_soil_v(l)
1547       surf_t_soil       => t_soil_v(l)
1548    ENDIF
1549
1550!
1551!-- Index offset of surface element point with respect to adjoining
1552!-- atmospheric grid point
1553    k_off = surf%koff
1554    j_off = surf%joff
1555    i_off = surf%ioff
1556
1557    !$OMP PARALLEL PRIVATE (m, i, j, k, lambda_h_sat, ke, lambda_soil, lambda_surface,             &
1558    !$OMP&                  c_surface_tmp, f1,m_total, f2, e_s, e, f3, m_min, m_liq_max, q_s,      &
1559    !$OMP&                  f_qsws_veg, f_qsws_soil, f_qsws_liq, f_shf, f_qsws, e_s_dt, dq_s_dt,   &
1560    !$OMP&                  coef_1, coef_2, tend)
1561    !$OMP DO SCHEDULE (STATIC)
1562    DO  m = 1, surf%ns
1563
1564       i   = surf%i(m)           
1565       j   = surf%j(m)
1566       k   = surf%k(m)
1567
1568!
1569!--    Define heat conductivity between surface and soil depending on surface
1570!--    type. For vegetation, a skin layer parameterization is used. The new
1571!--    parameterization uses a combination of two conductivities: a constant
1572!--    conductivity for the skin layer, and a conductivity according to the
1573!--    uppermost soil layer. For bare soil and pavements, no skin layer is
1574!--    applied. In these cases, the temperature is assumed to be constant
1575!--    between the surface and the first soil layer. The heat conductivity is
1576!--    then derived from the soil/pavement properties.
1577!--    For water surfaces, the conductivity is already set to 1E10.
1578!--    Moreover, the heat capacity is set. For bare soil the heat capacity is
1579!--    the capacity of the uppermost soil layer, for pavement it is that of
1580!--    the material involved.
1581
1582!
1583!--    for vegetation type surfaces, the thermal conductivity of the soil is
1584!--    needed
1585
1586       IF ( surf%vegetation_surface(m) )  THEN
1587
1588          lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(nzb_soil,m)) *      &
1589                         lambda_h_water ** surf_m_soil%var_2d(nzb_soil,m)
1590                         
1591          ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(nzb_soil,m) /   &
1592                                                     surf%m_sat(nzb_soil,m) ) )                   
1593                         
1594          lambda_soil = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry )   &
1595                           * ddz_soil(nzb_soil) * 2.0_wp
1596
1597!
1598!--       When bare soil is set without a thermal conductivity (no skin layer),
1599!--       a heat capacity is that of the soil layer, otherwise it is a
1600!--       combination of the conductivities from the skin and the soil layer
1601          IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN
1602            surf%c_surface(m) = (rho_c_soil * (1.0_wp - surf%m_sat(nzb_soil,m))&
1603                              + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) ) &
1604                              * dz_soil(nzb_soil) * 0.5_wp   
1605            lambda_surface = lambda_soil
1606
1607          ELSE IF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m))&
1608          THEN
1609             lambda_surface = surf%lambda_surface_s(m) * lambda_soil           &
1610                              / ( surf%lambda_surface_s(m) + lambda_soil )
1611          ELSE
1612
1613             lambda_surface = surf%lambda_surface_u(m) * lambda_soil           &
1614                              / ( surf%lambda_surface_u(m) + lambda_soil )
1615          ENDIF
1616       ELSE
1617          lambda_surface = surf%lambda_surface_s(m)
1618       ENDIF
1619
1620!
1621!--    Set heat capacity of the skin/surface. It is ususally zero when a skin
1622!--    layer is used, and non-zero otherwise.
1623       c_surface_tmp = surf%c_surface(m) 
1624
1625!
1626!--    First step: calculate aerodyamic resistance. As pt, us, ts
1627!--    are not available for the prognostic time step, data from the last
1628!--    time step is used here. Note that this formulation is the
1629!--    equivalent to the ECMWF formulation using drag coefficients
1630!        IF ( bulk_cloud_model )  THEN
1631!           pt1 = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
1632!           qv1 = q(k,j,i) - ql(k,j,i)
1633!        ELSEIF ( cloud_droplets ) THEN
1634!           pt1 = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
1635!           qv1 = q(k,j,i)
1636!        ELSE
1637!           pt1 = pt(k,j,i)
1638!           IF ( humidity )  THEN
1639!              qv1 = q(k,j,i)
1640!           ELSE
1641!              qv1 = 0.0_wp
1642!           ENDIF
1643!        ENDIF
1644!
1645!--     Calculation of r_a for vertical surfaces
1646!--
1647!--     heat transfer coefficient for forced convection along vertical walls
1648!--     follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
1649!--           
1650!--       H = httc (Tsfc - Tair)
1651!--       httc = rw * (11.8 + 4.2 * Ueff) - 4.0
1652!--           
1653!--             rw: wall patch roughness relative to 1.0 for concrete
1654!--             Ueff: effective wind speed
1655!--             - 4.0 is a reduction of Rowley et al (1930) formulation based on
1656!--             Cole and Sturrock (1977)
1657!--           
1658!--             Ucan: Canyon wind speed
1659!--             wstar: convective velocity
1660!--             Qs: surface heat flux
1661!--             zH: height of the convective layer
1662!--             wstar = (g/Tcan*Qs*zH)**(1./3.)
1663               
1664!--    Effective velocity components must always
1665!--    be defined at scalar grid point. The wall normal component is
1666!--    obtained by simple linear interpolation. ( An alternative would
1667!--    be an logarithmic interpolation. )
1668!--    A roughness lenght of 0.001 is assumed for concrete (the inverse,
1669!--    1000 is used in the nominator for scaling)
1670!--    To do: detailed investigation which approach gives more reliable results!
1671!--    Please note, in case of very small friction velocity, e.g. in little
1672!--    holes, the resistance can become negative. For this reason, limit r_a
1673!--    to positive values.
1674       IF ( horizontal  .OR.  .NOT. aero_resist_kray )  THEN
1675          surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) /            &
1676                             ( surf%ts(m) * surf%us(m) + 1.0E-20_wp ) )
1677       ELSE
1678          surf%r_a(m) = rho_cp / ( surf%z0(m) * 1000.0_wp                      &
1679                        * ( 11.8_wp + 4.2_wp *                                 &
1680                        SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
1681                                   ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
1682                                   ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2,  &
1683                              0.01_wp ) )                                      &
1684                           )  - 4.0_wp  ) 
1685       ENDIF
1686!
1687!--    Make sure that the resistance does not drop to zero for neutral
1688!--    stratification. Also, set a maximum resistance to avoid the breakdown of
1689!--    MOST for locations with zero wind speed
1690       IF ( surf%r_a(m) <   1.0_wp )  surf%r_a(m) =   1.0_wp
1691       IF ( surf%r_a(m) > 300.0_wp )  surf%r_a(m) = 300.0_wp       
1692!
1693!--    Second step: calculate canopy resistance r_canopy
1694!--    f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
1695 
1696!--    f1: correction for incoming shortwave radiation (stomata close at
1697!--    night)
1698       f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) /          &
1699                        (0.81_wp * (0.004_wp * surf%rad_sw_in(m)               &
1700                         + 1.0_wp)) )
1701
1702!
1703!--    f2: correction for soil moisture availability to plants (the
1704!--    integrated soil moisture must thus be considered here)
1705!--    f2 = 0 for very dry soils
1706       m_total = 0.0_wp
1707       DO  ks = nzb_soil, nzt_soil
1708           m_total = m_total + surf%root_fr(ks,m)                              &
1709                     * MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(ks,m) )
1710       ENDDO 
1711
1712!
1713!--    The calculation of f2 is based on only one wilting point value for all
1714!--    soil layers. The value at k=nzb_soil is used here as a proxy but might
1715!--    need refinement in the future.
1716       IF ( m_total > surf%m_wilt(nzb_soil,m)  .AND.                           &
1717            m_total < surf%m_fc(nzb_soil,m) )  THEN
1718          f2 = ( m_total - surf%m_wilt(nzb_soil,m) ) /                         &
1719               ( surf%m_fc(nzb_soil,m) - surf%m_wilt(nzb_soil,m) )
1720       ELSEIF ( m_total >= surf%m_fc(nzb_soil,m) )  THEN
1721          f2 = 1.0_wp
1722       ELSE
1723          f2 = 1.0E-20_wp
1724       ENDIF
1725
1726!
1727!--    Calculate water vapour pressure at saturation and convert to hPa
1728!--    The magnus formula is limited to temperatures up to 333.15 K to
1729!--    avoid negative values of q_s
1730       e_s = 0.01_wp * magnus( MIN(surf_t_surface%var_1d(m), 333.15_wp) )
1731
1732!
1733!--    f3: correction for vapour pressure deficit
1734       IF ( surf%g_d(m) /= 0.0_wp )  THEN
1735!
1736!--       Calculate vapour pressure
1737          e  = surf%qv1(m) * surface_pressure / ( surf%qv1(m) + rd_d_rv )
1738          f3 = EXP ( - surf%g_d(m) * (e_s - e) )
1739       ELSE
1740          f3 = 1.0_wp
1741       ENDIF
1742!
1743!--    Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1744!--    this calculation is obsolete, as r_canopy is not used below.
1745!--    To do: check for very dry soil -> r_canopy goes to infinity
1746       surf%r_canopy(m) = surf%r_canopy_min(m) /                               &
1747                              ( surf%lai(m) * f1 * f2 * f3 + 1.0E-20_wp )
1748!
1749!--    Third step: calculate bare soil resistance r_soil.
1750       m_min = surf%c_veg(m) * surf%m_wilt(nzb_soil,m) +                       &
1751                         ( 1.0_wp - surf%c_veg(m) ) * surf%m_res(nzb_soil,m)
1752
1753
1754       f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) /                       &
1755            ( surf%m_fc(nzb_soil,m) - m_min )
1756       f2 = MAX( f2, 1.0E-20_wp )
1757       f2 = MIN( f2, 1.0_wp     )
1758
1759       surf%r_soil(m) = surf%r_soil_min(m) / f2
1760       
1761!
1762!--    Calculate the maximum possible liquid water amount on plants and
1763!--    bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
1764!--    assumed, while paved surfaces might hold up 1 mm of water. The
1765!--    liquid water fraction for paved surfaces is calculated after
1766!--    Noilhan & Planton (1989), while the ECMWF formulation is used for
1767!--    vegetated surfaces and bare soils.
1768       IF ( surf%pavement_surface(m) )  THEN
1769          m_liq_max = m_max_depth * 5.0_wp
1770          surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq%var_1d(m) / m_liq_max)**0.67 )
1771       ELSE
1772          m_liq_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m)              &
1773                      + ( 1.0_wp - surf%c_veg(m) ) )
1774          surf%c_liq(m) = MIN( 1.0_wp, surf_m_liq%var_1d(m) / m_liq_max )
1775       ENDIF
1776!
1777!--    Calculate saturation water vapor mixing ratio
1778       q_s = rd_d_rv * e_s / ( surface_pressure - e_s )
1779!
1780!--    In case of dewfall, set evapotranspiration to zero
1781!--    All super-saturated water is then removed from the air
1782       IF ( humidity  .AND.  q_s <= surf%qv1(m) )  THEN
1783          surf%r_canopy(m) = 0.0_wp
1784          surf%r_soil(m)   = 0.0_wp
1785       ENDIF
1786
1787!
1788!--    Calculate coefficients for the total evapotranspiration
1789!--    In case of water surface, set vegetation and soil fluxes to zero.
1790!--    For pavements, only evaporation of liquid water is possible.
1791       IF ( surf%water_surface(m) )  THEN
1792          f_qsws_veg  = 0.0_wp
1793          f_qsws_soil = 0.0_wp
1794          f_qsws_liq  = rho_lv / surf%r_a(m)
1795       ELSEIF ( surf%pavement_surface(m) )  THEN
1796          f_qsws_veg  = 0.0_wp
1797          f_qsws_soil = 0.0_wp
1798          f_qsws_liq  = rho_lv * surf%c_liq(m) / surf%r_a(m)
1799       ELSE
1800          f_qsws_veg  = rho_lv * surf%c_veg(m) *                               &
1801                            ( 1.0_wp        - surf%c_liq(m)    ) /             &
1802                            ( surf%r_a(m) + surf%r_canopy(m) )
1803          f_qsws_soil = rho_lv * (1.0_wp    - surf%c_veg(m)    ) /             &
1804                            ( surf%r_a(m) + surf%r_soil(m)   )
1805          f_qsws_liq  = rho_lv * surf%c_veg(m) * surf%c_liq(m)   /             &
1806                              surf%r_a(m)
1807       ENDIF
1808
1809       f_shf  = rho_cp / surf%r_a(m)
1810       f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1811!
1812!--    Calculate derivative of q_s for Taylor series expansion
1813       e_s_dt = e_s * ( 17.62_wp / ( surf_t_surface%var_1d(m) - 29.65_wp) -   &
1814                        17.62_wp*( surf_t_surface%var_1d(m) - 273.15_wp)      &
1815                       / ( surf_t_surface%var_1d(m) - 29.65_wp)**2 )
1816
1817       dq_s_dt = rd_d_rv * e_s_dt / ( surface_pressure - e_s_dt )
1818!
1819!--    Calculate net radiation radiation without longwave outgoing flux because
1820!--    it has a dependency on surface temperature and thus enters the prognostic
1821!--    equations directly
1822       surf%rad_net_l(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)              &
1823                           + surf%rad_lw_in(m)
1824!
1825!--    Calculate new skin temperature
1826       IF ( humidity )  THEN
1827!
1828!--       Numerator of the prognostic equation
1829          coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m)             &
1830                   * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)             &
1831                   + f_shf * surf%pt1(m) + f_qsws * ( surf%qv1(m) - q_s        &
1832                   + dq_s_dt * surf_t_surface%var_1d(m) ) + lambda_surface     &
1833                   * surf_t_soil%var_2d(nzb_soil,m)
1834
1835!
1836!--       Denominator of the prognostic equation
1837          coef_2 = surf%rad_lw_out_change_0(m) + f_qsws * dq_s_dt              &
1838                   + lambda_surface + f_shf / exner(nzb)
1839       ELSE
1840!
1841!--       Numerator of the prognostic equation
1842          coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m)             &
1843                   * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)             &
1844                   + f_shf * surf%pt1(m)  + lambda_surface                     &
1845                   * surf_t_soil%var_2d(nzb_soil,m)
1846!
1847!--       Denominator of the prognostic equation
1848          coef_2 = surf%rad_lw_out_change_0(m) + lambda_surface + f_shf / exner(nzb)
1849
1850       ENDIF
1851
1852       tend = 0.0_wp
1853
1854!
1855!--    Implicit solution when the surface layer has no heat capacity,
1856!--    otherwise use RK3 scheme.
1857       surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *&
1858                          surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2&
1859                                             * dt_3d * tsc(2) ) 
1860
1861!
1862!--    Add RK3 term
1863       IF ( c_surface_tmp /= 0.0_wp )  THEN
1864
1865          surf_t_surface_p%var_1d(m) = surf_t_surface_p%var_1d(m) + dt_3d *    &
1866                                       tsc(3) * surf_tt_surface_m%var_1d(m)
1867
1868!
1869!--       Calculate true tendency
1870          tend = ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) -     &
1871                   dt_3d * tsc(3) * surf_tt_surface_m%var_1d(m)) / (dt_3d  * tsc(2))
1872!
1873!--       Calculate t_surface tendencies for the next Runge-Kutta step
1874          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1875             IF ( intermediate_timestep_count == 1 )  THEN
1876                surf_tt_surface_m%var_1d(m) = tend
1877             ELSEIF ( intermediate_timestep_count <                            &
1878                      intermediate_timestep_count_max )  THEN
1879                surf_tt_surface_m%var_1d(m) = -9.5625_wp * tend +              &
1880                                               5.3125_wp * surf_tt_surface_m%var_1d(m)
1881             ENDIF
1882          ENDIF
1883       ENDIF
1884
1885!
1886!--    In case of fast changes in the skin temperature, it is possible to
1887!--    update the radiative fluxes independently from the prescribed
1888!--    radiation call frequency. This effectively prevents oscillations,
1889!--    especially when setting skip_time_do_radiation /= 0. The threshold
1890!--    value of 0.2 used here is just a first guess. This method should be
1891!--    revised in the future as tests have shown that the threshold is
1892!--    often reached, when no oscillations would occur (causes immense
1893!--    computing time for the radiation code).
1894       IF ( ABS( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )       &
1895            > 0.2_wp  .AND. &
1896            unscheduled_radiation_calls )  THEN
1897          force_radiation_call_l = .TRUE.
1898       ENDIF
1899
1900       surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exner(nzb)
1901
1902!
1903!--    Calculate fluxes
1904       surf%rad_net_l(m) = surf%rad_net_l(m) +                                 &
1905                            surf%rad_lw_out_change_0(m)                        &
1906                          * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)      &
1907                          - surf%rad_lw_out_change_0(m) * surf_t_surface_p%var_1d(m)
1908
1909       surf%rad_net(m) = surf%rad_net_l(m)
1910       surf%rad_lw_out(m) = surf%rad_lw_out(m) + surf%rad_lw_out_change_0(m) * &
1911                     ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )
1912
1913       surf%ghf(m) = lambda_surface * ( surf_t_surface_p%var_1d(m)             &
1914                                             - surf_t_soil%var_2d(nzb_soil,m) )
1915
1916       surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / c_p
1917
1918!
1919! update the 3d field of rad_lw_out array to have consistent output
1920       IF ( horizontal ) THEN
1921          IF ( radiation_scheme == 'rrtmg' ) THEN
1922             rad_lw_out(k+k_off,j+j_off,i+i_off) = surf%rad_lw_out(m)
1923          ELSE
1924             rad_lw_out(0,j+j_off,i+i_off) = surf%rad_lw_out(m)
1925          ENDIF
1926       ENDIF
1927
1928       IF ( humidity )  THEN
1929          surf%qsws(m)  = - f_qsws * ( surf%qv1(m) - q_s + dq_s_dt             &
1930                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
1931                            surf_t_surface_p%var_1d(m) )
1932
1933          surf%qsws_veg(m)  = - f_qsws_veg  * ( surf%qv1(m) - q_s              &
1934                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
1935                              * surf_t_surface_p%var_1d(m) )
1936
1937          surf%qsws_soil(m) = - f_qsws_soil * ( surf%qv1(m) - q_s              &
1938                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
1939                              * surf_t_surface_p%var_1d(m) )
1940
1941          surf%qsws_liq(m)  = - f_qsws_liq  * ( surf%qv1(m) - q_s              &
1942                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
1943                              * surf_t_surface_p%var_1d(m) )
1944       ENDIF
1945
1946!
1947!--    Calculate the true surface resistance. ABS is used here to avoid negative
1948!--    values that can occur for very small fluxes due to the artifical addition
1949!--    of 1.0E-20.
1950       IF ( .NOT.  humidity )  THEN
1951          surf%r_s(m) = 1.0E10_wp
1952       ELSE
1953          surf%r_s(m) = ABS(rho_lv / (f_qsws + 1.0E-20_wp) - surf%r_a(m))
1954       ENDIF
1955!
1956!--    Calculate change in liquid water reservoir due to dew fall or
1957!--    evaporation of liquid water
1958       IF ( humidity )  THEN
1959!
1960!--       If precipitation is activated, add rain water to qsws_liq
1961!--       and qsws_soil according the the vegetation coverage.
1962!--       precipitation_rate is given in mm.
1963          IF ( precipitation )  THEN
1964
1965!
1966!--          Add precipitation to liquid water reservoir, if possible.
1967!--          Otherwise, add the water to soil. In case of
1968!--          pavements, the exceeding water amount is explicitly removed
1969!--          (as fictive runoff by drainage systems)
1970             IF ( surf%pavement_surface(m) )  THEN
1971                IF ( surf_m_liq%var_1d(m) < m_liq_max )  THEN
1972                   surf%qsws_liq(m) = surf%qsws_liq(m)                         &
1973                                 + prr(k+k_off,j+j_off,i+i_off)                &
1974                                 * hyrho(k+k_off)                              &
1975                                 * 0.001_wp * rho_l * l_v
1976                ENDIF         
1977             ELSE
1978                IF ( surf_m_liq%var_1d(m) < m_liq_max )  THEN
1979                   surf%qsws_liq(m) = surf%qsws_liq(m)                         &
1980                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
1981                                 * hyrho(k+k_off)                              &
1982                                 * 0.001_wp * rho_l * l_v
1983                   surf%qsws_soil(m) = surf%qsws_soil(m) + ( 1.0_wp -          &
1984                                 surf%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
1985                                 * hyrho(k+k_off)                              &
1986                                 * 0.001_wp * rho_l * l_v                                 
1987                ELSE
1988
1989!--                Add precipitation to bare soil according to the bare soil
1990!--                coverage.
1991                   surf%qsws_soil(m) = surf%qsws_soil(m)                       &
1992                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
1993                                 * hyrho(k+k_off)                              &
1994                                 * 0.001_wp * rho_l * l_v
1995
1996                ENDIF
1997             ENDIF
1998
1999          ENDIF
2000
2001!
2002!--       If the air is saturated, check the reservoir water level
2003          IF ( surf%qsws(m) < 0.0_wp )  THEN
2004!
2005!--          Check if reservoir is full (avoid values > m_liq_max)
2006!--          In that case, qsws_liq goes to qsws_soil for pervious surfaces. In
2007!--          this case qsws_veg is zero anyway (because c_liq = 1),       
2008!--          so that tend is zero and no further check is needed
2009             IF ( surf_m_liq%var_1d(m) == m_liq_max )  THEN
2010                IF ( .NOT. surf%pavement_surface(m))  THEN
2011                   surf%qsws_soil(m) = surf%qsws_soil(m) + surf%qsws_liq(m)
2012                ENDIF
2013                surf%qsws_liq(m)  = 0.0_wp
2014             ENDIF
2015
2016!
2017!--          In case qsws_veg becomes negative (unphysical behavior),
2018!--          let the water enter the liquid water reservoir as dew on the
2019!--          plant
2020             IF ( surf%qsws_veg(m) < 0.0_wp )  THEN
2021                surf%qsws_liq(m) = surf%qsws_liq(m) + surf%qsws_veg(m)
2022                surf%qsws_veg(m) = 0.0_wp
2023             ENDIF
2024          ENDIF                   
2025 
2026          surf%qsws(m) = surf%qsws(m) / l_v
2027 
2028          tend = - surf%qsws_liq(m) * drho_l_lv
2029          surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d *              &
2030                                        ( tsc(2) * tend +                      &
2031                                          tsc(3) * surf_tm_liq_m%var_1d(m) )
2032!
2033!--       Check if reservoir is overfull -> reduce to maximum
2034!--       (conservation of water is violated here)
2035          surf_m_liq_p%var_1d(m) = MIN( surf_m_liq_p%var_1d(m),m_liq_max )
2036
2037!
2038!--       Check if reservoir is empty (avoid values < 0.0)
2039!--       (conservation of water is violated here)
2040          surf_m_liq_p%var_1d(m) = MAX( surf_m_liq_p%var_1d(m), 0.0_wp )
2041!
2042!--       Calculate m_liq tendencies for the next Runge-Kutta step
2043          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2044             IF ( intermediate_timestep_count == 1 )  THEN
2045                surf_tm_liq_m%var_1d(m) = tend
2046             ELSEIF ( intermediate_timestep_count <                            &
2047                      intermediate_timestep_count_max )  THEN
2048                surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend +                  &
2049                                           5.3125_wp * surf_tm_liq_m%var_1d(m)
2050             ENDIF
2051          ENDIF
2052
2053       ENDIF
2054
2055    ENDDO
2056    !$OMP END PARALLEL
2057
2058!
2059!-- Make a logical OR for all processes. Force radiation call if at
2060!-- least one processor reached the threshold change in skin temperature
2061    IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count       &
2062         == intermediate_timestep_count_max-1 )  THEN
2063#if defined( __parallel )
2064       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2065       CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,       &
2066                           1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
2067#else
2068       force_radiation_call = force_radiation_call_l
2069#endif
2070       force_radiation_call_l = .FALSE.
2071    ENDIF
2072
2073!
2074!-- Calculate surface water vapor mixing ratio
2075    IF ( humidity )  THEN
2076       CALL calc_q_surface
2077    ENDIF
2078
2079!
2080!-- Calculate new roughness lengths (for water surfaces only)
2081    IF ( horizontal  .AND.  .NOT. constant_roughness )  CALL calc_z0_water_surface
2082   
2083    IF ( debug_output_timestep )  THEN
2084       WRITE( debug_string, * ) 'lsm_energy_balance', horizontal, l
2085       CALL debug_message( debug_string, 'end' )
2086    ENDIF
2087
2088    CONTAINS
2089!------------------------------------------------------------------------------!
2090! Description:
2091! ------------
2092!> Calculation of mixing ratio of the skin layer (surface). It is assumend
2093!> that the skin is always saturated.
2094!------------------------------------------------------------------------------!
2095    SUBROUTINE calc_q_surface
2096
2097       IMPLICIT NONE
2098
2099       REAL(wp) ::  e_s           !< saturation water vapor pressure
2100       REAL(wp) ::  q_s           !< saturation mixing ratio
2101       REAL(wp) ::  resistance    !< aerodynamic and soil resistance term
2102
2103
2104       !$OMP PARALLEL PRIVATE (m, i, j, k, e_s, q_s, resistance)
2105       !$OMP DO SCHEDULE (STATIC)
2106       DO  m = 1, surf%ns
2107
2108          i   = surf%i(m)           
2109          j   = surf%j(m)
2110          k   = surf%k(m)
2111!
2112!--       Calculate water vapour pressure at saturation and convert to hPa
2113          e_s = 0.01_wp * magnus( MIN(surf_t_surface_p%var_1d(m), 333.15_wp) )                   
2114
2115!
2116!--       Calculate mixing ratio at saturation
2117          q_s = rd_d_rv * e_s / ( surface_pressure - e_s )
2118
2119          resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) + 1E-5_wp )
2120
2121!
2122!--       Calculate mixing ratio at surface
2123          IF ( bulk_cloud_model )  THEN
2124             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
2125                                        ( 1.0_wp - resistance ) *              &
2126                                        ( q(k,j,i) - ql(k,j,i) )
2127          ELSE
2128             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
2129                                        ( 1.0_wp - resistance ) *              &
2130                                          q(k,j,i)
2131          ENDIF
2132         
2133          surf%q_surface(m) = q(k+k_off,j+j_off,i+i_off)
2134!
2135!--       Update virtual potential temperature
2136          surf%vpt_surface(m) = surf%pt_surface(m) *                           &
2137                                  ( 1.0_wp + 0.61_wp * surf%q_surface(m) )
2138
2139       
2140                     
2141       ENDDO
2142       !$OMP END PARALLEL
2143 
2144    END SUBROUTINE calc_q_surface
2145       
2146 END SUBROUTINE lsm_energy_balance
2147   
2148   
2149
2150!------------------------------------------------------------------------------!
2151! Description:
2152! ------------
2153!> Header output for land surface model
2154!------------------------------------------------------------------------------!
2155    SUBROUTINE lsm_header ( io )
2156
2157
2158       IMPLICIT NONE
2159
2160       CHARACTER (LEN=86) ::  t_soil_chr          !< String for soil temperature profile
2161       CHARACTER (LEN=86) ::  roots_chr           !< String for root profile
2162       CHARACTER (LEN=86) ::  vertical_index_chr  !< String for the vertical index
2163       CHARACTER (LEN=86) ::  m_soil_chr          !< String for soil moisture
2164       CHARACTER (LEN=86) ::  soil_depth_chr      !< String for soil depth
2165       CHARACTER (LEN=10) ::  coor_chr            !< Temporary string
2166   
2167       INTEGER(iwp) ::  i                         !< Loop index over soil layers
2168 
2169       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
2170 
2171       t_soil_chr = ''
2172       m_soil_chr    = ''
2173       soil_depth_chr  = '' 
2174       roots_chr        = '' 
2175       vertical_index_chr   = ''
2176
2177       i = 1
2178       DO i = nzb_soil, nzt_soil
2179          WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i)
2180          t_soil_chr = TRIM( t_soil_chr ) // ' ' // TRIM( coor_chr )
2181
2182          WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i)
2183          m_soil_chr = TRIM( m_soil_chr ) // ' ' // TRIM( coor_chr )
2184
2185          WRITE (coor_chr,'(F10.2,7X)')  - zs(i)
2186          soil_depth_chr = TRIM( soil_depth_chr ) // ' '  // TRIM( coor_chr )
2187
2188          WRITE (coor_chr,'(F10.2,7X)')  root_fraction(i)
2189          roots_chr = TRIM( roots_chr ) // ' '  // TRIM( coor_chr )
2190
2191          WRITE (coor_chr,'(I10,7X)')  i
2192          vertical_index_chr = TRIM( vertical_index_chr ) // ' '  //           &
2193                               TRIM( coor_chr )
2194       ENDDO
2195
2196!
2197!--    Write land surface model header
2198       WRITE( io,  1 )
2199       IF ( conserve_water_content )  THEN
2200          WRITE( io, 2 )
2201       ELSE
2202          WRITE( io, 3 )
2203       ENDIF
2204
2205       IF ( vegetation_type_f%from_file )  THEN
2206          WRITE( io, 5 )
2207       ELSE
2208          WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ),        &
2209                         TRIM (soil_type_name(soil_type) )
2210       ENDIF
2211       WRITE( io, 6 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
2212                        TRIM( m_soil_chr ), TRIM( roots_chr ),                 &
2213                        TRIM( vertical_index_chr )
2214
22151   FORMAT (//' Land surface model information:'/                              &
2216              ' ------------------------------'/)
22172   FORMAT ('    --> Soil bottom is closed (water content is conserved',       &
2218            ', default)')
22193   FORMAT ('    --> Soil bottom is open (water content is not conserved)')         
22204   FORMAT ('    --> Land surface type  : ',A,/                                &
2221            '    --> Soil porosity type : ',A)
22225   FORMAT ('    --> Land surface type  : read from file' /                    &
2223            '    --> Soil porosity type : read from file' )
22246   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
2225            '       Height:        ',A,'  m'/                                  &
2226            '       Temperature:   ',A,'  K'/                                  &
2227            '       Moisture:      ',A,'  m**3/m**3'/                          &
2228            '       Root fraction: ',A,'  '/                                   &
2229            '       Grid point:    ',A)
2230
2231
2232    END SUBROUTINE lsm_header
2233
2234
2235!------------------------------------------------------------------------------!
2236! Description:
2237! ------------
2238!> Initialization of the land surface model
2239!------------------------------------------------------------------------------!
2240    SUBROUTINE lsm_init
2241
2242       USE control_parameters,                                                 &
2243           ONLY:  message_string
2244
2245       USE indices,                                                            &
2246           ONLY:  nx, ny, topo_min_level
2247
2248       USE pmc_interface,                                                      &
2249           ONLY:  nested_run
2250   
2251       IMPLICIT NONE
2252
2253       LOGICAL      ::  init_msoil_from_parent   !< flag controlling initialization of soil moisture in nested child domains
2254       LOGICAL      ::  init_tsoil_from_parent   !< flag controlling initialization of soil temperature in nested child domains
2255
2256       INTEGER(iwp) ::  i                       !< running index
2257       INTEGER(iwp) ::  j                       !< running index
2258       INTEGER(iwp) ::  k                       !< running index
2259       INTEGER(iwp) ::  kn                      !< running index
2260       INTEGER(iwp) ::  ko                      !< running index
2261       INTEGER(iwp) ::  kroot                   !< running index
2262       INTEGER(iwp) ::  kzs                     !< running index
2263       INTEGER(iwp) ::  l                       !< running index surface facing
2264       INTEGER(iwp) ::  m                       !< running index
2265       INTEGER(iwp) ::  st                      !< soil-type index
2266       INTEGER(iwp) ::  n_soil_layers_total     !< temperature variable, stores the total number of soil layers + 4
2267
2268       REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
2269       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_soil_init !< temporary array used for averaging soil profiles
2270
2271       IF ( debug_output )  CALL debug_message( 'lsm_init', 'start' )
2272!
2273!--    If no cloud physics is used, rho_surface has not been calculated before
2274       IF (  .NOT.  bulk_cloud_model  .AND.  .NOT.  cloud_droplets )  THEN
2275          CALL calc_mean_profile( pt, 4 )
2276          rho_surface = hyp(nzb) / ( r_d * hom(topo_min_level+1,1,4,0) * exner(nzb) )
2277       ENDIF
2278
2279!
2280!--    Calculate frequently used parameters
2281       rho_cp    = c_p * rho_surface
2282       rho_lv    = rho_surface * l_v
2283       drho_l_lv = 1.0_wp / (rho_l * l_v)
2284
2285!
2286!--    Set initial values for prognostic quantities
2287!--    Horizontal surfaces
2288       tt_surface_h_m%var_1d = 0.0_wp
2289       tt_soil_h_m%var_2d    = 0.0_wp
2290       tm_soil_h_m%var_2d    = 0.0_wp
2291       tm_liq_h_m%var_1d     = 0.0_wp
2292       surf_lsm_h%c_liq      = 0.0_wp
2293
2294       surf_lsm_h%ghf = 0.0_wp
2295
2296       surf_lsm_h%qsws_liq  = 0.0_wp
2297       surf_lsm_h%qsws_soil = 0.0_wp
2298       surf_lsm_h%qsws_veg  = 0.0_wp
2299
2300       surf_lsm_h%r_a        = 50.0_wp
2301       surf_lsm_h%r_s        = 50.0_wp
2302       surf_lsm_h%r_canopy   = 0.0_wp
2303       surf_lsm_h%r_soil     = 0.0_wp
2304!
2305!--    Do the same for vertical surfaces
2306       DO  l = 0, 3
2307          tt_surface_v_m(l)%var_1d = 0.0_wp
2308          tt_soil_v_m(l)%var_2d    = 0.0_wp
2309          tm_soil_v_m(l)%var_2d    = 0.0_wp
2310          tm_liq_v_m(l)%var_1d     = 0.0_wp
2311          surf_lsm_v(l)%c_liq      = 0.0_wp
2312
2313          surf_lsm_v(l)%ghf = 0.0_wp
2314
2315          surf_lsm_v(l)%qsws_liq  = 0.0_wp
2316          surf_lsm_v(l)%qsws_soil = 0.0_wp
2317          surf_lsm_v(l)%qsws_veg  = 0.0_wp
2318
2319          surf_lsm_v(l)%r_a        = 50.0_wp
2320          surf_lsm_v(l)%r_s        = 50.0_wp
2321          surf_lsm_v(l)%r_canopy   = 0.0_wp
2322          surf_lsm_v(l)%r_soil     = 0.0_wp
2323       ENDDO
2324
2325!
2326!--    Set initial values for prognostic soil quantities
2327       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2328          t_soil_h%var_2d = 0.0_wp
2329          m_soil_h%var_2d = 0.0_wp
2330          m_liq_h%var_1d  = 0.0_wp
2331
2332          DO  l = 0, 3
2333             t_soil_v(l)%var_2d = 0.0_wp
2334             m_soil_v(l)%var_2d = 0.0_wp
2335             m_liq_v(l)%var_1d  = 0.0_wp
2336          ENDDO
2337       ENDIF
2338!
2339!--    Allocate 3D soil model arrays
2340!--    First, for horizontal surfaces
2341       ALLOCATE ( surf_lsm_h%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
2342       ALLOCATE ( surf_lsm_h%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2343       ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
2344       ALLOCATE ( surf_lsm_h%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns))
2345       ALLOCATE ( surf_lsm_h%l_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
2346       ALLOCATE ( surf_lsm_h%m_fc(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
2347       ALLOCATE ( surf_lsm_h%m_res(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
2348       ALLOCATE ( surf_lsm_h%m_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
2349       ALLOCATE ( surf_lsm_h%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_h%ns)      )
2350       ALLOCATE ( surf_lsm_h%n_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
2351       ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2352       ALLOCATE ( surf_lsm_h%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2353       ALLOCATE ( surf_lsm_h%root_fr(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
2354   
2355       surf_lsm_h%lambda_h     = 0.0_wp
2356!
2357!--    If required, allocate humidity-related variables for the soil model
2358       IF ( humidity )  THEN
2359          ALLOCATE ( surf_lsm_h%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2360          ALLOCATE ( surf_lsm_h%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  ) 
2361
2362          surf_lsm_h%lambda_w = 0.0_wp 
2363       ENDIF
2364!
2365!--    For vertical surfaces
2366       DO  l = 0, 3
2367          ALLOCATE ( surf_lsm_v(l)%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
2368          ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2369          ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
2370          ALLOCATE ( surf_lsm_v(l)%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns))
2371          ALLOCATE ( surf_lsm_v(l)%l_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
2372          ALLOCATE ( surf_lsm_v(l)%m_fc(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
2373          ALLOCATE ( surf_lsm_v(l)%m_res(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
2374          ALLOCATE ( surf_lsm_v(l)%m_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
2375          ALLOCATE ( surf_lsm_v(l)%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)      )
2376          ALLOCATE ( surf_lsm_v(l)%n_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
2377          ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) ) 
2378          ALLOCATE ( surf_lsm_v(l)%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) ) 
2379          ALLOCATE ( surf_lsm_v(l)%root_fr(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )
2380
2381          surf_lsm_v(l)%lambda_h     = 0.0_wp 
2382         
2383!
2384!--       If required, allocate humidity-related variables for the soil model
2385          IF ( humidity )  THEN
2386             ALLOCATE ( surf_lsm_v(l)%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2387             ALLOCATE ( surf_lsm_v(l)%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  ) 
2388
2389             surf_lsm_v(l)%lambda_w = 0.0_wp 
2390          ENDIF     
2391       ENDDO
2392!
2393!--    Allocate albedo type and emissivity for vegetation, water and pavement
2394!--    fraction.
2395!--    Set default values at each surface element.
2396       ALLOCATE ( surf_lsm_h%albedo_type(0:2,1:surf_lsm_h%ns) )
2397       ALLOCATE ( surf_lsm_h%emissivity(0:2,1:surf_lsm_h%ns) )
2398!
2399!--    Initialize albedo type according to its default type, in order to set values
2400!--    independent on default albedo_type in radiation model.
2401       surf_lsm_h%albedo_type(ind_veg_wall,:)  =                               &
2402                             INT( vegetation_pars(ind_v_at,vegetation_type) )
2403       surf_lsm_h%albedo_type(ind_wat_win,:)   =                               &
2404                             INT( water_pars(ind_w_at,water_type)           )
2405       surf_lsm_h%albedo_type(ind_pav_green,:) =                               &
2406                             INT( pavement_pars(ind_p_at,pavement_type)     )
2407       surf_lsm_h%emissivity  = emissivity
2408       DO  l = 0, 3
2409          ALLOCATE ( surf_lsm_v(l)%albedo_type(0:2,1:surf_lsm_v(l)%ns) )
2410          ALLOCATE ( surf_lsm_v(l)%emissivity(0:2,1:surf_lsm_v(l)%ns)  )
2411!
2412!--       Initialize albedo type according to its default type, in order to
2413!--       set values independent on default albedo_type in radiation model.
2414          surf_lsm_v(l)%albedo_type(ind_veg_wall,:)  =                         &
2415                             INT( vegetation_pars(ind_v_at,vegetation_type) )
2416          surf_lsm_v(l)%albedo_type(ind_wat_win,:)   =                         &
2417                             INT( water_pars(ind_w_at,water_type)           )
2418          surf_lsm_v(l)%albedo_type(ind_pav_green,:) =                         &
2419                             INT( pavement_pars(ind_p_at,pavement_type)     )
2420          surf_lsm_v(l)%emissivity  = emissivity
2421       ENDDO
2422!
2423!--    Allocate arrays for relative surface fraction.
2424!--    0 - vegetation fraction, 2 - water fraction, 1 - pavement fraction
2425       ALLOCATE( surf_lsm_h%frac(0:2,1:surf_lsm_h%ns) )
2426       surf_lsm_h%frac = 0.0_wp
2427       DO  l = 0, 3
2428          ALLOCATE( surf_lsm_v(l)%frac(0:2,1:surf_lsm_v(l)%ns) )
2429          surf_lsm_v(l)%frac = 0.0_wp
2430       ENDDO
2431!
2432!--    For vertical walls only - allocate special flag indicating if any building is on
2433!--    top of any natural surfaces. Used for initialization only.
2434       DO  l = 0, 3
2435          ALLOCATE( surf_lsm_v(l)%building_covered(1:surf_lsm_v(l)%ns) )
2436       ENDDO
2437!
2438!--    Allocate arrays for the respective types and their names on the surface
2439!--    elements. This will be required to treat deposition of chemical species.
2440       ALLOCATE( surf_lsm_h%pavement_type(1:surf_lsm_h%ns)   )
2441       ALLOCATE( surf_lsm_h%vegetation_type(1:surf_lsm_h%ns) )
2442       ALLOCATE( surf_lsm_h%water_type(1:surf_lsm_h%ns)      )
2443       
2444       surf_lsm_h%pavement_type   = 0
2445       surf_lsm_h%vegetation_type = 0
2446       surf_lsm_h%water_type      = 0
2447       
2448       ALLOCATE( surf_lsm_h%pavement_type_name(1:surf_lsm_h%ns)   )
2449       ALLOCATE( surf_lsm_h%vegetation_type_name(1:surf_lsm_h%ns) )
2450       ALLOCATE( surf_lsm_h%water_type_name(1:surf_lsm_h%ns)      )
2451       
2452       surf_lsm_h%pavement_type_name   = 'none'
2453       surf_lsm_h%vegetation_type_name = 'none'
2454       surf_lsm_h%water_type_name      = 'none'
2455       
2456       DO  l = 0, 3
2457          ALLOCATE( surf_lsm_v(l)%pavement_type(1:surf_lsm_v(l)%ns)   )
2458          ALLOCATE( surf_lsm_v(l)%vegetation_type(1:surf_lsm_v(l)%ns) )
2459          ALLOCATE( surf_lsm_v(l)%water_type(1:surf_lsm_v(l)%ns)      )
2460         
2461          surf_lsm_v(l)%pavement_type   = 0
2462          surf_lsm_v(l)%vegetation_type = 0
2463          surf_lsm_v(l)%water_type      = 0
2464       
2465          ALLOCATE( surf_lsm_v(l)%pavement_type_name(1:surf_lsm_v(l)%ns)   )
2466          ALLOCATE( surf_lsm_v(l)%vegetation_type_name(1:surf_lsm_v(l)%ns) )
2467          ALLOCATE( surf_lsm_v(l)%water_type_name(1:surf_lsm_v(l)%ns)      )
2468       
2469          surf_lsm_v(l)%pavement_type_name   = 'none'
2470          surf_lsm_v(l)%vegetation_type_name = 'none'
2471          surf_lsm_v(l)%water_type_name      = 'none'       
2472       ENDDO
2473       
2474!
2475!--    Set flag parameter for the prescribed surface type depending on user
2476!--    input. Set surface fraction to 1 for the respective type.
2477       SELECT CASE ( TRIM( surface_type ) )
2478         
2479          CASE ( 'vegetation' )
2480         
2481             surf_lsm_h%vegetation_surface = .TRUE.
2482             surf_lsm_h%frac(ind_veg_wall,:) = 1.0_wp
2483             DO  l = 0, 3
2484                surf_lsm_v(l)%vegetation_surface = .TRUE.
2485                surf_lsm_v(l)%frac(ind_veg_wall,:) = 1.0_wp
2486             ENDDO
2487   
2488          CASE ( 'water' )
2489             
2490             surf_lsm_h%water_surface = .TRUE.
2491             surf_lsm_h%frac(ind_wat_win,:) = 1.0_wp
2492!
2493!--          Note, vertical water surface does not really make sense.
2494             DO  l = 0, 3 
2495                surf_lsm_v(l)%water_surface   = .TRUE.
2496                surf_lsm_v(l)%frac(ind_wat_win,:) = 1.0_wp
2497             ENDDO
2498
2499          CASE ( 'pavement' )
2500             
2501             surf_lsm_h%pavement_surface = .TRUE.
2502                surf_lsm_h%frac(ind_pav_green,:) = 1.0_wp
2503             DO  l = 0, 3
2504                surf_lsm_v(l)%pavement_surface   = .TRUE.
2505                surf_lsm_v(l)%frac(ind_pav_green,:) = 1.0_wp
2506             ENDDO
2507
2508          CASE ( 'netcdf' )
2509
2510             DO  m = 1, surf_lsm_h%ns
2511                i = surf_lsm_h%i(m)
2512                j = surf_lsm_h%j(m)
2513                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )    &
2514                   surf_lsm_h%vegetation_surface(m) = .TRUE.
2515                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )      &
2516                   surf_lsm_h%pavement_surface(m) = .TRUE.
2517                IF ( water_type_f%var(j,i)      /= water_type_f%fill )         &
2518                   surf_lsm_h%water_surface(m) = .TRUE.
2519!
2520!--             Check if at least one type is set.
2521                IF ( .NOT. surf_lsm_h%vegetation_surface(m)  .AND.             &
2522                     .NOT. surf_lsm_h%pavement_surface(m)    .AND.             &
2523                     .NOT. surf_lsm_h%water_surface(m) )  THEN
2524                   WRITE( message_string, * ) 'Horizontal surface element ' // &
2525                                       ' at i, j = ',  i, j,                   &
2526                                       ' is neither a vegetation, ' //         &
2527                                       'pavement, nor a water surface.'
2528                   CALL message( 'land_surface_model_mod', 'PA0619',          &
2529                                  2, 2, myid, 6, 0 )
2530                ENDIF
2531               
2532             ENDDO
2533!
2534!--          For vertical surfaces some special checks and treatment are
2535!--          required for correct initialization.
2536             DO  l = 0, 3
2537                DO  m = 1, surf_lsm_v(l)%ns
2538!
2539!--                Only for vertical surfaces. Check if at the grid point where
2540!--                the wall is defined (i+ioff, j+joff) is any building.
2541!--                This case, no natural surfaces properties will be defined at
2542!--                at this grid point, leading to problems in the initialization.
2543!--                To overcome this, define a special flag which
2544!--                indicates that a building is defined at the wall grid point 
2545!--                and take the surface properties from the adjoining grid
2546!--                point, i.e. without offset values.
2547!--                Further, there can occur a special case where elevation
2548!--                changes are larger than building heights. This case, (j,i)
2549!--                and (j+joff,i+ioff) grid points may be both covered by
2550!--                buildings, but vertical, but vertically natural walls may
2551!--                be located between the buildings. This case, it is not
2552!--                guaranteed that information about natural surface types is
2553!--                given, neither at (j,i) nor at (j+joff,i+ioff), again leading
2554!--                to non-initialized surface properties.
2555                   surf_lsm_v(l)%building_covered(m) = .FALSE.
2556!
2557!--                Wall grid point is building-covered. This case, set
2558!--                flag indicating that surface properties are initialized
2559!--                from neighboring reference grid point, which is not
2560!--                building_covered.
2561                   IF ( building_type_f%from_file )  THEN
2562                      i = surf_lsm_v(l)%i(m)
2563                      j = surf_lsm_v(l)%j(m)
2564                      IF ( building_type_f%var(j+surf_lsm_v(l)%joff,           &
2565                                               i+surf_lsm_v(l)%ioff) /=        &
2566                           building_type_f%fill )                              &
2567                         surf_lsm_v(l)%building_covered(m) = .TRUE.
2568!
2569!--                   Wall grid point as well as neighboring reference grid
2570!--                   point are both building-covered. This case, surface
2571!--                   properties are not necessarily defined (not covered by
2572!--                   checks for static input file) at this surface. Hence,
2573!--                   initialize surface properties by simply setting
2574!--                   vegetation_type_f to bare-soil bulk parametrization.
2575!--                   soil_type_f as well as surface_fractions_f will be set
2576!--                   also.                     
2577                      IF ( building_type_f%var(j+surf_lsm_v(l)%joff,           &
2578                                               i+surf_lsm_v(l)%ioff) /=        &
2579                           building_type_f%fill  .AND.                         &
2580                           building_type_f%var(j,i) /= building_type_f%fill )  &
2581                      THEN
2582                         vegetation_type_f%var(j,i)                 = 1 ! bare soil
2583                         soil_type_f%var_2d(j,i)                    = 1 
2584!                         
2585!--                      Set surface_fraction if provided in static input,
2586!--                      else, in case no tiles are used, this will be done
2587!--                      on basis of the prescribed types (vegetation/pavement/
2588!--                      water_type).
2589                         IF ( surface_fraction_f%from_file )  THEN
2590                            surface_fraction_f%frac(ind_veg_wall,j,i)  = 1.0_wp
2591                            surface_fraction_f%frac(ind_pav_green,j,i) = 0.0_wp
2592                            surface_fraction_f%frac(ind_wat_win,j,i)   = 0.0_wp
2593                         ENDIF
2594                      ENDIF
2595                     
2596                   ENDIF
2597!
2598!--                Normally proceed with setting surface types.
2599                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2600                                            surf_lsm_v(l)%building_covered(m) )
2601                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2602                                            surf_lsm_v(l)%building_covered(m) )
2603                   IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
2604                      surf_lsm_v(l)%vegetation_surface(m) = .TRUE.
2605                   IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )   &
2606                      surf_lsm_v(l)%pavement_surface(m) = .TRUE.
2607                   IF ( water_type_f%var(j,i)      /= water_type_f%fill )      &
2608                      surf_lsm_v(l)%water_surface(m) = .TRUE.
2609!
2610!--                Check if at least one type is set.
2611                   IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m)  .AND.       &
2612                        .NOT. surf_lsm_v(l)%pavement_surface(m)    .AND.       &
2613                        .NOT. surf_lsm_v(l)%water_surface(m) )  THEN
2614                      WRITE( message_string, * ) 'Vertical surface element ' //&
2615                                       ' at i, j = ',  i, j,                   &
2616                                       ' is neither a vegetation, ' //         &
2617                                       'pavement, nor a water surface.'
2618                      CALL message( 'land_surface_model_mod', 'PA0619',        &
2619                                     2, 2, myid, 6, 0 )
2620                   ENDIF
2621                ENDDO
2622             ENDDO
2623
2624       END SELECT
2625!
2626!--    In case of netcdf input file, further initialize surface fractions.
2627!--    At the moment only 1 surface is given at a location, so that the fraction
2628!--    is either 0 or 1. This will be revised later. If surface fraction
2629!--    is not given in static input file, relative fractions will be derived
2630!--    from given surface type. In this case, only 1 type is given at a certain
2631!--    location (already checked). 
2632       IF ( input_pids_static  .AND.  surface_fraction_f%from_file )  THEN
2633          DO  m = 1, surf_lsm_h%ns
2634             i = surf_lsm_h%i(m)
2635             j = surf_lsm_h%j(m)
2636!
2637!--          0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction             
2638             surf_lsm_h%frac(ind_veg_wall,m)  =                                &
2639                                    surface_fraction_f%frac(ind_veg_wall,j,i)         
2640             surf_lsm_h%frac(ind_pav_green,m) =                                &
2641                                    surface_fraction_f%frac(ind_pav_green,j,i)       
2642             surf_lsm_h%frac(ind_wat_win,m)   =                                &
2643                                    surface_fraction_f%frac(ind_wat_win,j,i)
2644
2645          ENDDO
2646          DO  l = 0, 3
2647             DO  m = 1, surf_lsm_v(l)%ns
2648                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
2649                                                surf_lsm_v(l)%building_covered(m) ) 
2650                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
2651                                                surf_lsm_v(l)%building_covered(m) ) 
2652!
2653!--             0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction       
2654                surf_lsm_v(l)%frac(ind_veg_wall,m)  =                          &
2655                                    surface_fraction_f%frac(ind_veg_wall,j,i)         
2656                surf_lsm_v(l)%frac(ind_pav_green,m) =                          &
2657                                    surface_fraction_f%frac(ind_pav_green,j,i)       
2658                surf_lsm_v(l)%frac(ind_wat_win,m)   =                          &
2659                                    surface_fraction_f%frac(ind_wat_win,j,i)
2660
2661             ENDDO
2662          ENDDO
2663       ELSEIF ( input_pids_static  .AND.  .NOT. surface_fraction_f%from_file ) &
2664       THEN
2665          DO  m = 1, surf_lsm_h%ns
2666             i = surf_lsm_h%i(m)
2667             j = surf_lsm_h%j(m)
2668
2669             IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )       &       
2670                surf_lsm_h%frac(ind_veg_wall,m)  = 1.0_wp
2671             IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill   )       &       
2672                surf_lsm_h%frac(ind_pav_green,m) = 1.0_wp 
2673             IF ( water_type_f%var(j,i)      /= water_type_f%fill      )       &       
2674                surf_lsm_h%frac(ind_wat_win,m)   = 1.0_wp       
2675          ENDDO
2676          DO  l = 0, 3
2677             DO  m = 1, surf_lsm_v(l)%ns
2678                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
2679                                                surf_lsm_v(l)%building_covered(m) ) 
2680                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
2681                                                surf_lsm_v(l)%building_covered(m) ) 
2682     
2683                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )    &       
2684                   surf_lsm_v(l)%frac(ind_veg_wall,m)  = 1.0_wp
2685                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill   )    &       
2686                   surf_lsm_v(l)%frac(ind_pav_green,m) = 1.0_wp 
2687                IF ( water_type_f%var(j,i)      /= water_type_f%fill      )    &       
2688                   surf_lsm_v(l)%frac(ind_wat_win,m)   = 1.0_wp     
2689             ENDDO
2690          ENDDO
2691       ENDIF
2692!
2693!--    Level 1, initialization of soil parameters.
2694!--    It is possible to overwrite each parameter by setting the respecticy
2695!--    NAMELIST variable to a value /= 9999999.9.
2696       IF ( soil_type /= 0 )  THEN 
2697 
2698          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
2699             alpha_vangenuchten = soil_pars(0,soil_type)
2700          ENDIF
2701
2702          IF ( l_vangenuchten == 9999999.9_wp )  THEN
2703             l_vangenuchten = soil_pars(1,soil_type)
2704          ENDIF
2705
2706          IF ( n_vangenuchten == 9999999.9_wp )  THEN
2707             n_vangenuchten = soil_pars(2,soil_type)           
2708          ENDIF
2709
2710          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
2711             hydraulic_conductivity = soil_pars(3,soil_type)           
2712          ENDIF
2713
2714          IF ( saturation_moisture == 9999999.9_wp )  THEN
2715             saturation_moisture = soil_pars(4,soil_type)           
2716          ENDIF
2717
2718          IF ( field_capacity == 9999999.9_wp )  THEN
2719             field_capacity = soil_pars(5,soil_type)           
2720          ENDIF
2721
2722          IF ( wilting_point == 9999999.9_wp )  THEN
2723             wilting_point = soil_pars(6,soil_type)           
2724          ENDIF
2725
2726          IF ( residual_moisture == 9999999.9_wp )  THEN
2727             residual_moisture = soil_pars(7,soil_type)       
2728          ENDIF
2729
2730       ENDIF
2731!
2732!--    Map values to the respective 2D/3D arrays
2733!--    Horizontal surfaces
2734       surf_lsm_h%alpha_vg      = alpha_vangenuchten
2735       surf_lsm_h%l_vg          = l_vangenuchten
2736       surf_lsm_h%n_vg          = n_vangenuchten 
2737       surf_lsm_h%gamma_w_sat   = hydraulic_conductivity
2738       surf_lsm_h%m_sat         = saturation_moisture
2739       surf_lsm_h%m_fc          = field_capacity
2740       surf_lsm_h%m_wilt        = wilting_point
2741       surf_lsm_h%m_res         = residual_moisture
2742       surf_lsm_h%r_soil_min    = min_soil_resistance
2743!
2744!--    Vertical surfaces
2745       DO  l = 0, 3
2746          surf_lsm_v(l)%alpha_vg      = alpha_vangenuchten
2747          surf_lsm_v(l)%l_vg          = l_vangenuchten
2748          surf_lsm_v(l)%n_vg          = n_vangenuchten 
2749          surf_lsm_v(l)%gamma_w_sat   = hydraulic_conductivity
2750          surf_lsm_v(l)%m_sat         = saturation_moisture
2751          surf_lsm_v(l)%m_fc          = field_capacity
2752          surf_lsm_v(l)%m_wilt        = wilting_point
2753          surf_lsm_v(l)%m_res         = residual_moisture
2754          surf_lsm_v(l)%r_soil_min    = min_soil_resistance
2755       ENDDO
2756!
2757!--    Level 2, initialization of soil parameters via soil_type read from file.
2758!--    Soil parameters are initialized for each (y,x)-grid point
2759!--    individually using default paramter settings according to the given
2760!--    soil type.
2761       IF ( soil_type_f%from_file )  THEN
2762!
2763!--       Level of detail = 1, i.e. a homogeneous soil distribution along the
2764!--       vertical dimension is assumed.
2765          IF ( soil_type_f%lod == 1 )  THEN
2766!
2767!--          Horizontal surfaces
2768             DO  m = 1, surf_lsm_h%ns
2769                i = surf_lsm_h%i(m)
2770                j = surf_lsm_h%j(m)
2771             
2772                st = soil_type_f%var_2d(j,i)
2773                IF ( st /= soil_type_f%fill )  THEN
2774                   surf_lsm_h%alpha_vg(:,m)    = soil_pars(0,st)
2775                   surf_lsm_h%l_vg(:,m)        = soil_pars(1,st)
2776                   surf_lsm_h%n_vg(:,m)        = soil_pars(2,st)
2777                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars(3,st)
2778                   surf_lsm_h%m_sat(:,m)       = soil_pars(4,st)
2779                   surf_lsm_h%m_fc(:,m)        = soil_pars(5,st)
2780                   surf_lsm_h%m_wilt(:,m)      = soil_pars(6,st)
2781                   surf_lsm_h%m_res(:,m)       = soil_pars(7,st)
2782                ENDIF
2783             ENDDO
2784!
2785!--          Vertical surfaces ( assumes the soil type given at respective (x,y)
2786             DO  l = 0, 3
2787                DO  m = 1, surf_lsm_v(l)%ns
2788                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2789                                                   surf_lsm_v(l)%building_covered(m) ) 
2790                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2791                                                   surf_lsm_v(l)%building_covered(m) ) 
2792
2793                   st = soil_type_f%var_2d(j,i)
2794                   IF ( st /= soil_type_f%fill )  THEN
2795                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars(0,st)
2796                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars(1,st)
2797                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars(2,st)
2798                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars(3,st)
2799                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars(4,st)
2800                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars(5,st)
2801                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars(6,st)
2802                      surf_lsm_v(l)%m_res(:,m)       = soil_pars(7,st)
2803                   ENDIF
2804                ENDDO
2805             ENDDO
2806!
2807!--       Level of detail = 2, i.e. soil type and thus the soil parameters
2808!--       can be heterogeneous along the vertical dimension.
2809          ELSE
2810!
2811!--          Horizontal surfaces
2812             DO  m = 1, surf_lsm_h%ns
2813                i = surf_lsm_h%i(m)
2814                j = surf_lsm_h%j(m)
2815             
2816                DO  k = nzb_soil, nzt_soil
2817                   st = soil_type_f%var_3d(k,j,i)
2818                   IF ( st /= soil_type_f%fill )  THEN
2819                      surf_lsm_h%alpha_vg(k,m)    = soil_pars(0,st)
2820                      surf_lsm_h%l_vg(k,m)        = soil_pars(1,st)
2821                      surf_lsm_h%n_vg(k,m)        = soil_pars(2,st)
2822                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars(3,st)
2823                      surf_lsm_h%m_sat(k,m)       = soil_pars(4,st)
2824                      surf_lsm_h%m_fc(k,m)        = soil_pars(5,st)
2825                      surf_lsm_h%m_wilt(k,m)      = soil_pars(6,st)
2826                      surf_lsm_h%m_res(k,m)       = soil_pars(7,st)
2827                   ENDIF
2828                ENDDO
2829             ENDDO
2830!
2831!--          Vertical surfaces ( assumes the soil type given at respective (x,y)
2832             DO  l = 0, 3
2833                DO  m = 1, surf_lsm_v(l)%ns
2834                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2835                                                   surf_lsm_v(l)%building_covered(m) ) 
2836                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2837                                                   surf_lsm_v(l)%building_covered(m) ) 
2838
2839                   DO  k = nzb_soil, nzt_soil
2840                      st = soil_type_f%var_3d(k,j,i)
2841                      IF ( st /= soil_type_f%fill )  THEN
2842                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars(0,st)
2843                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars(1,st)
2844                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars(2,st)
2845                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars(3,st)
2846                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars(4,st)
2847                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars(5,st)
2848                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars(6,st)
2849                         surf_lsm_v(l)%m_res(k,m)       = soil_pars(7,st)
2850                      ENDIF
2851                   ENDDO
2852                ENDDO
2853             ENDDO
2854          ENDIF
2855       ENDIF
2856!
2857!--    Level 3, initialization of single soil parameters at single z,x,y
2858!--    position via soil_pars read from file.
2859       IF ( soil_pars_f%from_file )  THEN
2860!
2861!--       Level of detail = 1, i.e. a homogeneous vertical distribution of soil
2862!--       parameters is assumed.
2863!--       Horizontal surfaces
2864          IF ( soil_pars_f%lod == 1 )  THEN
2865!
2866!--          Horizontal surfaces
2867             DO  m = 1, surf_lsm_h%ns
2868                i = surf_lsm_h%i(m)
2869                j = surf_lsm_h%j(m)
2870
2871                IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )              &
2872                   surf_lsm_h%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
2873                IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )              &
2874                   surf_lsm_h%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
2875                IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )              &
2876                   surf_lsm_h%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
2877                IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )              &
2878                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
2879                IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )              &
2880                   surf_lsm_h%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
2881                IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )              &
2882                   surf_lsm_h%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
2883                IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )              &
2884                   surf_lsm_h%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
2885                IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )              &
2886                   surf_lsm_h%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)
2887
2888             ENDDO
2889!
2890!--          Vertical surfaces
2891             DO  l = 0, 3
2892                DO  m = 1, surf_lsm_v(l)%ns
2893                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2894                                                   surf_lsm_v(l)%building_covered(m) ) 
2895                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2896                                                   surf_lsm_v(l)%building_covered(m) ) 
2897
2898                   IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )           &
2899                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
2900                   IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )           &
2901                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
2902                   IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )           &
2903                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
2904                   IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )           &
2905                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
2906                   IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )           &
2907                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
2908                   IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )           &
2909                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
2910                   IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )           &
2911                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
2912                   IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )           &
2913                      surf_lsm_v(l)%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)
2914
2915                ENDDO
2916             ENDDO
2917!
2918!--       Level of detail = 2, i.e. soil parameters can be set at each soil
2919!--       layer individually.
2920          ELSE
2921!
2922!--          Horizontal surfaces
2923             DO  m = 1, surf_lsm_h%ns
2924                i = surf_lsm_h%i(m)
2925                j = surf_lsm_h%j(m)
2926
2927                DO  k = nzb_soil, nzt_soil
2928                   IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
2929                      surf_lsm_h%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
2930                   IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
2931                      surf_lsm_h%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
2932                   IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
2933                      surf_lsm_h%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
2934                   IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
2935                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
2936                   IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
2937                      surf_lsm_h%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
2938                   IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
2939                      surf_lsm_h%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
2940                   IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
2941                      surf_lsm_h%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
2942                   IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
2943                      surf_lsm_h%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
2944                ENDDO
2945
2946             ENDDO
2947!
2948!--          Vertical surfaces
2949             DO  l = 0, 3
2950                DO  m = 1, surf_lsm_v(l)%ns
2951                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2952                                                   surf_lsm_v(l)%building_covered(m) ) 
2953                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2954                                                   surf_lsm_v(l)%building_covered(m) ) 
2955
2956                   DO  k = nzb_soil, nzt_soil
2957                      IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
2958                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
2959                      IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
2960                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
2961                      IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
2962                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
2963                      IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
2964                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
2965                      IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
2966                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
2967                      IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
2968                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
2969                      IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
2970                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
2971                      IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
2972                         surf_lsm_v(l)%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
2973                   ENDDO
2974
2975                ENDDO
2976             ENDDO
2977
2978          ENDIF
2979       ENDIF
2980
2981!
2982!--    Level 1, initialization of vegetation parameters. A horizontally
2983!--    homogeneous distribution is assumed here.
2984       IF ( vegetation_type /= 0 )  THEN
2985
2986          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
2987             min_canopy_resistance = vegetation_pars(ind_v_rc_min,vegetation_type)
2988          ENDIF
2989
2990          IF ( leaf_area_index == 9999999.9_wp )  THEN
2991             leaf_area_index = vegetation_pars(ind_v_rc_lai,vegetation_type)         
2992          ENDIF
2993
2994          IF ( vegetation_coverage == 9999999.9_wp )  THEN
2995             vegetation_coverage = vegetation_pars(ind_v_c_veg,vegetation_type)     
2996          ENDIF
2997
2998          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
2999              canopy_resistance_coefficient= vegetation_pars(ind_v_gd,vegetation_type)     
3000          ENDIF
3001
3002          IF ( z0_vegetation == 9999999.9_wp )  THEN
3003             z0_vegetation  = vegetation_pars(ind_v_z0,vegetation_type) 
3004          ENDIF
3005
3006          IF ( z0h_vegetation == 9999999.9_wp )  THEN
3007             z0h_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
3008          ENDIF
3009         
3010          IF ( z0q_vegetation == 9999999.9_wp )  THEN
3011             z0q_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
3012          ENDIF
3013         
3014          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
3015             lambda_surface_stable = vegetation_pars(ind_v_lambda_s,vegetation_type) 
3016          ENDIF
3017
3018          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
3019             lambda_surface_unstable = vegetation_pars(ind_v_lambda_u,vegetation_type)           
3020          ENDIF
3021
3022          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
3023             f_shortwave_incoming = vegetation_pars(ind_v_f_sw_in,vegetation_type)       
3024          ENDIF
3025
3026          IF ( c_surface == 9999999.9_wp )  THEN
3027             c_surface = vegetation_pars(ind_v_c_surf,vegetation_type)       
3028          ENDIF
3029
3030          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
3031             albedo_type = INT(vegetation_pars(ind_v_at,vegetation_type))       
3032          ENDIF
3033   
3034          IF ( emissivity == 9999999.9_wp )  THEN
3035             emissivity = vegetation_pars(ind_v_emis,vegetation_type)     
3036          ENDIF
3037
3038       ENDIF
3039!
3040!--    Map values onto horizontal elemements
3041       DO  m = 1, surf_lsm_h%ns
3042          IF ( surf_lsm_h%vegetation_surface(m) )  THEN
3043             surf_lsm_h%r_canopy_min(m)     = min_canopy_resistance
3044             surf_lsm_h%lai(m)              = leaf_area_index
3045             surf_lsm_h%c_veg(m)            = vegetation_coverage
3046             surf_lsm_h%g_d(m)              = canopy_resistance_coefficient
3047             surf_lsm_h%z0(m)               = z0_vegetation
3048             surf_lsm_h%z0h(m)              = z0h_vegetation
3049             surf_lsm_h%z0q(m)              = z0q_vegetation
3050             surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
3051             surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
3052             surf_lsm_h%f_sw_in(m)          = f_shortwave_incoming
3053             surf_lsm_h%c_surface(m)        = c_surface
3054             surf_lsm_h%albedo_type(ind_veg_wall,m) = albedo_type
3055             surf_lsm_h%emissivity(ind_veg_wall,m)  = emissivity
3056             
3057             surf_lsm_h%vegetation_type(m)      = vegetation_type
3058             surf_lsm_h%vegetation_type_name(m) = vegetation_type_name(vegetation_type)
3059          ELSE
3060             surf_lsm_h%lai(m)   = 0.0_wp
3061             surf_lsm_h%c_veg(m) = 0.0_wp
3062             surf_lsm_h%g_d(m)   = 0.0_wp
3063          ENDIF
3064 
3065       ENDDO
3066!
3067!--    Map values onto vertical elements, even though this does not make
3068!--    much sense.
3069       DO  l = 0, 3
3070          DO  m = 1, surf_lsm_v(l)%ns
3071             IF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
3072                surf_lsm_v(l)%r_canopy_min(m)     = min_canopy_resistance
3073                surf_lsm_v(l)%lai(m)              = leaf_area_index
3074                surf_lsm_v(l)%c_veg(m)            = vegetation_coverage
3075                surf_lsm_v(l)%g_d(m)              = canopy_resistance_coefficient
3076                surf_lsm_v(l)%z0(m)               = z0_vegetation
3077                surf_lsm_v(l)%z0h(m)              = z0h_vegetation
3078                surf_lsm_v(l)%z0q(m)              = z0q_vegetation
3079                surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
3080                surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
3081                surf_lsm_v(l)%f_sw_in(m)          = f_shortwave_incoming
3082                surf_lsm_v(l)%c_surface(m)        = c_surface
3083                surf_lsm_v(l)%albedo_type(ind_veg_wall,m) = albedo_type
3084                surf_lsm_v(l)%emissivity(ind_veg_wall,m)  = emissivity
3085               
3086                surf_lsm_v(l)%vegetation_type(m)      = vegetation_type
3087                surf_lsm_v(l)%vegetation_type_name(m) = vegetation_type_name(vegetation_type)
3088             ELSE
3089                surf_lsm_v(l)%lai(m)   = 0.0_wp
3090                surf_lsm_v(l)%c_veg(m) = 0.0_wp
3091                surf_lsm_v(l)%g_d(m)   = 0.0_wp
3092             ENDIF
3093          ENDDO
3094       ENDDO
3095
3096!
3097!--    Level 2, initialization of vegation parameters via vegetation_type read
3098!--    from file. Vegetation parameters are initialized for each (y,x)-grid point
3099!--    individually using default paramter settings according to the given
3100!--    vegetation type.
3101       IF ( vegetation_type_f%from_file )  THEN
3102!
3103!--       Horizontal surfaces
3104          DO  m = 1, surf_lsm_h%ns
3105             i = surf_lsm_h%i(m)
3106             j = surf_lsm_h%j(m)
3107             
3108             st = vegetation_type_f%var(j,i)
3109             IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
3110                surf_lsm_h%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
3111                surf_lsm_h%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
3112                surf_lsm_h%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
3113                surf_lsm_h%g_d(m)              = vegetation_pars(ind_v_gd,st)
3114                surf_lsm_h%z0(m)               = vegetation_pars(ind_v_z0,st)
3115                surf_lsm_h%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
3116                surf_lsm_h%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
3117                surf_lsm_h%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
3118                surf_lsm_h%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
3119                surf_lsm_h%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
3120                surf_lsm_h%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
3121                surf_lsm_h%albedo_type(ind_veg_wall,m) = INT( vegetation_pars(ind_v_at,st) )
3122                surf_lsm_h%emissivity(ind_veg_wall,m)  = vegetation_pars(ind_v_emis,st)
3123               
3124                surf_lsm_h%vegetation_type(m)      = st
3125                surf_lsm_h%vegetation_type_name(m) = vegetation_type_name(st)
3126             ENDIF
3127          ENDDO
3128!
3129!--       Vertical surfaces
3130          DO  l = 0, 3
3131             DO  m = 1, surf_lsm_v(l)%ns
3132                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3133                                                surf_lsm_v(l)%building_covered(m) ) 
3134                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
3135                                                surf_lsm_v(l)%building_covered(m) ) 
3136             
3137                st = vegetation_type_f%var(j,i)
3138                IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
3139                   surf_lsm_v(l)%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
3140                   surf_lsm_v(l)%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
3141                   surf_lsm_v(l)%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
3142                   surf_lsm_v(l)%g_d(m)              = vegetation_pars(ind_v_gd,st)
3143                   surf_lsm_v(l)%z0(m)               = vegetation_pars(ind_v_z0,st)
3144                   surf_lsm_v(l)%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
3145                   surf_lsm_v(l)%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
3146                   surf_lsm_v(l)%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
3147                   surf_lsm_v(l)%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
3148                   surf_lsm_v(l)%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
3149                   surf_lsm_v(l)%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
3150                   surf_lsm_v(l)%albedo_type(ind_veg_wall,m) = INT( vegetation_pars(ind_v_at,st) )
3151                   surf_lsm_v(l)%emissivity(ind_veg_wall,m)  = vegetation_pars(ind_v_emis,st)
3152                   
3153                   surf_lsm_v(l)%vegetation_type(m)      = st
3154                   surf_lsm_v(l)%vegetation_type_name(m) = vegetation_type_name(st)
3155                ENDIF
3156             ENDDO
3157          ENDDO
3158       ENDIF
3159!
3160!--    Level 3, initialization of vegation parameters at single (x,y)
3161!--    position via vegetation_pars read from file.
3162       IF ( vegetation_pars_f%from_file )  THEN
3163!
3164!--       Horizontal surfaces
3165          DO  m = 1, surf_lsm_h%ns
3166
3167             i = surf_lsm_h%i(m)
3168             j = surf_lsm_h%j(m)
3169!
3170!--          If surface element is not a vegetation surface and any value in
3171!--          vegetation_pars is given, neglect this information and give an
3172!--          informative message that this value will not be used.   
3173             IF ( .NOT. surf_lsm_h%vegetation_surface(m)  .AND.                &
3174                   ANY( vegetation_pars_f%pars_xy(:,j,i) /=                    &
3175                   vegetation_pars_f%fill ) )  THEN
3176                WRITE( message_string, * )                                     &
3177                                 'surface element at grid point (j,i) = (',    &
3178                                 j, i, ') is not a vegetation surface, ',      &
3179                                 'so that information given in ',              &
3180                                 'vegetation_pars at this point is neglected.' 
3181                CALL message( 'land_surface_model_mod', 'PA0436', 0, 0, myid, 6, 0 )
3182             ELSE
3183
3184                IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=            &
3185                     vegetation_pars_f%fill )                                  &
3186                   surf_lsm_h%r_canopy_min(m)  =                               &
3187                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
3188                IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=            &
3189                     vegetation_pars_f%fill )                                  &
3190                   surf_lsm_h%lai(m)           =                               &
3191                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
3192                IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=             &
3193                     vegetation_pars_f%fill )                                  &
3194                   surf_lsm_h%c_veg(m)         =                               &
3195                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
3196                IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=                &
3197                     vegetation_pars_f%fill )                                  &
3198                   surf_lsm_h%g_d(m)           =                               &
3199                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
3200                IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=                &
3201                     vegetation_pars_f%fill )                                  &
3202                   surf_lsm_h%z0(m)            =                               &
3203                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
3204                IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=              &
3205                     vegetation_pars_f%fill )  THEN
3206                   surf_lsm_h%z0h(m)           =                               &
3207                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3208                   surf_lsm_h%z0q(m)           =                               &
3209                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3210                ENDIF
3211                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=          &
3212                     vegetation_pars_f%fill )                                  &
3213                   surf_lsm_h%lambda_surface_s(m) =                            &
3214                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
3215                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=          &
3216                     vegetation_pars_f%fill )                                  &
3217                   surf_lsm_h%lambda_surface_u(m) =                            &
3218                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
3219                IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=           &
3220                     vegetation_pars_f%fill )                                  &
3221                   surf_lsm_h%f_sw_in(m)          =                            &
3222                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
3223                IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=            &
3224                     vegetation_pars_f%fill )                                  &
3225                   surf_lsm_h%c_surface(m)        =                            &
3226                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
3227                IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=                &
3228                     vegetation_pars_f%fill )                                  &
3229                   surf_lsm_h%albedo_type(ind_veg_wall,m) =                    &
3230                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
3231                IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=              &
3232                     vegetation_pars_f%fill )                                  &
3233                   surf_lsm_h%emissivity(ind_veg_wall,m)  =                    &
3234                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
3235             ENDIF
3236          ENDDO
3237!
3238!--       Vertical surfaces
3239          DO  l = 0, 3
3240             DO  m = 1, surf_lsm_v(l)%ns
3241                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3242                                                surf_lsm_v(l)%building_covered(m) ) 
3243                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3244                                                surf_lsm_v(l)%building_covered(m) ) 
3245!
3246!--             If surface element is not a vegetation surface and any value in
3247!--             vegetation_pars is given, neglect this information and give an
3248!--             informative message that this value will not be used.   
3249                IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m)  .AND.          &
3250                      ANY( vegetation_pars_f%pars_xy(:,j,i) /=                 &
3251                      vegetation_pars_f%fill ) )  THEN
3252                   WRITE( message_string, * )                                  &
3253                                 'surface element at grid point (j,i) = (',    &
3254                                 j, i, ') is not a vegetation surface, ',      &
3255                                 'so that information given in ',              &
3256                                 'vegetation_pars at this point is neglected.' 
3257                   CALL message( 'land_surface_model_mod', 'PA0436', 0, 0, myid, 6, 0 )
3258                ELSE
3259
3260                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=         &
3261                        vegetation_pars_f%fill )                               &
3262                      surf_lsm_v(l)%r_canopy_min(m)  =                         &
3263                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
3264                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=         &
3265                        vegetation_pars_f%fill )                               &
3266                      surf_lsm_v(l)%lai(m)           =                         &
3267                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
3268                   IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=          &
3269                        vegetation_pars_f%fill )                               &
3270                      surf_lsm_v(l)%c_veg(m)         =                         &
3271                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
3272                   IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=             &
3273                        vegetation_pars_f%fill )                               &
3274                     surf_lsm_v(l)%g_d(m)            =                         &
3275                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
3276                   IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=             &
3277                        vegetation_pars_f%fill )                               &
3278                      surf_lsm_v(l)%z0(m)            =                         &
3279                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
3280                   IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=           &
3281                        vegetation_pars_f%fill )  THEN
3282                      surf_lsm_v(l)%z0h(m)           =                         &
3283                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3284                      surf_lsm_v(l)%z0q(m)           =                         &
3285                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3286                   ENDIF
3287                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=       &
3288                        vegetation_pars_f%fill )                               &
3289                      surf_lsm_v(l)%lambda_surface_s(m)  =                     &
3290                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
3291                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=       &
3292                        vegetation_pars_f%fill )                               &
3293                      surf_lsm_v(l)%lambda_surface_u(m)  =                     &
3294                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
3295                   IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=        &
3296                        vegetation_pars_f%fill )                               &
3297                      surf_lsm_v(l)%f_sw_in(m)           =                     &
3298                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
3299                   IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=         &
3300                        vegetation_pars_f%fill )                               &
3301                      surf_lsm_v(l)%c_surface(m)         =                     &
3302                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
3303                   IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=             &
3304                        vegetation_pars_f%fill )                               &
3305                      surf_lsm_v(l)%albedo_type(ind_veg_wall,m) =              &
3306                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
3307                   IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=           &
3308                        vegetation_pars_f%fill )                               &
3309                      surf_lsm_v(l)%emissivity(ind_veg_wall,m)  =              &
3310                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
3311                ENDIF
3312
3313             ENDDO
3314          ENDDO
3315       ENDIF 
3316
3317!
3318!--    Level 1, initialization of water parameters. A horizontally
3319!--    homogeneous distribution is assumed here.
3320       IF ( water_type /= 0 )  THEN
3321
3322          IF ( water_temperature == 9999999.9_wp )  THEN
3323             water_temperature = water_pars(ind_w_temp,water_type)       
3324          ENDIF
3325
3326          IF ( z0_water == 9999999.9_wp )  THEN
3327             z0_water = water_pars(ind_w_z0,water_type)       
3328          ENDIF       
3329
3330          IF ( z0h_water == 9999999.9_wp )  THEN
3331             z0h_water = water_pars(ind_w_z0h,water_type)       
3332          ENDIF 
3333         
3334          IF ( z0q_water == 9999999.9_wp )  THEN
3335             z0q_water = water_pars(ind_w_z0h,water_type)       
3336          ENDIF
3337
3338          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
3339             albedo_type = INT(water_pars(ind_w_at,water_type))       
3340          ENDIF
3341   
3342          IF ( emissivity == 9999999.9_wp )  THEN
3343             emissivity = water_pars(ind_w_emis,water_type)       
3344          ENDIF
3345
3346       ENDIF 
3347!
3348!--    Map values onto horizontal elemements
3349       DO  m = 1, surf_lsm_h%ns
3350          IF ( surf_lsm_h%water_surface(m) )  THEN
3351             IF ( TRIM( initializing_actions ) /= 'read_restart_data' )        &
3352                t_soil_h%var_2d(:,m)        = water_temperature
3353             surf_lsm_h%z0(m)               = z0_water
3354             surf_lsm_h%z0h(m)              = z0h_water
3355             surf_lsm_h%z0q(m)              = z0q_water
3356             surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
3357             surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp               
3358             surf_lsm_h%c_surface(m)        = 0.0_wp
3359             surf_lsm_h%albedo_type(ind_wat_win,m) = albedo_type
3360             surf_lsm_h%emissivity(ind_wat_win,m)  = emissivity
3361             
3362             surf_lsm_h%water_type(m)      = water_type
3363             surf_lsm_h%water_type_name(m) = water_type_name(water_type)
3364          ENDIF
3365       ENDDO
3366!
3367!--    Map values onto vertical elements, even though this does not make
3368!--    much sense.
3369       DO  l = 0, 3
3370          DO  m = 1, surf_lsm_v(l)%ns
3371             IF ( surf_lsm_v(l)%water_surface(m) )  THEN
3372                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
3373                   t_soil_v(l)%var_2d(:,m)           = water_temperature
3374                surf_lsm_v(l)%z0(m)               = z0_water
3375                surf_lsm_v(l)%z0h(m)              = z0h_water
3376                surf_lsm_v(l)%z0q(m)              = z0q_water
3377                surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp
3378                surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp               
3379                surf_lsm_v(l)%c_surface(m)        = 0.0_wp
3380                surf_lsm_v(l)%albedo_type(ind_wat_win,m) = albedo_type
3381                surf_lsm_v(l)%emissivity(ind_wat_win,m)  = emissivity
3382               
3383                surf_lsm_v(l)%water_type(m)      = water_type
3384                surf_lsm_v(l)%water_type_name(m) = water_type_name(water_type)
3385             ENDIF
3386          ENDDO
3387       ENDDO
3388!
3389!
3390!--    Level 2, initialization of water parameters via water_type read
3391!--    from file. Water surfaces are initialized for each (y,x)-grid point
3392!--    individually using default paramter settings according to the given
3393!--    water type.
3394!--    Note, parameter 3/4 of water_pars are albedo and emissivity,
3395!--    whereas paramter 3/4 of water_pars_f are heat conductivities!
3396       IF ( water_type_f%from_file )  THEN
3397!
3398!--       Horizontal surfaces
3399          DO  m = 1, surf_lsm_h%ns
3400             i = surf_lsm_h%i(m)
3401             j = surf_lsm_h%j(m)
3402             
3403             st = water_type_f%var(j,i)
3404             IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
3405                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
3406                   t_soil_h%var_2d(:,m) = water_pars(ind_w_temp,st)
3407                surf_lsm_h%z0(m)     = water_pars(ind_w_z0,st)
3408                surf_lsm_h%z0h(m)    = water_pars(ind_w_z0h,st)
3409                surf_lsm_h%z0q(m)    = water_pars(ind_w_z0h,st)
3410                surf_lsm_h%lambda_surface_s(m) = water_pars(ind_w_lambda_s,st)
3411                surf_lsm_h%lambda_surface_u(m) = water_pars(ind_w_lambda_u,st)             
3412                surf_lsm_h%c_surface(m)        = 0.0_wp
3413                surf_lsm_h%albedo_type(ind_wat_win,m) = INT( water_pars(ind_w_at,st) )
3414                surf_lsm_h%emissivity(ind_wat_win,m)  = water_pars(ind_w_emis,st)
3415               
3416                surf_lsm_h%water_type(m)      = st
3417                surf_lsm_h%water_type_name(m) = water_type_name(st)
3418             ENDIF
3419          ENDDO
3420!
3421!--       Vertical surfaces
3422          DO  l = 0, 3
3423             DO  m = 1, surf_lsm_v(l)%ns
3424                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3425                                                surf_lsm_v(l)%building_covered(m) ) 
3426                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3427                                                surf_lsm_v(l)%building_covered(m) ) 
3428             
3429                st = water_type_f%var(j,i)
3430                IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
3431                   IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  &
3432                      t_soil_v(l)%var_2d(:,m) = water_pars(ind_w_temp,st)
3433                   surf_lsm_v(l)%z0(m)     = water_pars(ind_w_z0,st)
3434                   surf_lsm_v(l)%z0h(m)    = water_pars(ind_w_z0h,st)
3435                   surf_lsm_v(l)%z0q(m)    = water_pars(ind_w_z0h,st)
3436                   surf_lsm_v(l)%lambda_surface_s(m) =                         &
3437                                                   water_pars(ind_w_lambda_s,st)
3438                   surf_lsm_v(l)%lambda_surface_u(m) =                         &
3439                                                   water_pars(ind_w_lambda_u,st)           
3440                   surf_lsm_v(l)%c_surface(m)     = 0.0_wp
3441                   surf_lsm_v(l)%albedo_type(ind_wat_win,m) =                  &
3442                                                  INT( water_pars(ind_w_at,st) )
3443                   surf_lsm_v(l)%emissivity(ind_wat_win,m)  =                  &
3444                                                  water_pars(ind_w_emis,st)
3445                                                 
3446                   surf_lsm_v(l)%water_type(m)      = st
3447                   surf_lsm_v(l)%water_type_name(m) = water_type_name(st)
3448                ENDIF
3449             ENDDO
3450          ENDDO
3451       ENDIF     
3452
3453!
3454!--    Level 3, initialization of water parameters at single (x,y)
3455!--    position via water_pars read from file.
3456       IF ( water_pars_f%from_file )  THEN
3457!
3458!--       Horizontal surfaces
3459          DO  m = 1, surf_lsm_h%ns
3460             i = surf_lsm_h%i(m)
3461             j = surf_lsm_h%j(m)
3462!
3463!--          If surface element is not a water surface and any value in
3464!--          water_pars is given, neglect this information and give an
3465!--          informative message that this value will not be used.   
3466             IF ( .NOT. surf_lsm_h%water_surface(m)  .AND.                     &
3467                   ANY( water_pars_f%pars_xy(:,j,i) /= water_pars_f%fill ) )  THEN
3468                WRITE( message_string, * )                                     &
3469                              'surface element at grid point (j,i) = (',       &
3470                              j, i, ') is not a water surface, ',              &
3471                              'so that information given in ',                 &
3472                              'water_pars at this point is neglected.' 
3473                CALL message( 'land_surface_model_mod', 'PA0645', 0, 0, myid, 6, 0 )
3474             ELSE
3475                IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                   &
3476                     water_pars_f%fill  .AND.                                  &
3477                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
3478                      t_soil_h%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
3479
3480                IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /= water_pars_f%fill ) &
3481                   surf_lsm_h%z0(m)     = water_pars_f%pars_xy(ind_w_z0,j,i)
3482
3483                IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /= water_pars_f%fill )&
3484                THEN
3485                   surf_lsm_h%z0h(m)    = water_pars_f%pars_xy(ind_w_z0h,j,i)
3486                   surf_lsm_h%z0q(m)    = water_pars_f%pars_xy(ind_w_z0h,j,i)
3487                ENDIF
3488                IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /=               &
3489                     water_pars_f%fill )                                       &
3490                   surf_lsm_h%lambda_surface_s(m) =                            &
3491                                        water_pars_f%pars_xy(ind_w_lambda_s,j,i)
3492
3493                IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /=               &
3494                      water_pars_f%fill )                                      &
3495                   surf_lsm_h%lambda_surface_u(m) =                            &
3496                                        water_pars_f%pars_xy(ind_w_lambda_u,j,i)     
3497       
3498                IF ( water_pars_f%pars_xy(ind_w_at,j,i) /=                     &
3499                     water_pars_f%fill )                                       &
3500                   surf_lsm_h%albedo_type(ind_wat_win,m) =                     &
3501                                       INT( water_pars_f%pars_xy(ind_w_at,j,i) )
3502
3503                IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /=                   &
3504                     water_pars_f%fill )                                       &
3505                   surf_lsm_h%emissivity(ind_wat_win,m) =                      &
3506                                          water_pars_f%pars_xy(ind_w_emis,j,i) 
3507             ENDIF
3508          ENDDO
3509!
3510!--       Vertical surfaces
3511          DO  l = 0, 3
3512             DO  m = 1, surf_lsm_v(l)%ns
3513                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3514                                                surf_lsm_v(l)%building_covered(m) ) 
3515                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3516                                                surf_lsm_v(l)%building_covered(m) ) 
3517!
3518!--             If surface element is not a water surface and any value in
3519!--             water_pars is given, neglect this information and give an
3520!--             informative message that this value will not be used.   
3521                IF ( .NOT. surf_lsm_v(l)%water_surface(m)  .AND.               &
3522                      ANY( water_pars_f%pars_xy(:,j,i) /=                      &
3523                      water_pars_f%fill ) )  THEN
3524                   WRITE( message_string, * )                                  &
3525                              'surface element at grid point (j,i) = (',       &
3526                              j, i, ') is not a water surface, ',              &
3527                              'so that information given in ',                 &
3528                              'water_pars at this point is neglected.' 
3529                   CALL message( 'land_surface_model_mod', 'PA0645',           &
3530                                  0, 0, myid, 6, 0 )
3531                ELSE
3532
3533                   IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                &
3534                     water_pars_f%fill  .AND.                                  &
3535                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
3536                      t_soil_v(l)%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
3537
3538                   IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /=                  &
3539                        water_pars_f%fill )                                    &
3540                      surf_lsm_v(l)%z0(m)   = water_pars_f%pars_xy(ind_w_z0,j,i)
3541
3542                   IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /=                 &
3543                       water_pars_f%fill )  THEN
3544                      surf_lsm_v(l)%z0h(m)  = water_pars_f%pars_xy(ind_w_z0h,j,i)
3545                      surf_lsm_v(l)%z0q(m)  = water_pars_f%pars_xy(ind_w_z0h,j,i)
3546                   ENDIF
3547
3548                   IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /=            &
3549                        water_pars_f%fill )                                    &
3550                      surf_lsm_v(l)%lambda_surface_s(m) =                      &
3551                                      water_pars_f%pars_xy(ind_w_lambda_s,j,i)
3552
3553                   IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /=            &
3554                        water_pars_f%fill )                                    &
3555                      surf_lsm_v(l)%lambda_surface_u(m) =                      &
3556                                      water_pars_f%pars_xy(ind_w_lambda_u,j,i)   
3557 
3558                   IF ( water_pars_f%pars_xy(ind_w_at,j,i) /=                  &
3559                        water_pars_f%fill )                                    &
3560                      surf_lsm_v(l)%albedo_type(ind_wat_win,m) =               &
3561                                      INT( water_pars_f%pars_xy(ind_w_at,j,i) )
3562
3563                   IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /=                &
3564                        water_pars_f%fill )                                    &
3565                      surf_lsm_v(l)%emissivity(ind_wat_win,m)  =               &
3566                                      water_pars_f%pars_xy(ind_w_emis,j,i) 
3567                ENDIF
3568             ENDDO
3569          ENDDO
3570
3571       ENDIF
3572!
3573!--    Initialize pavement-type surfaces, level 1
3574       IF ( pavement_type /= 0 )  THEN 
3575
3576!
3577!--       When a pavement_type is used, overwrite a possible setting of
3578!--       the pavement depth as it is already defined by the pavement type
3579          pavement_depth_level = 0
3580
3581          IF ( z0_pavement == 9999999.9_wp )  THEN
3582             z0_pavement  = pavement_pars(ind_p_z0,pavement_type) 
3583          ENDIF
3584
3585          IF ( z0h_pavement == 9999999.9_wp )  THEN
3586             z0h_pavement = pavement_pars(ind_p_z0h,pavement_type)
3587          ENDIF
3588         
3589          IF ( z0q_pavement == 9999999.9_wp )  THEN
3590             z0q_pavement = pavement_pars(ind_p_z0h,pavement_type)
3591          ENDIF
3592
3593          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
3594             pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type)
3595          ENDIF
3596
3597          IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
3598             pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)
3599          ENDIF   
3600   
3601          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
3602             albedo_type = INT(pavement_pars(ind_p_at,pavement_type))       
3603          ENDIF
3604   
3605          IF ( emissivity == 9999999.9_wp )  THEN
3606             emissivity = pavement_pars(ind_p_emis,pavement_type)       
3607          ENDIF
3608
3609!
3610!--       If the depth level of the pavement is not set, determine it from
3611!--       lookup table.
3612          IF ( pavement_depth_level == 0 )  THEN
3613             DO  k = nzb_soil, nzt_soil 
3614                IF ( pavement_subsurface_pars_1(k,pavement_type) == 9999999.9_wp &
3615                .OR. pavement_subsurface_pars_2(k,pavement_type) == 9999999.9_wp)&
3616                THEN
3617                   nzt_pavement = k-1
3618                   EXIT
3619                ENDIF
3620             ENDDO
3621          ELSE
3622             nzt_pavement = pavement_depth_level
3623          ENDIF
3624
3625       ENDIF
3626!
3627!--    Level 1 initialization of pavement type surfaces. Horizontally
3628!--    homogeneous characteristics are assumed
3629       surf_lsm_h%nzt_pavement = pavement_depth_level
3630       DO  m = 1, surf_lsm_h%ns
3631          IF ( surf_lsm_h%pavement_surface(m) )  THEN
3632             surf_lsm_h%nzt_pavement(m)        = nzt_pavement
3633             surf_lsm_h%z0(m)                  = z0_pavement
3634             surf_lsm_h%z0h(m)                 = z0h_pavement
3635             surf_lsm_h%z0q(m)                 = z0q_pavement
3636             surf_lsm_h%lambda_surface_s(m)    = pavement_heat_conduct         &
3637                                                  * ddz_soil(nzb_soil)         &
3638                                                  * 2.0_wp   
3639             surf_lsm_h%lambda_surface_u(m)    = pavement_heat_conduct         &
3640                                                  * ddz_soil(nzb_soil)         &
3641                                                  * 2.0_wp           
3642             surf_lsm_h%c_surface(m)           = pavement_heat_capacity        &
3643                                                        * dz_soil(nzb_soil)    &
3644                                                        * 0.25_wp                                   
3645
3646             surf_lsm_h%albedo_type(ind_pav_green,m) = albedo_type
3647             surf_lsm_h%emissivity(ind_pav_green,m)  = emissivity     
3648             
3649             surf_lsm_h%pavement_type(m)      = pavement_type
3650             surf_lsm_h%pavement_type_name(m) = pavement_type_name(pavement_type)
3651     
3652             IF ( pavement_type /= 0 )  THEN
3653                DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
3654                   surf_lsm_h%lambda_h_def(k,m)    =                           &
3655                                     pavement_subsurface_pars_1(k,pavement_type)                       
3656                   surf_lsm_h%rho_c_total_def(k,m) =                           &
3657                                     pavement_subsurface_pars_2(k,pavement_type) 
3658                ENDDO
3659             ELSE
3660                surf_lsm_h%lambda_h_def(:,m)     = pavement_heat_conduct
3661                surf_lsm_h%rho_c_total_def(:,m)  = pavement_heat_capacity
3662             ENDIF       
3663          ENDIF
3664       ENDDO                               
3665
3666       DO  l = 0, 3
3667          surf_lsm_v(l)%nzt_pavement = pavement_depth_level
3668          DO  m = 1, surf_lsm_v(l)%ns
3669             IF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
3670                surf_lsm_v(l)%nzt_pavement