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

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

Move checks for correct dimensions in static input file; improve checks concerning buildings; check whether at least one surface type is set at a natural-type surface element

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