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

Last change on this file since 2735 was 2735, checked in by suehring, 6 years ago

Output of resistance also urban-type surfaces

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