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

Last change on this file since 3152 was 3152, checked in by suehring, 3 years ago

Further adjustments for surface structure

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