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

Last change on this file since 4495 was 4495, checked in by raasch, 14 months ago

restart data handling with MPI-IO added, first part

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