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

Last change on this file since 4110 was 4110, checked in by suehring, 5 years ago

last changes documented

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