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

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

minor bugfixes in land surface model

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