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

Last change on this file since 3987 was 3987, checked in by kanani, 2 years ago

clean up location, debug and error messages

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