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

Last change on this file since 4188 was 4188, checked in by suehring, 2 years ago

Minor adjustments in error messages and error numbers. Some typos are corrected.

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