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

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

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

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