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

Last change on this file since 4312 was 4312, checked in by suehring, 17 months ago

Checks and initialization for relative surface fractions revised

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