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

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

bugfix for lsm data output

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