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

Last change on this file since 3856 was 3856, checked in by suehring, 4 years ago

Bugfix in lsm_init in case no surface fractions are provided

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