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

Last change on this file since 2608 was 2608, checked in by schwenkel, 4 years ago

Inital revision of diagnostic_quantities_mod allows unified calculation of magnus equation and saturion mixing ratio

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