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

Last change on this file since 2921 was 2921, checked in by Giersch, 4 years ago

further inipar parameter has been added to restart data, bugfix in spinup mechanism

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