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

Last change on this file since 2894 was 2894, checked in by Giersch, 7 years ago

Reading/Writing? data in case of restart runs revised

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