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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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