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

Last change on this file since 3943 was 3943, checked in by maronga, 2 years ago

bugfixes in urban surface model; output of greenz roof transpiration added/corrected; minor formatting improvements

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