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

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

Land-surface model: Revise limitation for soil moisture in case it exceeds its saturation value; Revise initialization of soil moisture and temperature in a nested run in case dynamic input information is available. This case, the soil within the child domains can be initialized separately; As part of this revision, migrate the netcdf input of soil temperature / moisture to this module, as well as the routine to inter/extrapolate soil profiles between different grids.; Plant-canopy: Check if any LAD is prescribed when plant-canopy model is applied, in order to avoid crashes in the radiation transfer model; Surface-layer fluxes: Initialization of Obukhov length also at vertical surfaces (if allocated); Urban-surface model: Add checks to ensure that relative fractions of walls, windowns and green surfaces sum-u to one; Revise message calls dealing with local checks

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