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

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

radiation: Take external radiation input from root domain dynamic file if no dynamic input file is provided for each nested domain; radiation: Combine MPI_ALLREDUCE calls to reduce latency; land_surface + plant_canopy: Give specific error numbers; land-surface: Adjust error messages for local checks; init_3d_model: Deallocate temporary string array for netcdf-data input since it may be re-used in other modules for different input files

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