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

Last change on this file since 4281 was 4261, checked in by scharf, 19 months ago

bugfix for rev. 4258

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