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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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