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

Last change on this file since 3677 was 3677, checked in by moh.hefny, 3 years ago

update rad_lw_out according to radiation model used

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