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

Last change on this file since 3786 was 3786, checked in by raasch, 4 years ago

unused variables removed, interoperable C datatypes introduced in particle type to avoid compiler warnings

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