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

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

Initialization of soil temperature and moisture via dynamic input file only for vegetation and pavement surfaces

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