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

Last change on this file since 3597 was 3597, checked in by maronga, 6 years ago

revised calculation of near surface air potential temperature

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