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

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

Enable initialization of z0q for vegetation, pavement and water surfaces via namelist input

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