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

Last change on this file since 3248 was 3248, checked in by sward, 6 years ago

Minor format changes

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