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

Last change on this file since 2172 was 2150, checked in by scharf, 7 years ago

last commit documented

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