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

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

In a nested child domain, distinguish between soil moisture and temperature initialization in case the parent domain is initialized via the dynamic input file; in the offline nesting, add a safety factor for the calculation of bulk Richardson number in order to avoid division by zero which can potentially happen if 3D buildings are located directly at the lateral model boundaries

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