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

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

changes documented

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