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

Last change on this file since 2575 was 2575, checked in by maronga, 4 years ago

bugfix in radiation model and improvements in land surface scheme

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