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

Last change on this file since 2516 was 2516, checked in by suehring, 4 years ago

document changes

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