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

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

renamed all Fortran NAMELISTS

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