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

Last change on this file since 2968 was 2968, checked in by suehring, 7 years ago

Bugfixes in initalization of land-surface model as well as timeseries data output in case of elevated model surfaces.

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