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

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

bugfixes for spinup mechanism to work with lsm+usm+radiation

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