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

Last change on this file since 2797 was 2797, checked in by suehring, 5 years ago

Output of ground-heat flux at natural- and urban-type surfaces in one output variable; enable restart data of _av variables that belong to both land- and urban-surface model

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