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

Last change on this file since 2729 was 2729, checked in by maronga, 6 years ago

input of deep soil temperature separated from prognostic soil temperature

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