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

Last change on this file since 3668 was 3668, checked in by maronga, 6 years ago

removed most_methods circular and lookup. added improved version of palm_csd

  • Property svn:keywords set to Id
File size: 327.4 KB
Line 
1!> @file land_surface_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: land_surface_model_mod.f90 3668 2019-01-14 12:49:24Z maronga $
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       ! update the 3d field of rad_lw_ou array
2153       IF ( horizontal ) rad_lw_out(k+k_off,j+j_off,i+i_off) = surf%rad_lw_out(m)
2154
2155       IF ( humidity )  THEN
2156          surf%qsws(m)  = - f_qsws * ( surf%qv1(m) - q_s + dq_s_dt             &
2157                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
2158                            surf_t_surface_p%var_1d(m) )
2159
2160          surf%qsws_veg(m)  = - f_qsws_veg  * ( surf%qv1(m) - q_s              &
2161                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
2162                              * surf_t_surface_p%var_1d(m) )
2163
2164          surf%qsws_soil(m) = - f_qsws_soil * ( surf%qv1(m) - q_s              &
2165                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
2166                              * surf_t_surface_p%var_1d(m) )
2167
2168          surf%qsws_liq(m)  = - f_qsws_liq  * ( 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       ENDIF
2172
2173!
2174!--    Calculate the true surface resistance. ABS is used here to avoid negative
2175!--    values that can occur for very small fluxes due to the artifical addition
2176!--    of 1.0E-20.
2177       IF ( .NOT.  humidity )  THEN
2178          surf%r_s(m) = 1.0E10_wp
2179       ELSE
2180          surf%r_s(m) = ABS(rho_lv / (f_qsws + 1.0E-20_wp) - surf%r_a(m))
2181       ENDIF
2182!
2183!--    Calculate change in liquid water reservoir due to dew fall or
2184!--    evaporation of liquid water
2185       IF ( humidity )  THEN
2186!
2187!--       If precipitation is activated, add rain water to qsws_liq
2188!--       and qsws_soil according the the vegetation coverage.
2189!--       precipitation_rate is given in mm.
2190          IF ( precipitation )  THEN
2191
2192!
2193!--          Add precipitation to liquid water reservoir, if possible.
2194!--          Otherwise, add the water to soil. In case of
2195!--          pavements, the exceeding water amount is explicitly removed
2196!--          (as fictive runoff by drainage systems)
2197             IF ( surf%pavement_surface(m) )  THEN
2198                IF ( surf_m_liq%var_1d(m) < m_liq_max )  THEN
2199                   surf%qsws_liq(m) = surf%qsws_liq(m)                         &
2200                                 + prr(k+k_off,j+j_off,i+i_off)                &
2201                                 * hyrho(k+k_off)                              &
2202                                 * 0.001_wp * rho_l * l_v
2203                ENDIF         
2204             ELSE
2205                IF ( surf_m_liq%var_1d(m) < m_liq_max )  THEN
2206                   surf%qsws_liq(m) = surf%qsws_liq(m)                         &
2207                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
2208                                 * hyrho(k+k_off)                              &
2209                                 * 0.001_wp * rho_l * l_v
2210                   surf%qsws_soil(m) = surf%qsws_soil(m) + ( 1.0_wp -          &
2211                                 surf%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
2212                                 * hyrho(k+k_off)                              &
2213                                 * 0.001_wp * rho_l * l_v                                 
2214                ELSE
2215
2216!--                Add precipitation to bare soil according to the bare soil
2217!--                coverage.
2218                   surf%qsws_soil(m) = surf%qsws_soil(m)                       &
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
2223                ENDIF
2224             ENDIF
2225
2226          ENDIF
2227
2228!
2229!--       If the air is saturated, check the reservoir water level
2230          IF ( surf%qsws(m) < 0.0_wp )  THEN
2231!
2232!--          Check if reservoir is full (avoid values > m_liq_max)
2233!--          In that case, qsws_liq goes to qsws_soil for pervious surfaces. In
2234!--          this case qsws_veg is zero anyway (because c_liq = 1),       
2235!--          so that tend is zero and no further check is needed
2236             IF ( surf_m_liq%var_1d(m) == m_liq_max )  THEN
2237                IF ( .NOT. surf%pavement_surface(m))  THEN
2238                   surf%qsws_soil(m) = surf%qsws_soil(m) + surf%qsws_liq(m)
2239                ENDIF
2240                surf%qsws_liq(m)  = 0.0_wp
2241             ENDIF
2242
2243!
2244!--          In case qsws_veg becomes negative (unphysical behavior),
2245!--          let the water enter the liquid water reservoir as dew on the
2246!--          plant
2247             IF ( surf%qsws_veg(m) < 0.0_wp )  THEN
2248                surf%qsws_liq(m) = surf%qsws_liq(m) + surf%qsws_veg(m)
2249                surf%qsws_veg(m) = 0.0_wp
2250             ENDIF
2251          ENDIF                   
2252 
2253          surf%qsws(m) = surf%qsws(m) / l_v
2254 
2255          tend = - surf%qsws_liq(m) * drho_l_lv
2256          surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d *              &
2257                                        ( tsc(2) * tend +                      &
2258                                          tsc(3) * surf_tm_liq_m%var_1d(m) )
2259!
2260!--       Check if reservoir is overfull -> reduce to maximum
2261!--       (conservation of water is violated here)
2262          surf_m_liq_p%var_1d(m) = MIN( surf_m_liq_p%var_1d(m),m_liq_max )
2263
2264!
2265!--       Check if reservoir is empty (avoid values < 0.0)
2266!--       (conservation of water is violated here)
2267          surf_m_liq_p%var_1d(m) = MAX( surf_m_liq_p%var_1d(m), 0.0_wp )
2268!
2269!--       Calculate m_liq tendencies for the next Runge-Kutta step
2270          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2271             IF ( intermediate_timestep_count == 1 )  THEN
2272                surf_tm_liq_m%var_1d(m) = tend
2273             ELSEIF ( intermediate_timestep_count <                            &
2274                      intermediate_timestep_count_max )  THEN
2275                surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend +                  &
2276                                           5.3125_wp * surf_tm_liq_m%var_1d(m)
2277             ENDIF
2278          ENDIF
2279
2280       ENDIF
2281
2282    ENDDO
2283
2284!
2285!-- Make a logical OR for all processes. Force radiation call if at
2286!-- least one processor reached the threshold change in skin temperature
2287    IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count       &
2288         == intermediate_timestep_count_max-1 )  THEN
2289#if defined( __parallel )
2290       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2291       CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,       &
2292                           1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
2293#else
2294       force_radiation_call = force_radiation_call_l
2295#endif
2296       force_radiation_call_l = .FALSE.
2297    ENDIF
2298
2299!
2300!-- Calculate surface water vapor mixing ratio
2301    IF ( humidity )  THEN
2302       CALL calc_q_surface
2303    ENDIF
2304
2305!
2306!-- Calculate new roughness lengths (for water surfaces only)
2307    IF ( horizontal  .AND.  .NOT. constant_roughness )  CALL calc_z0_water_surface
2308   
2309   
2310    CONTAINS
2311!------------------------------------------------------------------------------!
2312! Description:
2313! ------------
2314!> Calculation of mixing ratio of the skin layer (surface). It is assumend
2315!> that the skin is always saturated.
2316!------------------------------------------------------------------------------!
2317    SUBROUTINE calc_q_surface
2318
2319       IMPLICIT NONE
2320
2321       REAL(wp) ::  e_s           !< saturation water vapor pressure
2322       REAL(wp) ::  q_s           !< saturation mixing ratio
2323       REAL(wp) ::  resistance    !< aerodynamic and soil resistance term
2324
2325       DO  m = 1, surf%ns
2326
2327          i   = surf%i(m)           
2328          j   = surf%j(m)
2329          k   = surf%k(m)
2330!
2331!--       Calculate water vapour pressure at saturation and convert to hPa
2332          e_s = 0.01_wp * magnus( MIN(surf_t_surface_p%var_1d(m), 333.15_wp) )                   
2333
2334!
2335!--       Calculate mixing ratio at saturation
2336          q_s = rd_d_rv * e_s / ( surface_pressure - e_s )
2337
2338          resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) + 1E-5_wp )
2339
2340!
2341!--       Calculate mixing ratio at surface
2342          IF ( bulk_cloud_model )  THEN
2343             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
2344                                        ( 1.0_wp - resistance ) *              &
2345                                        ( q(k,j,i) - ql(k,j,i) )
2346          ELSE
2347             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
2348                                        ( 1.0_wp - resistance ) *              &
2349                                          q(k,j,i)
2350          ENDIF
2351         
2352          surf%q_surface(m) = q(k+k_off,j+j_off,i+i_off)
2353!
2354!--       Update virtual potential temperature
2355          surf%vpt_surface(m) = surf%pt_surface(m) *                           &
2356                                  ( 1.0_wp + 0.61_wp * surf%q_surface(m) )
2357
2358       
2359                     
2360       ENDDO
2361 
2362    END SUBROUTINE calc_q_surface
2363       
2364 END SUBROUTINE lsm_energy_balance
2365   
2366   
2367
2368!------------------------------------------------------------------------------!
2369! Description:
2370! ------------
2371!> Header output for land surface model
2372!------------------------------------------------------------------------------!
2373    SUBROUTINE lsm_header ( io )
2374
2375
2376       IMPLICIT NONE
2377
2378       CHARACTER (LEN=86) ::  t_soil_chr          !< String for soil temperature profile
2379       CHARACTER (LEN=86) ::  roots_chr           !< String for root profile
2380       CHARACTER (LEN=86) ::  vertical_index_chr  !< String for the vertical index
2381       CHARACTER (LEN=86) ::  m_soil_chr          !< String for soil moisture
2382       CHARACTER (LEN=86) ::  soil_depth_chr      !< String for soil depth
2383       CHARACTER (LEN=10) ::  coor_chr            !< Temporary string
2384   
2385       INTEGER(iwp) ::  i                         !< Loop index over soil layers
2386 
2387       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
2388 
2389       t_soil_chr = ''
2390       m_soil_chr    = ''
2391       soil_depth_chr  = '' 
2392       roots_chr        = '' 
2393       vertical_index_chr   = ''
2394
2395       i = 1
2396       DO i = nzb_soil, nzt_soil
2397          WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i)
2398          t_soil_chr = TRIM( t_soil_chr ) // ' ' // TRIM( coor_chr )
2399
2400          WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i)
2401          m_soil_chr = TRIM( m_soil_chr ) // ' ' // TRIM( coor_chr )
2402
2403          WRITE (coor_chr,'(F10.2,7X)')  - zs(i)
2404          soil_depth_chr = TRIM( soil_depth_chr ) // ' '  // TRIM( coor_chr )
2405
2406          WRITE (coor_chr,'(F10.2,7X)')  root_fraction(i)
2407          roots_chr = TRIM( roots_chr ) // ' '  // TRIM( coor_chr )
2408
2409          WRITE (coor_chr,'(I10,7X)')  i
2410          vertical_index_chr = TRIM( vertical_index_chr ) // ' '  //           &
2411                               TRIM( coor_chr )
2412       ENDDO
2413
2414!
2415!--    Write land surface model header
2416       WRITE( io,  1 )
2417       IF ( conserve_water_content )  THEN
2418          WRITE( io, 2 )
2419       ELSE
2420          WRITE( io, 3 )
2421       ENDIF
2422
2423       IF ( vegetation_type_f%from_file )  THEN
2424          WRITE( io, 5 )
2425       ELSE
2426          WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ),        &
2427                         TRIM (soil_type_name(soil_type) )
2428       ENDIF
2429       WRITE( io, 6 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
2430                        TRIM( m_soil_chr ), TRIM( roots_chr ),                 &
2431                        TRIM( vertical_index_chr )
2432
24331   FORMAT (//' Land surface model information:'/                              &
2434              ' ------------------------------'/)
24352   FORMAT ('    --> Soil bottom is closed (water content is conserved',       &
2436            ', default)')
24373   FORMAT ('    --> Soil bottom is open (water content is not conserved)')         
24384   FORMAT ('    --> Land surface type  : ',A,/                                &
2439            '    --> Soil porosity type : ',A)
24405   FORMAT ('    --> Land surface type  : read from file' /                    &
2441            '    --> Soil porosity type : read from file' )
24426   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
2443            '       Height:        ',A,'  m'/                                  &
2444            '       Temperature:   ',A,'  K'/                                  &
2445            '       Moisture:      ',A,'  m**3/m**3'/                          &
2446            '       Root fraction: ',A,'  '/                                   &
2447            '       Grid point:    ',A)
2448
2449
2450    END SUBROUTINE lsm_header
2451
2452
2453!------------------------------------------------------------------------------!
2454! Description:
2455! ------------
2456!> Initialization of the land surface model
2457!------------------------------------------------------------------------------!
2458    SUBROUTINE lsm_init
2459
2460       USE control_parameters,                                                 &
2461           ONLY:  message_string
2462
2463       USE indices,                                                            &
2464           ONLY:  nx, ny, topo_min_level
2465
2466       USE pmc_interface,                                                      &
2467           ONLY:  nested_run
2468   
2469       IMPLICIT NONE
2470
2471       LOGICAL      ::  init_soil_from_parent   !< flag controlling initialization of soil in child domains
2472
2473       INTEGER(iwp) ::  i                       !< running index
2474       INTEGER(iwp) ::  i_off                   !< index offset of surface element, seen from atmospheric grid point
2475       INTEGER(iwp) ::  j                       !< running index
2476       INTEGER(iwp) ::  j_off                   !< index offset of surface element, seen from atmospheric grid point
2477       INTEGER(iwp) ::  k                       !< running index
2478       INTEGER(iwp) ::  kn                      !< running index
2479       INTEGER(iwp) ::  ko                      !< running index
2480       INTEGER(iwp) ::  kroot                   !< running index
2481       INTEGER(iwp) ::  kzs                     !< running index
2482       INTEGER(iwp) ::  l                       !< running index surface facing
2483       INTEGER(iwp) ::  m                       !< running index
2484       INTEGER(iwp) ::  st                      !< soil-type index
2485       INTEGER(iwp) ::  n_soil_layers_total     !< temperature variable, stores the total number of soil layers + 4
2486
2487       REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
2488       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pr_soil_init !< temporary array used for averaging soil profiles
2489!
2490!--    If no cloud physics is used, rho_surface has not been calculated before
2491       IF (  .NOT.  bulk_cloud_model  .AND.  .NOT.  cloud_droplets )  THEN
2492          CALL calc_mean_profile( pt, 4 )
2493          rho_surface = hyp(nzb) / ( r_d * hom(topo_min_level+1,1,4,0) * exner(nzb) )
2494       ENDIF
2495
2496!
2497!--    Calculate frequently used parameters
2498       rho_cp    = c_p * rho_surface
2499       rho_lv    = rho_surface * l_v
2500       drho_l_lv = 1.0_wp / (rho_l * l_v)
2501
2502!
2503!--    Set initial values for prognostic quantities
2504!--    Horizontal surfaces
2505       tt_surface_h_m%var_1d = 0.0_wp
2506       tt_soil_h_m%var_2d    = 0.0_wp
2507       tm_soil_h_m%var_2d    = 0.0_wp
2508       tm_liq_h_m%var_1d     = 0.0_wp
2509       surf_lsm_h%c_liq      = 0.0_wp
2510
2511       surf_lsm_h%ghf = 0.0_wp
2512
2513       surf_lsm_h%qsws_liq  = 0.0_wp
2514       surf_lsm_h%qsws_soil = 0.0_wp
2515       surf_lsm_h%qsws_veg  = 0.0_wp
2516
2517       surf_lsm_h%r_a        = 50.0_wp
2518       surf_lsm_h%r_s        = 50.0_wp
2519       surf_lsm_h%r_canopy   = 0.0_wp
2520       surf_lsm_h%r_soil     = 0.0_wp
2521!
2522!--    Do the same for vertical surfaces
2523       DO  l = 0, 3
2524          tt_surface_v_m(l)%var_1d = 0.0_wp
2525          tt_soil_v_m(l)%var_2d    = 0.0_wp
2526          tm_soil_v_m(l)%var_2d    = 0.0_wp
2527          tm_liq_v_m(l)%var_1d     = 0.0_wp
2528          surf_lsm_v(l)%c_liq      = 0.0_wp
2529
2530          surf_lsm_v(l)%ghf = 0.0_wp
2531
2532          surf_lsm_v(l)%qsws_liq  = 0.0_wp
2533          surf_lsm_v(l)%qsws_soil = 0.0_wp
2534          surf_lsm_v(l)%qsws_veg  = 0.0_wp
2535
2536          surf_lsm_v(l)%r_a        = 50.0_wp
2537          surf_lsm_v(l)%r_s        = 50.0_wp
2538          surf_lsm_v(l)%r_canopy   = 0.0_wp
2539          surf_lsm_v(l)%r_soil     = 0.0_wp
2540       ENDDO
2541
2542!
2543!--    Set initial values for prognostic soil quantities
2544       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2545          t_soil_h%var_2d = 0.0_wp
2546          m_soil_h%var_2d = 0.0_wp
2547          m_liq_h%var_1d  = 0.0_wp
2548
2549          DO  l = 0, 3
2550             t_soil_v(l)%var_2d = 0.0_wp
2551             m_soil_v(l)%var_2d = 0.0_wp
2552             m_liq_v(l)%var_1d  = 0.0_wp
2553          ENDDO
2554       ENDIF
2555!
2556!--    Allocate 3D soil model arrays
2557!--    First, for horizontal surfaces
2558       ALLOCATE ( surf_lsm_h%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
2559       ALLOCATE ( surf_lsm_h%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2560       ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
2561       ALLOCATE ( surf_lsm_h%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns))
2562       ALLOCATE ( surf_lsm_h%l_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
2563       ALLOCATE ( surf_lsm_h%m_fc(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
2564       ALLOCATE ( surf_lsm_h%m_res(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
2565       ALLOCATE ( surf_lsm_h%m_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
2566       ALLOCATE ( surf_lsm_h%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_h%ns)      )
2567       ALLOCATE ( surf_lsm_h%n_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
2568       ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2569       ALLOCATE ( surf_lsm_h%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2570       ALLOCATE ( surf_lsm_h%root_fr(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
2571   
2572       surf_lsm_h%lambda_h     = 0.0_wp
2573!
2574!--    If required, allocate humidity-related variables for the soil model
2575       IF ( humidity )  THEN
2576          ALLOCATE ( surf_lsm_h%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
2577          ALLOCATE ( surf_lsm_h%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  ) 
2578
2579          surf_lsm_h%lambda_w = 0.0_wp 
2580       ENDIF
2581!
2582!--    For vertical surfaces
2583       DO  l = 0, 3
2584          ALLOCATE ( surf_lsm_v(l)%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
2585          ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2586          ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
2587          ALLOCATE ( surf_lsm_v(l)%lambda_h_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns))
2588          ALLOCATE ( surf_lsm_v(l)%l_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
2589          ALLOCATE ( surf_lsm_v(l)%m_fc(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
2590          ALLOCATE ( surf_lsm_v(l)%m_res(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
2591          ALLOCATE ( surf_lsm_v(l)%m_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
2592          ALLOCATE ( surf_lsm_v(l)%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)      )
2593          ALLOCATE ( surf_lsm_v(l)%n_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
2594          ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) ) 
2595          ALLOCATE ( surf_lsm_v(l)%rho_c_total_def(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) ) 
2596          ALLOCATE ( surf_lsm_v(l)%root_fr(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )
2597
2598          surf_lsm_v(l)%lambda_h     = 0.0_wp 
2599         
2600!
2601!--       If required, allocate humidity-related variables for the soil model
2602          IF ( humidity )  THEN
2603             ALLOCATE ( surf_lsm_v(l)%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
2604             ALLOCATE ( surf_lsm_v(l)%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  ) 
2605
2606             surf_lsm_v(l)%lambda_w = 0.0_wp 
2607          ENDIF     
2608       ENDDO
2609!
2610!--    Allocate albedo type and emissivity for vegetation, water and pavement
2611!--    fraction.
2612!--    Set default values at each surface element.
2613       ALLOCATE ( surf_lsm_h%albedo_type(0:2,1:surf_lsm_h%ns) )
2614       ALLOCATE ( surf_lsm_h%emissivity(0:2,1:surf_lsm_h%ns) )
2615       surf_lsm_h%albedo_type = 0     
2616       surf_lsm_h%emissivity  = emissivity
2617       DO  l = 0, 3
2618          ALLOCATE ( surf_lsm_v(l)%albedo_type(0:2,1:surf_lsm_v(l)%ns) )
2619          ALLOCATE ( surf_lsm_v(l)%emissivity(0:2,1:surf_lsm_v(l)%ns)  )
2620          surf_lsm_v(l)%albedo_type = 0
2621          surf_lsm_v(l)%emissivity  = emissivity
2622       ENDDO
2623!
2624!--    Allocate arrays for relative surface fraction.
2625!--    0 - vegetation fraction, 2 - water fraction, 1 - pavement fraction
2626       ALLOCATE( surf_lsm_h%frac(0:2,1:surf_lsm_h%ns) )
2627       surf_lsm_h%frac = 0.0_wp
2628       DO  l = 0, 3
2629          ALLOCATE( surf_lsm_v(l)%frac(0:2,1:surf_lsm_v(l)%ns) )
2630          surf_lsm_v(l)%frac = 0.0_wp
2631       ENDDO
2632!
2633!--    For vertical walls only - allocate special flag indicating if any building is on
2634!--    top of any natural surfaces. Used for initialization only.
2635       DO  l = 0, 3
2636          ALLOCATE( surf_lsm_v(l)%building_covered(1:surf_lsm_v(l)%ns) )
2637       ENDDO
2638!
2639!--    Allocate arrays for the respective types and their names on the surface
2640!--    elements. This will be required to treat deposition of chemical species.
2641       ALLOCATE( surf_lsm_h%pavement_type(1:surf_lsm_h%ns)   )
2642       ALLOCATE( surf_lsm_h%vegetation_type(1:surf_lsm_h%ns) )
2643       ALLOCATE( surf_lsm_h%water_type(1:surf_lsm_h%ns)      )
2644       
2645       surf_lsm_h%pavement_type   = 0
2646       surf_lsm_h%vegetation_type = 0
2647       surf_lsm_h%water_type      = 0
2648       
2649       ALLOCATE( surf_lsm_h%pavement_type_name(1:surf_lsm_h%ns)   )
2650       ALLOCATE( surf_lsm_h%vegetation_type_name(1:surf_lsm_h%ns) )
2651       ALLOCATE( surf_lsm_h%water_type_name(1:surf_lsm_h%ns)      )
2652       
2653       surf_lsm_h%pavement_type_name   = 'none'
2654       surf_lsm_h%vegetation_type_name = 'none'
2655       surf_lsm_h%water_type_name      = 'none'
2656       
2657       DO  l = 0, 3
2658          ALLOCATE( surf_lsm_v(l)%pavement_type(1:surf_lsm_v(l)%ns)   )
2659          ALLOCATE( surf_lsm_v(l)%vegetation_type(1:surf_lsm_v(l)%ns) )
2660          ALLOCATE( surf_lsm_v(l)%water_type(1:surf_lsm_v(l)%ns)      )
2661         
2662          surf_lsm_v(l)%pavement_type   = 0
2663          surf_lsm_v(l)%vegetation_type = 0
2664          surf_lsm_v(l)%water_type      = 0
2665       
2666          ALLOCATE( surf_lsm_v(l)%pavement_type_name(1:surf_lsm_v(l)%ns)   )
2667          ALLOCATE( surf_lsm_v(l)%vegetation_type_name(1:surf_lsm_v(l)%ns) )
2668          ALLOCATE( surf_lsm_v(l)%water_type_name(1:surf_lsm_v(l)%ns)      )
2669       
2670          surf_lsm_v(l)%pavement_type_name   = 'none'
2671          surf_lsm_v(l)%vegetation_type_name = 'none'
2672          surf_lsm_v(l)%water_type_name      = 'none'       
2673       ENDDO
2674       
2675!
2676!--    Set flag parameter for the prescribed surface type depending on user
2677!--    input. Set surface fraction to 1 for the respective type.
2678       SELECT CASE ( TRIM( surface_type ) )
2679         
2680          CASE ( 'vegetation' )
2681         
2682             surf_lsm_h%vegetation_surface = .TRUE.
2683             surf_lsm_h%frac(ind_veg_wall,:) = 1.0_wp
2684             DO  l = 0, 3
2685                surf_lsm_v(l)%vegetation_surface = .TRUE.
2686                surf_lsm_v(l)%frac(ind_veg_wall,:) = 1.0_wp
2687             ENDDO
2688   
2689          CASE ( 'water' )
2690             
2691             surf_lsm_h%water_surface = .TRUE.
2692             surf_lsm_h%frac(ind_wat_win,:) = 1.0_wp
2693!
2694!--          Note, vertical water surface does not really make sense.
2695             DO  l = 0, 3 
2696                surf_lsm_v(l)%water_surface   = .TRUE.
2697                surf_lsm_v(l)%frac(ind_wat_win,:) = 1.0_wp
2698             ENDDO
2699
2700          CASE ( 'pavement' )
2701             
2702             surf_lsm_h%pavement_surface = .TRUE.
2703                surf_lsm_h%frac(ind_pav_green,:) = 1.0_wp
2704             DO  l = 0, 3
2705                surf_lsm_v(l)%pavement_surface   = .TRUE.
2706                surf_lsm_v(l)%frac(ind_pav_green,:) = 1.0_wp
2707             ENDDO
2708
2709          CASE ( 'netcdf' )
2710
2711             DO  m = 1, surf_lsm_h%ns
2712                i = surf_lsm_h%i(m)
2713                j = surf_lsm_h%j(m)
2714                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )    &
2715                   surf_lsm_h%vegetation_surface(m) = .TRUE.
2716                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )      &
2717                   surf_lsm_h%pavement_surface(m) = .TRUE.
2718                IF ( water_type_f%var(j,i)      /= water_type_f%fill )         &
2719                   surf_lsm_h%water_surface(m) = .TRUE.
2720             ENDDO
2721!
2722!--          For vertical surfaces some special checks and treatment are
2723!--          required for correct initialization.
2724             DO  l = 0, 3
2725                DO  m = 1, surf_lsm_v(l)%ns
2726!
2727!--                Only for vertical surfaces. Check if at the grid point where
2728!--                the wall is defined (i+ioff, j+joff) is any building.
2729!--                This case, no natural surfaces properties will be defined at
2730!--                at this grid point, leading to problems in the initialization.
2731!--                To overcome this, define a special flag which
2732!--                indicates that a building is defined at the wall grid point 
2733!--                and take the surface properties from the adjoining grid
2734!--                point, i.e. without offset values.
2735!--                Further, there can occur a special case where elevation
2736!--                changes are larger than building heights. This case, (j,i)
2737!--                and (j+joff,i+ioff) grid points may be both covered by
2738!--                buildings, but vertical, but vertically natural walls may
2739!--                be located between the buildings. This case, it is not
2740!--                guaranteed that information about natural surface types is
2741!--                given, neither at (j,i) nor at (j+joff,i+ioff), again leading
2742!--                to non-initialized surface properties.
2743                   surf_lsm_v(l)%building_covered(m) = .FALSE.
2744!
2745!--                Wall grid point is building-covered. This case, set
2746!--                flag indicating that surface properties are initialized
2747!--                from neighboring reference grid point, which is not
2748!--                building_covered.
2749                   IF ( building_type_f%from_file )  THEN
2750                      i = surf_lsm_v(l)%i(m)
2751                      j = surf_lsm_v(l)%j(m)
2752                      IF ( building_type_f%var(j+surf_lsm_v(l)%joff,           &
2753                                               i+surf_lsm_v(l)%ioff) /=        &
2754                           building_type_f%fill )                              &
2755                         surf_lsm_v(l)%building_covered(m) = .TRUE.
2756!
2757!--                   Wall grid point as well as neighboring reference grid
2758!--                   point are both building-covered. This case, surface
2759!--                   properties are not necessarily defined (not covered by
2760!--                   checks for static input file) at this surface. Hence,
2761!--                   initialize surface properties by simply setting
2762!--                   vegetation_type_f to bare-soil bulk parametrization.
2763!--                   soil_type_f as well as surface_fractions_f will be set
2764!--                   also.                     
2765                      IF ( building_type_f%var(j+surf_lsm_v(l)%joff,           &
2766                                               i+surf_lsm_v(l)%ioff) /=        &
2767                           building_type_f%fill  .AND.                         &
2768                           building_type_f%var(j,i) /= building_type_f%fill )  &
2769                      THEN
2770                         vegetation_type_f%var(j,i)                 = 1 ! bare soil
2771                         soil_type_f%var_2d(j,i)                    = 1 
2772                         surface_fraction_f%frac(ind_veg_wall,j,i)  = 1.0_wp
2773                         surface_fraction_f%frac(ind_pav_green,j,i) = 0.0_wp
2774                         surface_fraction_f%frac(ind_wat_win,j,i)   = 0.0_wp
2775                      ENDIF
2776                     
2777                   ENDIF
2778!
2779!--                Normally proceed with setting surface types.
2780                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2781                                            surf_lsm_v(l)%building_covered(m) )
2782                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2783                                            surf_lsm_v(l)%building_covered(m) )
2784                   IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
2785                      surf_lsm_v(l)%vegetation_surface(m) = .TRUE.
2786                   IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )   &
2787                      surf_lsm_v(l)%pavement_surface(m) = .TRUE.
2788                   IF ( water_type_f%var(j,i)      /= water_type_f%fill )      &
2789                      surf_lsm_v(l)%water_surface(m) = .TRUE.
2790!
2791!--                Check if surface element has the flag %building_covered_all,
2792!--                this case, set vegetation_type appropriate
2793                ENDDO
2794             ENDDO
2795
2796       END SELECT
2797!
2798!--    In case of netcdf input file, further initialize surface fractions.
2799!--    At the moment only 1 surface is given at a location, so that the fraction
2800!--    is either 0 or 1. This will be revised later. If surface fraction
2801!--    is not given in static input file, relative fractions will be derived
2802!--    from given surface type. In this case, only 1 type is given at a certain
2803!--    location (already checked). 
2804       IF ( input_pids_static  .AND.  surface_fraction_f%from_file )  THEN
2805          DO  m = 1, surf_lsm_h%ns
2806             i = surf_lsm_h%i(m)
2807             j = surf_lsm_h%j(m)
2808!
2809!--          0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction             
2810             surf_lsm_h%frac(ind_veg_wall,m)  =                                &
2811                                    surface_fraction_f%frac(ind_veg_wall,j,i)         
2812             surf_lsm_h%frac(ind_pav_green,m) =                                &
2813                                    surface_fraction_f%frac(ind_pav_green,j,i)       
2814             surf_lsm_h%frac(ind_wat_win,m)   =                                &
2815                                    surface_fraction_f%frac(ind_wat_win,j,i)
2816
2817          ENDDO
2818          DO  l = 0, 3
2819             DO  m = 1, surf_lsm_v(l)%ns
2820                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
2821                                                surf_lsm_v(l)%building_covered(m) ) 
2822                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
2823                                                surf_lsm_v(l)%building_covered(m) ) 
2824!
2825!--             0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction       
2826                surf_lsm_v(l)%frac(ind_veg_wall,m)  =                          &
2827                                    surface_fraction_f%frac(ind_veg_wall,j,i)         
2828                surf_lsm_v(l)%frac(ind_pav_green,m) =                          &
2829                                    surface_fraction_f%frac(ind_pav_green,j,i)       
2830                surf_lsm_v(l)%frac(ind_wat_win,m)   =                          &
2831                                    surface_fraction_f%frac(ind_wat_win,j,i)
2832
2833             ENDDO
2834          ENDDO
2835       ELSEIF ( input_pids_static  .AND.  .NOT. surface_fraction_f%from_file ) &
2836       THEN
2837          DO  m = 1, surf_lsm_h%ns
2838             i = surf_lsm_h%i(m)
2839             j = surf_lsm_h%j(m)
2840
2841             IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )       &       
2842                surf_lsm_h%frac(ind_veg_wall,m)  = 1.0_wp
2843             IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill   )       &       
2844                surf_lsm_h%frac(ind_pav_green,m) = 1.0_wp 
2845             IF ( water_type_f%var(j,i)      /= water_type_f%fill      )       &       
2846                surf_lsm_h%frac(ind_wat_win,m)   = 1.0_wp       
2847          ENDDO
2848          DO  l = 0, 3
2849             DO  m = 1, surf_lsm_v(l)%ns
2850                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
2851                                                surf_lsm_v(l)%building_covered(m) ) 
2852                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
2853                                                surf_lsm_v(l)%building_covered(m) ) 
2854     
2855                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )    &       
2856                   surf_lsm_v(l)%frac(ind_veg_wall,m)  = 1.0_wp
2857                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill   )    &       
2858                   surf_lsm_v(l)%frac(ind_pav_green,m) = 1.0_wp 
2859                IF ( water_type_f%var(j,i)      /= water_type_f%fill      )    &       
2860                   surf_lsm_v(l)%frac(ind_wat_win,m)   = 1.0_wp     
2861             ENDDO
2862          ENDDO
2863       ENDIF
2864!
2865!--    Level 1, initialization of soil parameters.
2866!--    It is possible to overwrite each parameter by setting the respecticy
2867!--    NAMELIST variable to a value /= 9999999.9.
2868       IF ( soil_type /= 0 )  THEN 
2869 
2870          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
2871             alpha_vangenuchten = soil_pars(0,soil_type)
2872          ENDIF
2873
2874          IF ( l_vangenuchten == 9999999.9_wp )  THEN
2875             l_vangenuchten = soil_pars(1,soil_type)
2876          ENDIF
2877
2878          IF ( n_vangenuchten == 9999999.9_wp )  THEN
2879             n_vangenuchten = soil_pars(2,soil_type)           
2880          ENDIF
2881
2882          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
2883             hydraulic_conductivity = soil_pars(3,soil_type)           
2884          ENDIF
2885
2886          IF ( saturation_moisture == 9999999.9_wp )  THEN
2887             saturation_moisture = soil_pars(4,soil_type)           
2888          ENDIF
2889
2890          IF ( field_capacity == 9999999.9_wp )  THEN
2891             field_capacity = soil_pars(5,soil_type)           
2892          ENDIF
2893
2894          IF ( wilting_point == 9999999.9_wp )  THEN
2895             wilting_point = soil_pars(6,soil_type)           
2896          ENDIF
2897
2898          IF ( residual_moisture == 9999999.9_wp )  THEN
2899             residual_moisture = soil_pars(7,soil_type)       
2900          ENDIF
2901
2902       ENDIF
2903!
2904!--    Map values to the respective 2D/3D arrays
2905!--    Horizontal surfaces
2906       surf_lsm_h%alpha_vg      = alpha_vangenuchten
2907       surf_lsm_h%l_vg          = l_vangenuchten
2908       surf_lsm_h%n_vg          = n_vangenuchten 
2909       surf_lsm_h%gamma_w_sat   = hydraulic_conductivity
2910       surf_lsm_h%m_sat         = saturation_moisture
2911       surf_lsm_h%m_fc          = field_capacity
2912       surf_lsm_h%m_wilt        = wilting_point
2913       surf_lsm_h%m_res         = residual_moisture
2914       surf_lsm_h%r_soil_min    = min_soil_resistance
2915!
2916!--    Vertical surfaces
2917       DO  l = 0, 3
2918          surf_lsm_v(l)%alpha_vg      = alpha_vangenuchten
2919          surf_lsm_v(l)%l_vg          = l_vangenuchten
2920          surf_lsm_v(l)%n_vg          = n_vangenuchten 
2921          surf_lsm_v(l)%gamma_w_sat   = hydraulic_conductivity
2922          surf_lsm_v(l)%m_sat         = saturation_moisture
2923          surf_lsm_v(l)%m_fc          = field_capacity
2924          surf_lsm_v(l)%m_wilt        = wilting_point
2925          surf_lsm_v(l)%m_res         = residual_moisture
2926          surf_lsm_v(l)%r_soil_min    = min_soil_resistance
2927       ENDDO
2928!
2929!--    Level 2, initialization of soil parameters via soil_type read from file.
2930!--    Soil parameters are initialized for each (y,x)-grid point
2931!--    individually using default paramter settings according to the given
2932!--    soil type.
2933       IF ( soil_type_f%from_file )  THEN
2934!
2935!--       Level of detail = 1, i.e. a homogeneous soil distribution along the
2936!--       vertical dimension is assumed.
2937          IF ( soil_type_f%lod == 1 )  THEN
2938!
2939!--          Horizontal surfaces
2940             DO  m = 1, surf_lsm_h%ns
2941                i = surf_lsm_h%i(m)
2942                j = surf_lsm_h%j(m)
2943             
2944                st = soil_type_f%var_2d(j,i)
2945                IF ( st /= soil_type_f%fill )  THEN
2946                   surf_lsm_h%alpha_vg(:,m)    = soil_pars(0,st)
2947                   surf_lsm_h%l_vg(:,m)        = soil_pars(1,st)
2948                   surf_lsm_h%n_vg(:,m)        = soil_pars(2,st)
2949                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars(3,st)
2950                   surf_lsm_h%m_sat(:,m)       = soil_pars(4,st)
2951                   surf_lsm_h%m_fc(:,m)        = soil_pars(5,st)
2952                   surf_lsm_h%m_wilt(:,m)      = soil_pars(6,st)
2953                   surf_lsm_h%m_res(:,m)       = soil_pars(7,st)
2954                ENDIF
2955             ENDDO
2956!
2957!--          Vertical surfaces ( assumes the soil type given at respective (x,y)
2958             DO  l = 0, 3
2959                DO  m = 1, surf_lsm_v(l)%ns
2960                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
2961                                                   surf_lsm_v(l)%building_covered(m) ) 
2962                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
2963                                                   surf_lsm_v(l)%building_covered(m) ) 
2964
2965                   st = soil_type_f%var_2d(j,i)
2966                   IF ( st /= soil_type_f%fill )  THEN
2967                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars(0,st)
2968                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars(1,st)
2969                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars(2,st)
2970                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars(3,st)
2971                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars(4,st)
2972                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars(5,st)
2973                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars(6,st)
2974                      surf_lsm_v(l)%m_res(:,m)       = soil_pars(7,st)
2975                   ENDIF
2976                ENDDO
2977             ENDDO
2978!
2979!--       Level of detail = 2, i.e. soil type and thus the soil parameters
2980!--       can be heterogeneous along the vertical dimension.
2981          ELSE
2982!
2983!--          Horizontal surfaces
2984             DO  m = 1, surf_lsm_h%ns
2985                i = surf_lsm_h%i(m)
2986                j = surf_lsm_h%j(m)
2987             
2988                DO  k = nzb_soil, nzt_soil
2989                   st = soil_type_f%var_3d(k,j,i)
2990                   IF ( st /= soil_type_f%fill )  THEN
2991                      surf_lsm_h%alpha_vg(k,m)    = soil_pars(0,st)
2992                      surf_lsm_h%l_vg(k,m)        = soil_pars(1,st)
2993                      surf_lsm_h%n_vg(k,m)        = soil_pars(2,st)
2994                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars(3,st)
2995                      surf_lsm_h%m_sat(k,m)       = soil_pars(4,st)
2996                      surf_lsm_h%m_fc(k,m)        = soil_pars(5,st)
2997                      surf_lsm_h%m_wilt(k,m)      = soil_pars(6,st)
2998                      surf_lsm_h%m_res(k,m)       = soil_pars(7,st)
2999                   ENDIF
3000                ENDDO
3001             ENDDO
3002!
3003!--          Vertical surfaces ( assumes the soil type given at respective (x,y)
3004             DO  l = 0, 3
3005                DO  m = 1, surf_lsm_v(l)%ns
3006                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
3007                                                   surf_lsm_v(l)%building_covered(m) ) 
3008                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
3009                                                   surf_lsm_v(l)%building_covered(m) ) 
3010
3011                   DO  k = nzb_soil, nzt_soil
3012                      st = soil_type_f%var_3d(k,j,i)
3013                      IF ( st /= soil_type_f%fill )  THEN
3014                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars(0,st)
3015                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars(1,st)
3016                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars(2,st)
3017                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars(3,st)
3018                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars(4,st)
3019                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars(5,st)
3020                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars(6,st)
3021                         surf_lsm_v(l)%m_res(k,m)       = soil_pars(7,st)
3022                      ENDIF
3023                   ENDDO
3024                ENDDO
3025             ENDDO
3026          ENDIF
3027       ENDIF
3028!
3029!--    Level 3, initialization of single soil parameters at single z,x,y
3030!--    position via soil_pars read from file.
3031       IF ( soil_pars_f%from_file )  THEN
3032!
3033!--       Level of detail = 1, i.e. a homogeneous vertical distribution of soil
3034!--       parameters is assumed.
3035!--       Horizontal surfaces
3036          IF ( soil_pars_f%lod == 1 )  THEN
3037!
3038!--          Horizontal surfaces
3039             DO  m = 1, surf_lsm_h%ns
3040                i = surf_lsm_h%i(m)
3041                j = surf_lsm_h%j(m)
3042
3043                IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )              &
3044                   surf_lsm_h%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
3045                IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )              &
3046                   surf_lsm_h%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
3047                IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )              &
3048                   surf_lsm_h%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
3049                IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )              &
3050                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
3051                IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )              &
3052                   surf_lsm_h%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
3053                IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )              &
3054                   surf_lsm_h%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
3055                IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )              &
3056                   surf_lsm_h%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
3057                IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )              &
3058                   surf_lsm_h%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)
3059
3060             ENDDO
3061!
3062!--          Vertical surfaces
3063             DO  l = 0, 3
3064                DO  m = 1, surf_lsm_v(l)%ns
3065                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
3066                                                   surf_lsm_v(l)%building_covered(m) ) 
3067                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
3068                                                   surf_lsm_v(l)%building_covered(m) ) 
3069
3070                   IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )           &
3071                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
3072                   IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )           &
3073                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
3074                   IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )           &
3075                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
3076                   IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )           &
3077                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
3078                   IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )           &
3079                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
3080                   IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )           &
3081                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
3082                   IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )           &
3083                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
3084                   IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )           &
3085                      surf_lsm_v(l)%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)
3086
3087                ENDDO
3088             ENDDO
3089!
3090!--       Level of detail = 2, i.e. soil parameters can be set at each soil
3091!--       layer individually.
3092          ELSE
3093!
3094!--          Horizontal surfaces
3095             DO  m = 1, surf_lsm_h%ns
3096                i = surf_lsm_h%i(m)
3097                j = surf_lsm_h%j(m)
3098
3099                DO  k = nzb_soil, nzt_soil
3100                   IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
3101                      surf_lsm_h%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
3102                   IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
3103                      surf_lsm_h%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
3104                   IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
3105                      surf_lsm_h%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
3106                   IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
3107                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
3108                   IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
3109                      surf_lsm_h%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
3110                   IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
3111                      surf_lsm_h%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
3112                   IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
3113                      surf_lsm_h%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
3114                   IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
3115                      surf_lsm_h%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
3116                ENDDO
3117
3118             ENDDO
3119!
3120!--          Vertical surfaces
3121             DO  l = 0, 3
3122                DO  m = 1, surf_lsm_v(l)%ns
3123                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
3124                                                   surf_lsm_v(l)%building_covered(m) ) 
3125                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
3126                                                   surf_lsm_v(l)%building_covered(m) ) 
3127
3128                   DO  k = nzb_soil, nzt_soil
3129                      IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
3130                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
3131                      IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
3132                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
3133                      IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
3134                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
3135                      IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
3136                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
3137                      IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
3138                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
3139                      IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
3140                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
3141                      IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
3142                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
3143                      IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
3144                         surf_lsm_v(l)%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
3145                   ENDDO
3146
3147                ENDDO
3148             ENDDO
3149
3150          ENDIF
3151       ENDIF
3152
3153!
3154!--    Level 1, initialization of vegetation parameters. A horizontally
3155!--    homogeneous distribution is assumed here.
3156       IF ( vegetation_type /= 0 )  THEN
3157
3158          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
3159             min_canopy_resistance = vegetation_pars(ind_v_rc_min,vegetation_type)
3160          ENDIF
3161
3162          IF ( leaf_area_index == 9999999.9_wp )  THEN
3163             leaf_area_index = vegetation_pars(ind_v_rc_lai,vegetation_type)         
3164          ENDIF
3165
3166          IF ( vegetation_coverage == 9999999.9_wp )  THEN
3167             vegetation_coverage = vegetation_pars(ind_v_c_veg,vegetation_type)     
3168          ENDIF
3169
3170          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
3171              canopy_resistance_coefficient= vegetation_pars(ind_v_gd,vegetation_type)     
3172          ENDIF
3173
3174          IF ( z0_vegetation == 9999999.9_wp )  THEN
3175             z0_vegetation  = vegetation_pars(ind_v_z0,vegetation_type) 
3176          ENDIF
3177
3178          IF ( z0h_vegetation == 9999999.9_wp )  THEN
3179             z0h_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
3180          ENDIF
3181         
3182          IF ( z0q_vegetation == 9999999.9_wp )  THEN
3183             z0q_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
3184          ENDIF
3185         
3186          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
3187             lambda_surface_stable = vegetation_pars(ind_v_lambda_s,vegetation_type) 
3188          ENDIF
3189
3190          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
3191             lambda_surface_unstable = vegetation_pars(ind_v_lambda_u,vegetation_type)           
3192          ENDIF
3193
3194          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
3195             f_shortwave_incoming = vegetation_pars(ind_v_f_sw_in,vegetation_type)       
3196          ENDIF
3197
3198          IF ( c_surface == 9999999.9_wp )  THEN
3199             c_surface = vegetation_pars(ind_v_c_surf,vegetation_type)       
3200          ENDIF
3201
3202          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
3203             albedo_type = INT(vegetation_pars(ind_v_at,vegetation_type))       
3204          ENDIF
3205   
3206          IF ( emissivity == 9999999.9_wp )  THEN
3207             emissivity = vegetation_pars(ind_v_emis,vegetation_type)     
3208          ENDIF
3209
3210       ENDIF
3211!
3212!--    Map values onto horizontal elemements
3213       DO  m = 1, surf_lsm_h%ns
3214          IF ( surf_lsm_h%vegetation_surface(m) )  THEN
3215             surf_lsm_h%r_canopy_min(m)     = min_canopy_resistance
3216             surf_lsm_h%lai(m)              = leaf_area_index
3217             surf_lsm_h%c_veg(m)            = vegetation_coverage
3218             surf_lsm_h%g_d(m)              = canopy_resistance_coefficient
3219             surf_lsm_h%z0(m)               = z0_vegetation
3220             surf_lsm_h%z0h(m)              = z0h_vegetation
3221             surf_lsm_h%z0q(m)              = z0q_vegetation
3222             surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
3223             surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
3224             surf_lsm_h%f_sw_in(m)          = f_shortwave_incoming
3225             surf_lsm_h%c_surface(m)        = c_surface
3226             surf_lsm_h%albedo_type(ind_veg_wall,m) = albedo_type
3227             surf_lsm_h%emissivity(ind_veg_wall,m)  = emissivity
3228             
3229             surf_lsm_h%vegetation_type(m)      = vegetation_type
3230             surf_lsm_h%vegetation_type_name(m) = vegetation_type_name(vegetation_type)
3231          ELSE
3232             surf_lsm_h%lai(m)   = 0.0_wp
3233             surf_lsm_h%c_veg(m) = 0.0_wp
3234             surf_lsm_h%g_d(m)   = 0.0_wp
3235          ENDIF
3236 
3237       ENDDO
3238!
3239!--    Map values onto vertical elements, even though this does not make
3240!--    much sense.
3241       DO  l = 0, 3
3242          DO  m = 1, surf_lsm_v(l)%ns
3243             IF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
3244                surf_lsm_v(l)%r_canopy_min(m)     = min_canopy_resistance
3245                surf_lsm_v(l)%lai(m)              = leaf_area_index
3246                surf_lsm_v(l)%c_veg(m)            = vegetation_coverage
3247                surf_lsm_v(l)%g_d(m)              = canopy_resistance_coefficient
3248                surf_lsm_v(l)%z0(m)               = z0_vegetation
3249                surf_lsm_v(l)%z0h(m)              = z0h_vegetation
3250                surf_lsm_v(l)%z0q(m)              = z0q_vegetation
3251                surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
3252                surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
3253                surf_lsm_v(l)%f_sw_in(m)          = f_shortwave_incoming
3254                surf_lsm_v(l)%c_surface(m)        = c_surface
3255                surf_lsm_v(l)%albedo_type(ind_veg_wall,m) = albedo_type
3256                surf_lsm_v(l)%emissivity(ind_veg_wall,m)  = emissivity
3257               
3258                surf_lsm_v(l)%vegetation_type(m)      = vegetation_type
3259                surf_lsm_v(l)%vegetation_type_name(m) = vegetation_type_name(vegetation_type)
3260             ELSE
3261                surf_lsm_v(l)%lai(m)   = 0.0_wp
3262                surf_lsm_v(l)%c_veg(m) = 0.0_wp
3263                surf_lsm_v(l)%g_d(m)   = 0.0_wp
3264             ENDIF
3265          ENDDO
3266       ENDDO
3267
3268!
3269!--    Level 2, initialization of vegation parameters via vegetation_type read
3270!--    from file. Vegetation parameters are initialized for each (y,x)-grid point
3271!--    individually using default paramter settings according to the given
3272!--    vegetation type.
3273       IF ( vegetation_type_f%from_file )  THEN
3274!
3275!--       Horizontal surfaces
3276          DO  m = 1, surf_lsm_h%ns
3277             i = surf_lsm_h%i(m)
3278             j = surf_lsm_h%j(m)
3279             
3280             st = vegetation_type_f%var(j,i)
3281             IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
3282                surf_lsm_h%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
3283                surf_lsm_h%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
3284                surf_lsm_h%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
3285                surf_lsm_h%g_d(m)              = vegetation_pars(ind_v_gd,st)
3286                surf_lsm_h%z0(m)               = vegetation_pars(ind_v_z0,st)
3287                surf_lsm_h%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
3288                surf_lsm_h%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
3289                surf_lsm_h%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
3290                surf_lsm_h%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
3291                surf_lsm_h%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
3292                surf_lsm_h%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
3293                surf_lsm_h%albedo_type(ind_veg_wall,m) = INT( vegetation_pars(ind_v_at,st) )
3294                surf_lsm_h%emissivity(ind_veg_wall,m)  = vegetation_pars(ind_v_emis,st)
3295               
3296                surf_lsm_h%vegetation_type(m)      = st
3297                surf_lsm_h%vegetation_type_name(m) = vegetation_type_name(st)
3298             ENDIF
3299          ENDDO
3300!
3301!--       Vertical surfaces
3302          DO  l = 0, 3
3303             DO  m = 1, surf_lsm_v(l)%ns
3304                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3305                                                surf_lsm_v(l)%building_covered(m) ) 
3306                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
3307                                                surf_lsm_v(l)%building_covered(m) ) 
3308             
3309                st = vegetation_type_f%var(j,i)
3310                IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
3311                   surf_lsm_v(l)%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
3312                   surf_lsm_v(l)%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
3313                   surf_lsm_v(l)%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
3314                   surf_lsm_v(l)%g_d(m)              = vegetation_pars(ind_v_gd,st)
3315                   surf_lsm_v(l)%z0(m)               = vegetation_pars(ind_v_z0,st)
3316                   surf_lsm_v(l)%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
3317                   surf_lsm_v(l)%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
3318                   surf_lsm_v(l)%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
3319                   surf_lsm_v(l)%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
3320                   surf_lsm_v(l)%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
3321                   surf_lsm_v(l)%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
3322                   surf_lsm_v(l)%albedo_type(ind_veg_wall,m) = INT( vegetation_pars(ind_v_at,st) )
3323                   surf_lsm_v(l)%emissivity(ind_veg_wall,m)  = vegetation_pars(ind_v_emis,st)
3324                   
3325                   surf_lsm_v(l)%vegetation_type(m)      = st
3326                   surf_lsm_v(l)%vegetation_type_name(m) = vegetation_type_name(st)
3327                ENDIF
3328             ENDDO
3329          ENDDO
3330       ENDIF
3331!
3332!--    Level 3, initialization of vegation parameters at single (x,y)
3333!--    position via vegetation_pars read from file.
3334       IF ( vegetation_pars_f%from_file )  THEN
3335!
3336!--       Horizontal surfaces
3337          DO  m = 1, surf_lsm_h%ns
3338
3339             i = surf_lsm_h%i(m)
3340             j = surf_lsm_h%j(m)
3341!
3342!--          If surface element is not a vegetation surface and any value in
3343!--          vegetation_pars is given, neglect this information and give an
3344!--          informative message that this value will not be used.   
3345             IF ( .NOT. surf_lsm_h%vegetation_surface(m)  .AND.                &
3346                   ANY( vegetation_pars_f%pars_xy(:,j,i) /=                    &
3347                   vegetation_pars_f%fill ) )  THEN
3348                WRITE( message_string, * )                                     &
3349                                 'surface element at grid point (j,i) = (',    &
3350                                 j, i, ') is not a vegation surface, ',        &
3351                                 'so that information given in ',              &
3352                                 'vegetation_pars at this point is neglected.' 
3353                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
3354             ELSE
3355
3356                IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=            &
3357                     vegetation_pars_f%fill )                                  &
3358                   surf_lsm_h%r_canopy_min(m)  =                               &
3359                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
3360                IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=            &
3361                     vegetation_pars_f%fill )                                  &
3362                   surf_lsm_h%lai(m)           =                               &
3363                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
3364                IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=             &
3365                     vegetation_pars_f%fill )                                  &
3366                   surf_lsm_h%c_veg(m)         =                               &
3367                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
3368                IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=                &
3369                     vegetation_pars_f%fill )                                  &
3370                   surf_lsm_h%g_d(m)           =                               &
3371                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
3372                IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=                &
3373                     vegetation_pars_f%fill )                                  &
3374                   surf_lsm_h%z0(m)            =                               &
3375                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
3376                IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=              &
3377                     vegetation_pars_f%fill )  THEN
3378                   surf_lsm_h%z0h(m)           =                               &
3379                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3380                   surf_lsm_h%z0q(m)           =                               &
3381                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3382                ENDIF
3383                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=          &
3384                     vegetation_pars_f%fill )                                  &
3385                   surf_lsm_h%lambda_surface_s(m) =                            &
3386                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
3387                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=          &
3388                     vegetation_pars_f%fill )                                  &
3389                   surf_lsm_h%lambda_surface_u(m) =                            &
3390                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
3391                IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=           &
3392                     vegetation_pars_f%fill )                                  &
3393                   surf_lsm_h%f_sw_in(m)          =                            &
3394                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
3395                IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=            &
3396                     vegetation_pars_f%fill )                                  &
3397                   surf_lsm_h%c_surface(m)        =                            &
3398                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
3399                IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=                &
3400                     vegetation_pars_f%fill )                                  &
3401                   surf_lsm_h%albedo_type(ind_veg_wall,m) =                    &
3402                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
3403                IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=              &
3404                     vegetation_pars_f%fill )                                  &
3405                   surf_lsm_h%emissivity(ind_veg_wall,m)  =                    &
3406                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
3407             ENDIF
3408          ENDDO
3409!
3410!--       Vertical surfaces
3411          DO  l = 0, 3
3412             DO  m = 1, surf_lsm_v(l)%ns
3413                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3414                                                surf_lsm_v(l)%building_covered(m) ) 
3415                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3416                                                surf_lsm_v(l)%building_covered(m) ) 
3417!
3418!--             If surface element is not a vegetation surface and any value in
3419!--             vegetation_pars is given, neglect this information and give an
3420!--             informative message that this value will not be used.   
3421                IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m)  .AND.          &
3422                      ANY( vegetation_pars_f%pars_xy(:,j,i) /=                 &
3423                      vegetation_pars_f%fill ) )  THEN
3424                   WRITE( message_string, * )                                  &
3425                                 'surface element at grid point (j,i) = (',    &
3426                                 j, i, ') is not a vegation surface, ',        &
3427                                 'so that information given in ',              &
3428                                 'vegetation_pars at this point is neglected.' 
3429                   CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
3430                ELSE
3431
3432                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=         &
3433                        vegetation_pars_f%fill )                               &
3434                      surf_lsm_v(l)%r_canopy_min(m)  =                         &
3435                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
3436                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=         &
3437                        vegetation_pars_f%fill )                               &
3438                      surf_lsm_v(l)%lai(m)           =                         &
3439                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
3440                   IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=          &
3441                        vegetation_pars_f%fill )                               &
3442                      surf_lsm_v(l)%c_veg(m)         =                         &
3443                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
3444                   IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=             &
3445                        vegetation_pars_f%fill )                               &
3446                     surf_lsm_v(l)%g_d(m)            =                         &
3447                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
3448                   IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=             &
3449                        vegetation_pars_f%fill )                               &
3450                      surf_lsm_v(l)%z0(m)            =                         &
3451                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
3452                   IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=           &
3453                        vegetation_pars_f%fill )  THEN
3454                      surf_lsm_v(l)%z0h(m)           =                         &
3455                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3456                      surf_lsm_v(l)%z0q(m)           =                         &
3457                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
3458                   ENDIF
3459                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=       &
3460                        vegetation_pars_f%fill )                               &
3461                      surf_lsm_v(l)%lambda_surface_s(m)  =                     &
3462                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
3463                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=       &
3464                        vegetation_pars_f%fill )                               &
3465                      surf_lsm_v(l)%lambda_surface_u(m)  =                     &
3466                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
3467                   IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=        &
3468                        vegetation_pars_f%fill )                               &
3469                      surf_lsm_v(l)%f_sw_in(m)           =                     &
3470                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
3471                   IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=         &
3472                        vegetation_pars_f%fill )                               &
3473                      surf_lsm_v(l)%c_surface(m)         =                     &
3474                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
3475                   IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=             &
3476                        vegetation_pars_f%fill )                               &
3477                      surf_lsm_v(l)%albedo_type(ind_veg_wall,m) =              &
3478                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
3479                   IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=           &
3480                        vegetation_pars_f%fill )                               &
3481                      surf_lsm_v(l)%emissivity(ind_veg_wall,m)  =              &
3482                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
3483                ENDIF
3484
3485             ENDDO
3486          ENDDO
3487       ENDIF 
3488
3489!
3490!--    Level 1, initialization of water parameters. A horizontally
3491!--    homogeneous distribution is assumed here.
3492       IF ( water_type /= 0 )  THEN
3493
3494          IF ( water_temperature == 9999999.9_wp )  THEN
3495             water_temperature = water_pars(ind_w_temp,water_type)       
3496          ENDIF
3497
3498          IF ( z0_water == 9999999.9_wp )  THEN
3499             z0_water = water_pars(ind_w_z0,water_type)       
3500          ENDIF       
3501
3502          IF ( z0h_water == 9999999.9_wp )  THEN
3503             z0h_water = water_pars(ind_w_z0h,water_type)       
3504          ENDIF 
3505         
3506          IF ( z0q_water == 9999999.9_wp )  THEN
3507             z0q_water = water_pars(ind_w_z0h,water_type)       
3508          ENDIF
3509
3510          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
3511             albedo_type = INT(water_pars(ind_w_at,water_type))       
3512          ENDIF
3513   
3514          IF ( emissivity == 9999999.9_wp )  THEN
3515             emissivity = water_pars(ind_w_emis,water_type)       
3516          ENDIF
3517
3518       ENDIF 
3519!
3520!--    Map values onto horizontal elemements
3521       DO  m = 1, surf_lsm_h%ns
3522          IF ( surf_lsm_h%water_surface(m) )  THEN
3523             IF ( TRIM( initializing_actions ) /= 'read_restart_data' )        &
3524                t_soil_h%var_2d(:,m)        = water_temperature
3525             surf_lsm_h%z0(m)               = z0_water
3526             surf_lsm_h%z0h(m)              = z0h_water
3527             surf_lsm_h%z0q(m)              = z0q_water
3528             surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
3529             surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp               
3530             surf_lsm_h%c_surface(m)        = 0.0_wp
3531             surf_lsm_h%albedo_type(ind_wat_win,m) = albedo_type
3532             surf_lsm_h%emissivity(ind_wat_win,m)  = emissivity
3533             
3534             surf_lsm_h%water_type(m)      = water_type
3535             surf_lsm_h%water_type_name(m) = water_type_name(water_type)
3536          ENDIF
3537       ENDDO
3538!
3539!--    Map values onto vertical elements, even though this does not make
3540!--    much sense.
3541       DO  l = 0, 3
3542          DO  m = 1, surf_lsm_v(l)%ns
3543             IF ( surf_lsm_v(l)%water_surface(m) )  THEN
3544                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
3545                   t_soil_v(l)%var_2d(:,m)           = water_temperature
3546                surf_lsm_v(l)%z0(m)               = z0_water
3547                surf_lsm_v(l)%z0h(m)              = z0h_water
3548                surf_lsm_v(l)%z0q(m)              = z0q_water
3549                surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp
3550                surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp               
3551                surf_lsm_v(l)%c_surface(m)        = 0.0_wp
3552                surf_lsm_v(l)%albedo_type(ind_wat_win,m) = albedo_type
3553                surf_lsm_v(l)%emissivity(ind_wat_win,m)  = emissivity
3554               
3555                surf_lsm_v(l)%water_type(m)      = water_type
3556                surf_lsm_v(l)%water_type_name(m) = water_type_name(water_type)
3557             ENDIF
3558          ENDDO
3559       ENDDO
3560!
3561!
3562!--    Level 2, initialization of water parameters via water_type read
3563!--    from file. Water surfaces are initialized for each (y,x)-grid point
3564!--    individually using default paramter settings according to the given
3565!--    water type.
3566!--    Note, parameter 3/4 of water_pars are albedo and emissivity,
3567!--    whereas paramter 3/4 of water_pars_f are heat conductivities!
3568       IF ( water_type_f%from_file )  THEN
3569!
3570!--       Horizontal surfaces
3571          DO  m = 1, surf_lsm_h%ns
3572             i = surf_lsm_h%i(m)
3573             j = surf_lsm_h%j(m)
3574             
3575             st = water_type_f%var(j,i)
3576             IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
3577                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
3578                   t_soil_h%var_2d(:,m) = water_pars(ind_w_temp,st)
3579                surf_lsm_h%z0(m)     = water_pars(ind_w_z0,st)
3580                surf_lsm_h%z0h(m)    = water_pars(ind_w_z0h,st)
3581                surf_lsm_h%z0q(m)    = water_pars(ind_w_z0h,st)
3582                surf_lsm_h%lambda_surface_s(m) = water_pars(ind_w_lambda_s,st)
3583                surf_lsm_h%lambda_surface_u(m) = water_pars(ind_w_lambda_u,st)             
3584                surf_lsm_h%c_surface(m)        = 0.0_wp
3585                surf_lsm_h%albedo_type(ind_wat_win,m) = INT( water_pars(ind_w_at,st) )
3586                surf_lsm_h%emissivity(ind_wat_win,m)  = water_pars(ind_w_emis,st)
3587               
3588                surf_lsm_h%water_type(m)      = st
3589                surf_lsm_h%water_type_name(m) = water_type_name(st)
3590             ENDIF
3591          ENDDO
3592!
3593!--       Vertical surfaces
3594          DO  l = 0, 3
3595             DO  m = 1, surf_lsm_v(l)%ns
3596                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3597                                                surf_lsm_v(l)%building_covered(m) ) 
3598                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3599                                                surf_lsm_v(l)%building_covered(m) ) 
3600             
3601                st = water_type_f%var(j,i)
3602                IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
3603                   IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  &
3604                      t_soil_v(l)%var_2d(:,m) = water_pars(ind_w_temp,st)
3605                   surf_lsm_v(l)%z0(m)     = water_pars(ind_w_z0,st)
3606                   surf_lsm_v(l)%z0h(m)    = water_pars(ind_w_z0h,st)
3607                   surf_lsm_v(l)%z0q(m)    = water_pars(ind_w_z0h,st)
3608                   surf_lsm_v(l)%lambda_surface_s(m) =                         &
3609                                                   water_pars(ind_w_lambda_s,st)
3610                   surf_lsm_v(l)%lambda_surface_u(m) =                         &
3611                                                   water_pars(ind_w_lambda_u,st)           
3612                   surf_lsm_v(l)%c_surface(m)     = 0.0_wp
3613                   surf_lsm_v(l)%albedo_type(ind_wat_win,m) =                  &
3614                                                  INT( water_pars(ind_w_at,st) )
3615                   surf_lsm_v(l)%emissivity(ind_wat_win,m)  =                  &
3616                                                  water_pars(ind_w_emis,st)
3617                                                 
3618                   surf_lsm_v(l)%water_type(m)      = st
3619                   surf_lsm_v(l)%water_type_name(m) = water_type_name(st)
3620                ENDIF
3621             ENDDO
3622          ENDDO
3623       ENDIF     
3624
3625!
3626!--    Level 3, initialization of water parameters at single (x,y)
3627!--    position via water_pars read from file.
3628       IF ( water_pars_f%from_file )  THEN
3629!
3630!--       Horizontal surfaces
3631          DO  m = 1, surf_lsm_h%ns
3632             i = surf_lsm_h%i(m)
3633             j = surf_lsm_h%j(m)
3634!
3635!--          If surface element is not a water surface and any value in
3636!--          water_pars is given, neglect this information and give an
3637!--          informative message that this value will not be used.   
3638             IF ( .NOT. surf_lsm_h%water_surface(m)  .AND.                     &
3639                   ANY( water_pars_f%pars_xy(:,j,i) /= water_pars_f%fill ) )  THEN
3640                WRITE( message_string, * )                                     &
3641                              'surface element at grid point (j,i) = (',       &
3642                              j, i, ') is not a water surface, ',              &
3643                              'so that information given in ',                 &
3644                              'water_pars at this point is neglected.' 
3645                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
3646             ELSE
3647                IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                   &
3648                     water_pars_f%fill  .AND.                                  &
3649                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
3650                      t_soil_h%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
3651
3652                IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /= water_pars_f%fill ) &
3653                   surf_lsm_h%z0(m)     = water_pars_f%pars_xy(ind_w_z0,j,i)
3654
3655                IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /= water_pars_f%fill )&
3656                THEN
3657                   surf_lsm_h%z0h(m)    = water_pars_f%pars_xy(ind_w_z0h,j,i)
3658                   surf_lsm_h%z0q(m)    = water_pars_f%pars_xy(ind_w_z0h,j,i)
3659                ENDIF
3660                IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /=               &
3661                     water_pars_f%fill )                                       &
3662                   surf_lsm_h%lambda_surface_s(m) =                            &
3663                                        water_pars_f%pars_xy(ind_w_lambda_s,j,i)
3664
3665                IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /=               &
3666                      water_pars_f%fill )                                      &
3667                   surf_lsm_h%lambda_surface_u(m) =                            &
3668                                        water_pars_f%pars_xy(ind_w_lambda_u,j,i)     
3669       
3670                IF ( water_pars_f%pars_xy(ind_w_at,j,i) /=                     &
3671                     water_pars_f%fill )                                       &
3672                   surf_lsm_h%albedo_type(ind_wat_win,m) =                     &
3673                                       INT( water_pars_f%pars_xy(ind_w_at,j,i) )
3674
3675                IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /=                   &
3676                     water_pars_f%fill )                                       &
3677                   surf_lsm_h%emissivity(ind_wat_win,m) =                      &
3678                                          water_pars_f%pars_xy(ind_w_emis,j,i) 
3679             ENDIF
3680          ENDDO
3681!
3682!--       Vertical surfaces
3683          DO  l = 0, 3
3684             DO  m = 1, surf_lsm_v(l)%ns
3685                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
3686                                                surf_lsm_v(l)%building_covered(m) ) 
3687                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
3688                                                surf_lsm_v(l)%building_covered(m) ) 
3689!
3690!--             If surface element is not a water surface and any value in
3691!--             water_pars is given, neglect this information and give an
3692!--             informative message that this value will not be used.   
3693                IF ( .NOT. surf_lsm_v(l)%water_surface(m)  .AND.               &
3694                      ANY( water_pars_f%pars_xy(:,j,i) /=                      &
3695                      water_pars_f%fill ) )  THEN
3696                   WRITE( message_string, * )                                  &
3697                              'surface element at grid point (j,i) = (',       &
3698                              j, i, ') is not a water surface, ',              &
3699                              'so that information given in ',                 &
3700                              'water_pars at this point is neglected.' 
3701                   CALL message( 'land_surface_model_mod', 'PA0999',           &
3702                                  0, 0, 0, 6, 0 )
3703                ELSE
3704
3705                   IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                &
3706                     water_pars_f%fill  .AND.                                  &
3707                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
3708                      t_soil_v(l)%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
3709
3710                   IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /=                  &
3711                        water_pars_f%fill )                                    &
3712                      surf_lsm_v(l)%z0(m)   = water_pars_f%pars_xy(ind_w_z0,j,i)
3713
3714                   IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /=                 &
3715                       water_pars_f%fill )  THEN
3716                      surf_lsm_v(l)%z0h(m)  = water_pars_f%pars_xy(ind_w_z0h,j,i)
3717                      surf_lsm_v(l)%z0q(m)  = water_pars_f%pars_xy(ind_w_z0h,j,i)
3718                   ENDIF
3719
3720                   IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /=            &
3721                        water_pars_f%fill )                                    &
3722                      surf_lsm_v(l)%lambda_surface_s(m) =                      &
3723                                      water_pars_f%pars_xy(ind_w_lambda_s,j,i)
3724
3725                   IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /=            &
3726                        water_pars_f%fill )                                    &
3727                      surf_lsm_v(l)%lambda_surface_u(m) =                      &
3728                                      water_pars_f%pars_xy(ind_w_lambda_u,j,i)   
3729 
3730                   IF ( water_pars_f%pars_xy(ind_w_at,j,i) /=                  &
3731                        water_pars_f%fill )                                    &
3732                      surf_lsm_v(l)%albedo_type(ind_wat_win,m) =               &
3733                                      INT( water_pars_f%pars_xy(ind_w_at,j,i) )
3734
3735                   IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /=                &
3736                        water_pars_f%fill )                                    &
3737                      surf_lsm_v(l)%emissivity(ind_wat_win,m)  =               &
3738                                      water_pars_f%pars_xy(ind_w_emis,j,i) 
3739                ENDIF
3740             ENDDO
3741          ENDDO
3742
3743       ENDIF
3744!
3745!--    Initialize pavement-type surfaces, level 1
3746       IF ( pavement_type /= 0 )  THEN 
3747
3748!
3749!--       When a pavement_type is used, overwrite a possible setting of
3750!--       the pavement depth as it is already defined by the pavement type
3751          pavement_depth_level = 0
3752
3753          IF ( z0_pavement == 9999999.9_wp )  THEN
3754             z0_pavement  = pavement_pars(ind_p_z0,pavement_type) 
3755          ENDIF
3756
3757          IF ( z0h_pavement == 9999999.9_wp )  THEN
3758             z0h_pavement = pavement_pars(ind_p_z0h,pavement_type)
3759          ENDIF
3760         
3761          IF ( z0q_pavement == 9999999.9_wp )  THEN
3762             z0q_pavement = pavement_pars(ind_p_z0h,pavement_type)
3763          ENDIF
3764
3765          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
3766             pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type)
3767          ENDIF
3768
3769          IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
3770             pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)
3771          ENDIF   
3772   
3773          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
3774             albedo_type = INT(pavement_pars(ind_p_at,pavement_type))       
3775          ENDIF
3776   
3777          IF ( emissivity == 9999999.9_wp )  THEN
3778             emissivity = pavement_pars(ind_p_emis,pavement_type)       
3779          ENDIF
3780
3781!
3782!--       If the depth level of the pavement is not set, determine it from
3783!--       lookup table.
3784          IF ( pavement_depth_level == 0 )  THEN
3785             DO  k = nzb_soil, nzt_soil 
3786                IF ( pavement_subsurface_pars_1(k,pavement_type) == 9999999.9_wp &
3787                .OR. pavement_subsurface_pars_2(k,pavement_type) == 9999999.9_wp)&
3788                THEN
3789                   nzt_pavement = k-1
3790                   EXIT
3791                ENDIF
3792             ENDDO
3793          ELSE
3794             nzt_pavement = pavement_depth_level
3795          ENDIF
3796
3797       ENDIF
3798!
3799!--    Level 1 initialization of pavement type surfaces. Horizontally
3800!--    homogeneous characteristics are assumed
3801       surf_lsm_h%nzt_pavement = pavement_depth_level
3802       DO  m = 1, surf_lsm_h%ns
3803          IF ( surf_lsm_h%pavement_surface(m) )  THEN
3804             surf_lsm_h%nzt_pavement(m)        = nzt_pavement
3805             surf_lsm_h%z0(m)                  = z0_pavement
3806             surf_lsm_h%z0h(m)                 = z0h_pavement
3807             surf_lsm_h%z0q(m)                 = z0q_pavement
3808             surf_lsm_h%lambda_surface_s(m)    = pavement_heat_conduct         &
3809                                                  * ddz_soil(nzb_soil)         &
3810                                                  * 2.0_wp   
3811             surf_lsm_h%lambda_surface_u(m)    = pavement_heat_conduct         &
3812                                                  * ddz_soil(nzb_soil)         &
3813                                                  * 2.0_wp           
3814             surf_lsm_h%c_surface(m)           = pavement_heat_capacity        &
3815                                                        * dz_soil(nzb_soil)    &
3816                                                        * 0.25_wp                                   
3817
3818             surf_lsm_h%albedo_type(ind_pav_green,m) = albedo_type
3819             surf_lsm_h%emissivity(ind_pav_green,m)  = emissivity     
3820             
3821             surf_lsm_h%pavement_type(m)      = pavement_type
3822             surf_lsm_h%pavement_type_name(m) = pavement_type_name(pavement_type)
3823     
3824             IF ( pavement_type /= 0 )  THEN
3825                DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
3826                   surf_lsm_h%lambda_h_def(k,m)    =                           &
3827                                     pavement_subsurface_pars_1(k,pavement_type)                       
3828                   surf_lsm_h%rho_c_total_def(k,m) =                           &
3829                                     pavement_subsurface_pars_2(k,pavement_type) 
3830                ENDDO
3831             ELSE
3832                surf_lsm_h%lambda_h_def(:,m)     = pavement_heat_conduct
3833                surf_lsm_h%rho_c_total_def(:,m)  = pavement_heat_capacity
3834             ENDIF       
3835          ENDIF
3836       ENDDO                               
3837
3838       DO  l = 0, 3
3839          surf_lsm_v(l)%nzt_pavement = pavement_depth_level
3840          DO  m = 1, surf_lsm_v(l)%ns
3841             IF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
3842                surf_lsm_v(l)%nzt_pavement(m)        = nzt_pavement
3843                surf_lsm_v(l)%z0(m)                  = z0_pavement
3844                surf_lsm_v(l)%z0h(m)                 = z0h_pavement
3845                surf_lsm_v(l)%z0q(m)                 = z0q_pavement
3846                surf_lsm_v(l)%lambda_surface_s(m)    = pavement_heat_conduct   &
3847                                                  * ddz_soil(nzb_soil)         &
3848                                                  * 2.0_wp   
3849                surf_lsm_v(l)%lambda_surface_u(m)    = pavement_heat_conduct   &
3850                                                  * ddz_soil(nzb_soil)         &
3851                                                  * 2.0_wp           
3852                surf_lsm_v(l)%c_surface(m)           = pavement_heat_capacity  &
3853                                                        * dz_soil(nzb_soil)    &
3854                                                        * 0.25_wp                                     
3855
3856                surf_lsm_v(l)%albedo_type(ind_pav_green,m) = albedo_type
3857                surf_lsm_v(l)%emissivity(ind_pav_green,m)  = emissivity
3858               
3859                surf_lsm_v(l)%pavement_type(m)      = pavement_type
3860                surf_lsm_v(l)%pavement_type_name(m) = pavement_type_name