source: palm/trunk/SOURCE/land_surface_model_mod.f90

Last change on this file was 4893, checked in by raasch, 7 months ago

revised output of surface data via MPI-IO for better performance

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