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

Last change on this file since 4251 was 4251, checked in by maronga, 3 years ago

bugfix in vegetation_type look-up table regarding setting of default albedo_type

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