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

Last change on this file since 4581 was 4581, checked in by suehring, 12 months ago

mesoscale nesting: omit explicit pressure forcing via geostrophic wind components; chemistry: enable profile output of vertical fluxes; urban-surface: bugfix in initialization in case of cyclic_fill

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