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

Last change on this file since 2024 was 2001, checked in by knoop, 8 years ago

last commit documented

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