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

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

bugfix in land surface model

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