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

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

last commit documented

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