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

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

Additional namelist parameter to switch on/off the nesting of chemical species; Enable the input of soil data from dynamic input files independent on atmosphere in order to initialize soil properties in nested child domains from dynamic input; Revise error message number for static/dynamic input; Revise and add checks for static/dynamic input; Bugfix, add netcdf4_parallel directive for collective read operation

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