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

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

last commit documented

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