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

Last change on this file since 2724 was 2724, checked in by maronga, 7 years ago

some changes in spinup mechanism and additional check in land surface model

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