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

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

Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.

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