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

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

Bugfix initialization of %pt_surface array; Output of surface temperature also for default-type surfaces

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