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

Last change on this file since 3881 was 3881, checked in by suehring, 2 years ago

Make level 3 initialization of urban-surfaces consistent to input data standard; revise flagging of ground-floor level surfaces at buidlings; bugfixes in level 3 initialization of albedo; bugfix in OpenMP directive; check for zero output timestep in surface output; assure that Obukhov lenght does not become zero

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