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

Last change on this file since 3136 was 3136, checked in by suehring, 3 years ago

Limit roughness for heat and moisture in case it exceeds the model-assumed surface-layer height; Limit surface resistance to positive values.

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