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

Last change on this file since 3868 was 3868, checked in by suehring, 2 years ago

More strict limitation of roughness length when it is in the order of the vertical grid spacing

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