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

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

Introduction of addtional surface variables indicating type and name of the surface elements; Bugfix in LSM

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