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

Last change on this file since 4535 was 4535, checked in by raasch, 12 months ago

bugfix for restart data format query

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