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

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

Several critical bugfixes: initialization of surface temperature at water bodies; consider special case in initialization of vertical natural surface elements

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