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

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

Bugfix for r2547

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