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

Last change on this file since 3146 was 3146, checked in by maronga, 3 years ago

major bugfix in calculation of Obukhov length and subsequent handling of surface information

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