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

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

last changes documented

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