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

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

added restart with MPI-IO for reading local arrays

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