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

Last change on this file since 4356 was 4356, checked in by suehring, 16 months ago

Bugfix in message calls for local checks; error messages in init_grid slightly revised; bugfix in time_integration (uninitialized emission time index); read_restart_data (change from J.Resler): use dynamic array allocation rather than automatic arrays to avoid problems with stack memory

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