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

Last change on this file since 3014 was 3014, checked in by maronga, 3 years ago

series of bugfixes

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