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

Last change on this file since 4441 was 4441, checked in by suehring, 19 months ago

Change order of dimension in surface arrays %frac, %emissivity and %albedo to allow for better vectorization in the radiation interactions; Set back turbulent length scale to 8 x grid spacing in the parametrized mode for the synthetic turbulence generator (was accidentally changed in last commit)

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