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

Last change on this file since 3847 was 3847, checked in by suehring, 4 years ago

Adjust message-call for checks that are especially carried out locally.

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