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

Last change on this file since 2532 was 2532, checked in by scharf, 4 years ago

bugfixes in lsm_data_output_3d

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