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

Last change on this file since 4338 was 4338, checked in by suehring, 16 months ago

To avoid divisions by zero, add security factor in calculation of roughness length over water surfaces.

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