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

Last change on this file since 2765 was 2765, checked in by maronga, 4 years ago

major bugfix in calculation of aerodynamic resistance for vertical surface elements

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