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

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

bugfix for sea surfaces with constant_roughness = .F.

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