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

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

Bugfix in albedo initialization, caused crashes in rrtmg calls

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