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

Last change on this file since 3620 was 3620, checked in by moh.hefny, 6 years ago

fix date and time and update 3d rad_lw_out

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