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

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

deleting of deprecated files; headers updated where needed

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