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

Last change on this file since 1976 was 1976, checked in by maronga, 5 years ago

further modularization of land surface model (2D/3D output and restart data). Bugfix for restart runs without land surface model

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