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

Last change on this file since 3486 was 3486, checked in by maronga, 3 years ago

bugfix in land surface model for the case case of excessive precipitation/dew formation

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