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

Last change on this file since 1980 was 1979, checked in by maronga, 8 years ago

last commit documented

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