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

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

Limit soil moisture to its saturation value and give a respective warning rather than an error message. Further, adjust checks for number of given soil temperatures for the nested case where the child might have no dynamic input file.

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