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

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

bugfix in land surface model

  • Property svn:keywords set to Id
File size: 141.6 KB
Line 
1!> @file land_surface_model_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21! Bugfix: initial values of pave_surface and water_surface were not set.
22!
23! Former revisions:
24! -----------------
25! $Id: land_surface_model_mod.f90 1978 2016-07-29 12:08:31Z 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       water_surface = .FALSE.
2001       pave_surface  = .FALSE.
2002
2003#if ! defined( __nopointer )
2004!
2005!--    Initial assignment of the pointers
2006       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
2007       t_surface => t_surface_1; t_surface_p => t_surface_2
2008       m_soil    => m_soil_1;    m_soil_p    => m_soil_2
2009       m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
2010#endif
2011
2012
2013    END SUBROUTINE lsm_init_arrays
2014
2015
2016!------------------------------------------------------------------------------!
2017! Description:
2018! ------------
2019!> Parin for &lsmpar for land surface model
2020!------------------------------------------------------------------------------!
2021    SUBROUTINE lsm_parin
2022
2023
2024       IMPLICIT NONE
2025
2026       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
2027       
2028       NAMELIST /lsm_par/         alpha_vangenuchten, c_surface,               &
2029                                  canopy_resistance_coefficient,               &
2030                                  conserve_water_content,                      &
2031                                  f_shortwave_incoming, field_capacity,        & 
2032                                  hydraulic_conductivity,                      &
2033                                  lambda_surface_stable,                       &
2034                                  lambda_surface_unstable, leaf_area_index,    &
2035                                  l_vangenuchten, min_canopy_resistance,       &
2036                                  min_soil_resistance, n_vangenuchten,         &
2037                                  pave_depth, pave_heat_capacity,              &
2038                                  pave_heat_conductivity,                      &
2039                                  residual_moisture, root_fraction,            &
2040                                  saturation_moisture, skip_time_do_lsm,       &
2041                                  soil_moisture, soil_temperature, soil_type,  &
2042                                  vegetation_coverage, veg_type, wilting_point,& 
2043                                  zs, z0_eb, z0h_eb, z0q_eb
2044       
2045       line = ' '
2046       
2047!
2048!--    Try to find land surface model package
2049       REWIND ( 11 )
2050       line = ' '
2051       DO   WHILE ( INDEX( line, '&lsm_par' ) == 0 )
2052          READ ( 11, '(A)', END=10 )  line
2053       ENDDO
2054       BACKSPACE ( 11 )
2055
2056!
2057!--    Read user-defined namelist
2058       READ ( 11, lsm_par )
2059
2060!
2061!--    Set flag that indicates that the land surface model is switched on
2062       land_surface = .TRUE.
2063
2064 10    CONTINUE
2065       
2066
2067    END SUBROUTINE lsm_parin
2068
2069
2070!------------------------------------------------------------------------------!
2071! Description:
2072! ------------
2073!> Soil model as part of the land surface model. The model predicts soil
2074!> temperature and water content.
2075!------------------------------------------------------------------------------!
2076    SUBROUTINE lsm_soil_model
2077
2078
2079       IMPLICIT NONE
2080
2081       INTEGER(iwp) ::  i   !< running index
2082       INTEGER(iwp) ::  j   !< running index
2083       INTEGER(iwp) ::  k   !< running index
2084
2085       REAL(wp)     :: h_vg !< Van Genuchten coef. h
2086
2087       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
2088                                                 lambda_temp, & !< temp. lambda
2089                                                 tend           !< tendency
2090
2091       DO  i = nxlg, nxrg
2092          DO  j = nysg, nyng
2093
2094             IF ( pave_surface(j,i) )  THEN
2095                rho_c_total(nzb_soil,j,i) = pave_heat_capacity
2096                lambda_temp(nzb_soil)     = pave_heat_conductivity
2097             ENDIF
2098
2099             IF (  .NOT.  water_surface(j,i) )  THEN
2100                DO  k = nzb_soil, nzt_soil
2101
2102
2103                   IF ( pave_surface(j,i)  .AND.  zs(k) <= pave_depth )  THEN
2104                   
2105                      rho_c_total(k,j,i) = pave_heat_capacity
2106                      lambda_temp(k)     = pave_heat_conductivity   
2107
2108                   ELSE           
2109!
2110!--                   Calculate volumetric heat capacity of the soil, taking
2111!--                   into account water content
2112                      rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) &
2113                                           + rho_c_water * m_soil(k,j,i))
2114
2115!
2116!--                   Calculate soil heat conductivity at the center of the soil
2117!--                   layers
2118                      lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *    &
2119                                     lambda_h_water ** m_soil(k,j,i)
2120
2121                      ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i)             &
2122                                                     / m_sat(j,i)))
2123
2124                      lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +    &
2125                                       lambda_h_dry
2126                   ENDIF
2127
2128                ENDDO
2129
2130!
2131!--             Calculate soil heat conductivity (lambda_h) at the _stag level
2132!--             using linear interpolation. For pavement surface, the
2133!--             true pavement depth is considered
2134                DO  k = nzb_soil, nzt_soil-1
2135                   IF ( pave_surface(j,i)  .AND.  zs(k)   < pave_depth         &
2136                                           .AND.  zs(k+1) > pave_depth )  THEN
2137                      lambda_h(k,j,i) = ( pave_depth - zs(k) ) / dz_soil(k+1)  &
2138                                        * lambda_temp(k)                       &
2139                                        + ( 1.0_wp - ( pave_depth - zs(k) )    &
2140                                        / dz_soil(k+1) ) * lambda_temp(k+1)
2141                   ELSE
2142                      lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
2143                                        * 0.5_wp
2144                   ENDIF
2145                ENDDO
2146                lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
2147
2148
2149
2150
2151!
2152!--             Prognostic equation for soil temperature t_soil
2153                tend(:) = 0.0_wp
2154
2155                tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *          &
2156                          ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)  &
2157                            - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)    &
2158                            + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
2159
2160                DO  k = nzb_soil+1, nzt_soil
2161                   tend(k) = (1.0_wp/rho_c_total(k,j,i))                       &
2162                             * (   lambda_h(k,j,i)                             &
2163                                 * ( t_soil(k+1,j,i) - t_soil(k,j,i) )         &
2164                                 * ddz_soil(k+1)                               &
2165                                 - lambda_h(k-1,j,i)                           &
2166                                 * ( t_soil(k,j,i) - t_soil(k-1,j,i) )         &
2167                                 * ddz_soil(k)                                 &
2168                               ) * ddz_soil_stag(k)
2169
2170                ENDDO
2171
2172                t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)&
2173                                                  + dt_3d * ( tsc(2)           &
2174                                                  * tend(nzb_soil:nzt_soil)    & 
2175                                                  + tsc(3)                     &
2176                                                  * tt_soil_m(:,j,i) )   
2177
2178!
2179!--             Calculate t_soil tendencies for the next Runge-Kutta step
2180                IF ( timestep_scheme(1:5) == 'runge' )  THEN
2181                   IF ( intermediate_timestep_count == 1 )  THEN
2182                      DO  k = nzb_soil, nzt_soil
2183                         tt_soil_m(k,j,i) = tend(k)
2184                      ENDDO
2185                   ELSEIF ( intermediate_timestep_count <                      &
2186                            intermediate_timestep_count_max )  THEN
2187                      DO  k = nzb_soil, nzt_soil
2188                         tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
2189                                         * tt_soil_m(k,j,i)
2190                      ENDDO
2191                   ENDIF
2192                ENDIF
2193
2194
2195                DO  k = nzb_soil, nzt_soil
2196
2197!
2198!--                Calculate soil diffusivity at the center of the soil layers
2199                   lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat       &
2200                                     / m_sat(j,i) ) * ( MAX( m_soil(k,j,i),    &
2201                                     m_wilt(j,i) ) / m_sat(j,i) )**(           &
2202                                     b_ch + 2.0_wp )
2203
2204!
2205!--                Parametrization of Van Genuchten
2206                   IF ( soil_type /= 7 )  THEN
2207!
2208!--                   Calculate the hydraulic conductivity after Van Genuchten
2209!--                   (1980)
2210                      h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -    &
2211                                 MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**(      &
2212                                 n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp  &
2213                             )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i)
2214
2215
2216                      gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +         &
2217                                      ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**(  &
2218                                      1.0_wp - 1.0_wp / n_vg(j,i) ) - (        &
2219                                      alpha_vg(j,i) * h_vg )**( n_vg(j,i)      &
2220                                      - 1.0_wp) )**2 )                         &
2221                                      / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg    &
2222                                      )**n_vg(j,i) )**( ( 1.0_wp  - 1.0_wp     &
2223                                      / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) )
2224
2225!
2226!--                Parametrization of Clapp & Hornberger
2227                   ELSE
2228                      gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i)       &
2229                                      / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
2230                   ENDIF
2231
2232                ENDDO
2233
2234!
2235!--             Prognostic equation for soil moisture content. Only performed,
2236!--             when humidity is enabled in the atmosphere and the surface type
2237!--             is not pavement (implies dry soil below).
2238                IF ( humidity  .AND.  .NOT.  pave_surface(j,i) )  THEN
2239!
2240!--                Calculate soil diffusivity (lambda_w) at the _stag level
2241!--                using linear interpolation. To do: replace this with
2242!--                ECMWF-IFS Eq. 8.81
2243                   DO  k = nzb_soil, nzt_soil-1
2244                     
2245                      lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
2246                                        * 0.5_wp
2247                      gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )    &
2248                                        * 0.5_wp
2249
2250                   ENDDO
2251
2252!
2253!
2254!--                In case of a closed bottom (= water content is conserved),
2255!--                set hydraulic conductivity to zero to that no water will be
2256!--                lost in the bottom layer.
2257                   IF ( conserve_water_content )  THEN
2258                      gamma_w(nzt_soil,j,i) = 0.0_wp
2259                   ELSE
2260                      gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
2261                   ENDIF     
2262
2263!--                The root extraction (= root_extr * qsws_veg_eb / (rho_l     
2264!--                * l_v)) ensures the mass conservation for water. The         
2265!--                transpiration of plants equals the cumulative withdrawals by
2266!--                the roots in the soil. The scheme takes into account the
2267!--                availability of water in the soil layers as well as the root
2268!--                fraction in the respective layer. Layer with moisture below
2269!--                wilting point will not contribute, which reflects the
2270!--                preference of plants to take water from moister layers.
2271
2272!
2273!--                Calculate the root extraction (ECMWF 7.69, the sum of
2274!--                root_extr = 1). The energy balance solver guarantees a
2275!--                positive transpiration, so that there is no need for an
2276!--                additional check.
2277                   m_total = 0.0_wp
2278                   DO  k = nzb_soil, nzt_soil
2279                       IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
2280                          m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
2281                       ENDIF
2282                   ENDDO 
2283
2284                   IF ( m_total > 0.0_wp )  THEN
2285                      DO  k = nzb_soil, nzt_soil
2286                         IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
2287                            root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i)      &
2288                                                            / m_total
2289                         ELSE
2290                            root_extr(k) = 0.0_wp
2291                         ENDIF
2292                      ENDDO
2293                   ENDIF
2294
2295!
2296!--                Prognostic equation for soil water content m_soil.
2297                   tend(:) = 0.0_wp
2298
2299                   tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (               &
2300                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
2301                            * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
2302                               root_extr(nzb_soil) * qsws_veg_eb(j,i)          &
2303                               + qsws_soil_eb(j,i) ) * drho_l_lv )             &
2304                               * ddz_soil_stag(nzb_soil)
2305
2306                   DO  k = nzb_soil+1, nzt_soil-1
2307                      tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)          &
2308                                - m_soil(k,j,i) ) * ddz_soil(k+1)              &
2309                                - gamma_w(k,j,i)                               &
2310                                - lambda_w(k-1,j,i) * (m_soil(k,j,i) -         &
2311                                m_soil(k-1,j,i)) * ddz_soil(k)                 &
2312                                + gamma_w(k-1,j,i) - (root_extr(k)             &
2313                                * qsws_veg_eb(j,i) * drho_l_lv)                &
2314                                ) * ddz_soil_stag(k)
2315
2316                   ENDDO
2317                   tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                  &
2318                                           - lambda_w(nzt_soil-1,j,i)          &
2319                                           * (m_soil(nzt_soil,j,i)             &
2320                                           - m_soil(nzt_soil-1,j,i))           &
2321                                           * ddz_soil(nzt_soil)                &
2322                                           + gamma_w(nzt_soil-1,j,i) - (       &
2323                                             root_extr(nzt_soil)               &
2324                                           * qsws_veg_eb(j,i) * drho_l_lv  )   &
2325                                     ) * ddz_soil_stag(nzt_soil)             
2326
2327                   m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
2328                                                   + dt_3d * ( tsc(2) * tend(:)   &
2329                                                   + tsc(3) * tm_soil_m(:,j,i) )   
2330   
2331!
2332!--                Account for dry soils (find a better solution here!)
2333                   DO  k = nzb_soil, nzt_soil
2334                      IF ( m_soil_p(k,j,i) < 0.0_wp )  m_soil_p(k,j,i) = 0.0_wp
2335                   ENDDO
2336
2337!
2338!--                Calculate m_soil tendencies for the next Runge-Kutta step
2339                   IF ( timestep_scheme(1:5) == 'runge' )  THEN
2340                      IF ( intermediate_timestep_count == 1 )  THEN
2341                         DO  k = nzb_soil, nzt_soil
2342                            tm_soil_m(k,j,i) = tend(k)
2343                         ENDDO
2344                      ELSEIF ( intermediate_timestep_count <                   &
2345                               intermediate_timestep_count_max )  THEN
2346                         DO  k = nzb_soil, nzt_soil
2347                            tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp&
2348                                     * tm_soil_m(k,j,i)
2349                         ENDDO
2350                      ENDIF
2351                   ENDIF
2352
2353                ENDIF
2354
2355             ENDIF
2356
2357          ENDDO
2358       ENDDO
2359
2360    END SUBROUTINE lsm_soil_model
2361
2362 
2363!------------------------------------------------------------------------------!
2364! Description:
2365! ------------
2366!> Swapping of timelevels
2367!------------------------------------------------------------------------------!
2368    SUBROUTINE lsm_swap_timelevel ( mod_count )
2369
2370       IMPLICIT NONE
2371
2372       INTEGER, INTENT(IN) :: mod_count
2373
2374#if defined( __nopointer )
2375
2376       t_surface    = t_surface_p
2377       t_soil       = t_soil_p
2378       IF ( humidity )  THEN
2379          m_soil    = m_soil_p
2380          m_liq_eb  = m_liq_eb_p
2381       ENDIF
2382
2383#else
2384   
2385       SELECT CASE ( mod_count )
2386
2387          CASE ( 0 )
2388
2389             t_surface  => t_surface_1; t_surface_p  => t_surface_2
2390             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
2391             IF ( humidity )  THEN
2392                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
2393                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
2394             ENDIF
2395
2396
2397          CASE ( 1 )
2398
2399             t_surface  => t_surface_2; t_surface_p  => t_surface_1
2400             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
2401             IF ( humidity )  THEN
2402                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
2403                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
2404             ENDIF
2405
2406       END SELECT
2407#endif
2408
2409    END SUBROUTINE lsm_swap_timelevel
2410
2411
2412
2413
2414!------------------------------------------------------------------------------!
2415!
2416! Description:
2417! ------------
2418!> Subroutine for averaging 3D data
2419!------------------------------------------------------------------------------!
2420SUBROUTINE lsm_3d_data_averaging( mode, variable )
2421 
2422
2423    USE control_parameters
2424
2425    USE indices
2426
2427    USE kinds
2428
2429    IMPLICIT NONE
2430
2431    CHARACTER (LEN=*) ::  mode    !<
2432    CHARACTER (LEN=*) :: variable !<
2433
2434    INTEGER(iwp) ::  i !<
2435    INTEGER(iwp) ::  j !<
2436    INTEGER(iwp) ::  k !<
2437
2438    IF ( mode == 'allocate' )  THEN
2439
2440       SELECT CASE ( TRIM( variable ) )
2441
2442             CASE ( 'c_liq*' )
2443                IF ( .NOT. ALLOCATED( c_liq_av ) )  THEN
2444                   ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )
2445                ENDIF
2446                c_liq_av = 0.0_wp
2447
2448             CASE ( 'c_soil*' )
2449                IF ( .NOT. ALLOCATED( c_soil_av ) )  THEN
2450                   ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )
2451                ENDIF
2452                c_soil_av = 0.0_wp
2453
2454             CASE ( 'c_veg*' )
2455                IF ( .NOT. ALLOCATED( c_veg_av ) )  THEN
2456                   ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )
2457                ENDIF
2458                c_veg_av = 0.0_wp
2459
2460             CASE ( 'ghf_eb*' )
2461                IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
2462                   ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
2463                ENDIF
2464                ghf_eb_av = 0.0_wp
2465
2466             CASE ( 'lai*' )
2467                IF ( .NOT. ALLOCATED( lai_av ) )  THEN
2468                   ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )
2469                ENDIF
2470                lai_av = 0.0_wp
2471
2472             CASE ( 'm_liq_eb*' )
2473                IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
2474                   ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
2475                ENDIF
2476                m_liq_eb_av = 0.0_wp
2477
2478             CASE ( 'm_soil' )
2479                IF ( .NOT. ALLOCATED( m_soil_av ) )  THEN
2480                   ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
2481                ENDIF
2482                m_soil_av = 0.0_wp
2483
2484             CASE ( 'qsws_eb*' )
2485                IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
2486                   ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
2487                ENDIF
2488                qsws_eb_av = 0.0_wp
2489
2490             CASE ( 'qsws_liq_eb*' )
2491                IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
2492                   ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
2493                ENDIF
2494                qsws_liq_eb_av = 0.0_wp
2495
2496             CASE ( 'qsws_soil_eb*' )
2497                IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
2498                   ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
2499                ENDIF
2500                qsws_soil_eb_av = 0.0_wp
2501
2502             CASE ( 'qsws_veg_eb*' )
2503                IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
2504                   ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
2505                ENDIF
2506                qsws_veg_eb_av = 0.0_wp
2507
2508             CASE ( 'r_a*' )
2509                IF ( .NOT. ALLOCATED( r_a_av ) )  THEN
2510                   ALLOCATE( r_a_av(nysg:nyng,nxlg:nxrg) )
2511                ENDIF
2512                r_a_av = 0.0_wp
2513
2514             CASE ( 'r_s*' )
2515                IF ( .NOT. ALLOCATED( r_s_av ) )  THEN
2516                   ALLOCATE( r_s_av(nysg:nyng,nxlg:nxrg) )
2517                ENDIF
2518                r_s_av = 0.0_wp
2519
2520             CASE ( 'shf_eb*' )
2521                IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
2522                   ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
2523                ENDIF
2524                shf_eb_av = 0.0_wp
2525
2526             CASE ( 't_soil' )
2527                IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
2528                   ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
2529                ENDIF
2530                t_soil_av = 0.0_wp
2531
2532          CASE DEFAULT
2533             CONTINUE
2534
2535       END SELECT
2536
2537    ELSEIF ( mode == 'sum' )  THEN
2538
2539       SELECT CASE ( TRIM( variable ) )
2540
2541          CASE ( 'c_liq*' )
2542             DO  i = nxlg, nxrg
2543                DO  j = nysg, nyng
2544                   c_liq_av(j,i) = c_liq_av(j,i) + c_liq(j,i)
2545                ENDDO
2546             ENDDO
2547
2548          CASE ( 'c_soil*' )
2549             DO  i = nxlg, nxrg
2550                DO  j = nysg, nyng
2551                   c_soil_av(j,i) = c_soil_av(j,i) + (1.0 - c_veg(j,i))
2552                ENDDO
2553             ENDDO
2554
2555          CASE ( 'c_veg*' )
2556             DO  i = nxlg, nxrg
2557                DO  j = nysg, nyng
2558                   c_veg_av(j,i) = c_veg_av(j,i) + c_veg(j,i)
2559                ENDDO
2560             ENDDO
2561
2562          CASE ( 'ghf_eb*' )
2563             DO  i = nxlg, nxrg
2564                DO  j = nysg, nyng
2565                   ghf_eb_av(j,i) = ghf_eb_av(j,i) + ghf_eb(j,i)
2566                ENDDO
2567             ENDDO
2568
2569          CASE ( 'lai*' )
2570             DO  i = nxlg, nxrg
2571                DO  j = nysg, nyng
2572                   lai_av(j,i) = lai_av(j,i) + lai(j,i)
2573                ENDDO
2574             ENDDO
2575
2576          CASE ( 'm_liq_eb*' )
2577             DO  i = nxlg, nxrg
2578                DO  j = nysg, nyng
2579                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) + m_liq_eb(j,i)
2580                ENDDO
2581             ENDDO
2582
2583          CASE ( 'm_soil' )
2584             DO  i = nxlg, nxrg
2585                DO  j = nysg, nyng
2586                   DO  k = nzb_soil, nzt_soil
2587                      m_soil_av(k,j,i) = m_soil_av(k,j,i) + m_soil(k,j,i)
2588                   ENDDO
2589                ENDDO
2590             ENDDO
2591
2592          CASE ( 'qsws_eb*' )
2593             DO  i = nxlg, nxrg
2594                DO  j = nysg, nyng
2595                   qsws_eb_av(j,i) = qsws_eb_av(j,i) + qsws_eb(j,i)
2596                ENDDO
2597             ENDDO
2598
2599          CASE ( 'qsws_liq_eb*' )
2600             DO  i = nxlg, nxrg
2601                DO  j = nysg, nyng
2602                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) + qsws_liq_eb(j,i)
2603                ENDDO
2604             ENDDO
2605
2606          CASE ( 'qsws_soil_eb*' )
2607             DO  i = nxlg, nxrg
2608                DO  j = nysg, nyng
2609                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) + qsws_soil_eb(j,i)
2610                ENDDO
2611             ENDDO
2612
2613          CASE ( 'qsws_veg_eb*' )
2614             DO  i = nxlg, nxrg
2615                DO  j = nysg, nyng
2616                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) + qsws_veg_eb(j,i)
2617                ENDDO
2618             ENDDO
2619
2620          CASE ( 'r_a*' )
2621             DO  i = nxlg, nxrg
2622                DO  j = nysg, nyng
2623                   r_a_av(j,i) = r_a_av(j,i) + r_a(j,i)
2624                ENDDO
2625             ENDDO
2626
2627          CASE ( 'r_s*' )
2628             DO  i = nxlg, nxrg
2629                DO  j = nysg, nyng
2630                   r_s_av(j,i) = r_s_av(j,i) + r_s(j,i)
2631                ENDDO
2632             ENDDO
2633
2634          CASE ( 'shf_eb*' )
2635             DO  i = nxlg, nxrg
2636                DO  j = nysg, nyng
2637                   shf_eb_av(j,i) = shf_eb_av(j,i) + shf_eb(j,i)
2638                ENDDO
2639             ENDDO
2640
2641          CASE ( 't_soil' )
2642             DO  i = nxlg, nxrg
2643                DO  j = nysg, nyng
2644                   DO  k = nzb_soil, nzt_soil
2645                      t_soil_av(k,j,i) = t_soil_av(k,j,i) + t_soil(k,j,i)
2646                   ENDDO
2647                ENDDO
2648             ENDDO
2649
2650          CASE DEFAULT
2651             CONTINUE
2652
2653       END SELECT
2654
2655    ELSEIF ( mode == 'average' )  THEN
2656
2657       SELECT CASE ( TRIM( variable ) )
2658
2659          CASE ( 'c_liq*' )
2660             DO  i = nxlg, nxrg
2661                DO  j = nysg, nyng
2662                   c_liq_av(j,i) = c_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
2663                ENDDO
2664             ENDDO
2665
2666          CASE ( 'c_soil*' )
2667             DO  i = nxlg, nxrg
2668                DO  j = nysg, nyng
2669                   c_soil_av(j,i) = c_soil_av(j,i) / REAL( average_count_3d, KIND=wp )
2670                ENDDO
2671             ENDDO
2672
2673          CASE ( 'c_veg*' )
2674             DO  i = nxlg, nxrg
2675                DO  j = nysg, nyng
2676                   c_veg_av(j,i) = c_veg_av(j,i) / REAL( average_count_3d, KIND=wp )
2677                ENDDO
2678             ENDDO
2679
2680          CASE ( 'ghf_eb*' )
2681             DO  i = nxlg, nxrg
2682                DO  j = nysg, nyng
2683                   ghf_eb_av(j,i) = ghf_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
2684                ENDDO
2685             ENDDO
2686
2687         CASE ( 'lai*' )
2688             DO  i = nxlg, nxrg
2689                DO  j = nysg, nyng
2690                   lai_av(j,i) = lai_av(j,i) / REAL( average_count_3d, KIND=wp )
2691                ENDDO
2692             ENDDO
2693
2694          CASE ( 'm_liq_eb*' )
2695             DO  i = nxlg, nxrg
2696                DO  j = nysg, nyng
2697                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
2698                ENDDO
2699             ENDDO
2700
2701          CASE ( 'm_soil' )
2702             DO  i = nxlg, nxrg
2703                DO  j = nysg, nyng
2704                   DO  k = nzb_soil, nzt_soil
2705                      m_soil_av(k,j,i) = m_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2706                   ENDDO
2707                ENDDO
2708             ENDDO
2709
2710          CASE ( 'qsws_eb*' )
2711             DO  i = nxlg, nxrg
2712                DO  j = nysg, nyng
2713                   qsws_eb_av(j,i) = qsws_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
2714                ENDDO
2715             ENDDO
2716
2717          CASE ( 'qsws_liq_eb*' )
2718             DO  i = nxlg, nxrg
2719                DO  j = nysg, nyng
2720                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
2721                ENDDO
2722             ENDDO
2723
2724          CASE ( 'qsws_soil_eb*' )
2725             DO  i = nxlg, nxrg
2726                DO  j = nysg, nyng
2727                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
2728                ENDDO
2729             ENDDO
2730
2731          CASE ( 'qsws_veg_eb*' )
2732             DO  i = nxlg, nxrg
2733                DO  j = nysg, nyng
2734                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
2735                ENDDO
2736             ENDDO
2737
2738          CASE ( 'r_a*' )
2739             DO  i = nxlg, nxrg
2740                DO  j = nysg, nyng
2741                   r_a_av(j,i) = r_a_av(j,i) / REAL( average_count_3d, KIND=wp )
2742                ENDDO
2743             ENDDO
2744
2745          CASE ( 'r_s*' )
2746             DO  i = nxlg, nxrg
2747                DO  j = nysg, nyng
2748                   r_s_av(j,i) = r_s_av(j,i) / REAL( average_count_3d, KIND=wp )
2749                ENDDO
2750             ENDDO
2751
2752          CASE ( 't_soil' )
2753             DO  i = nxlg, nxrg
2754                DO  j = nysg, nyng
2755                   DO  k = nzb_soil, nzt_soil
2756                      t_soil_av(k,j,i) = t_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2757                   ENDDO
2758                ENDDO
2759             ENDDO
2760
2761       END SELECT
2762
2763    ENDIF
2764
2765END SUBROUTINE lsm_3d_data_averaging
2766
2767
2768!------------------------------------------------------------------------------!
2769!
2770! Description:
2771! ------------
2772!> Subroutine defining appropriate grid for netcdf variables.
2773!> It is called out from subroutine netcdf.
2774!------------------------------------------------------------------------------!
2775 SUBROUTINE lsm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
2776   
2777     IMPLICIT NONE
2778
2779     CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
2780     LOGICAL, INTENT(OUT)           ::  found       !<
2781     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
2782     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
2783     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
2784
2785     found  = .TRUE.
2786
2787!
2788!--  Check for the grid
2789     SELECT CASE ( TRIM( var ) )
2790
2791        CASE ( 'm_soil', 't_soil', 'm_soil_xy', 't_soil_xy', 'm_soil_xz',      &
2792               't_soil_xz', 'm_soil_yz', 't_soil_yz' )
2793           grid_x = 'x'
2794           grid_y = 'y'
2795           grid_z = 'zs'
2796
2797        CASE DEFAULT
2798           found  = .FALSE.
2799           grid_x = 'none'
2800           grid_y = 'none'
2801           grid_z = 'none'
2802     END SELECT
2803
2804 END SUBROUTINE lsm_define_netcdf_grid
2805
2806!------------------------------------------------------------------------------!
2807!
2808! Description:
2809! ------------
2810!> Subroutine defining 3D output variables
2811!------------------------------------------------------------------------------!
2812 SUBROUTINE lsm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
2813                                two_d, nzb_do, nzt_do )
2814 
2815    USE indices
2816
2817    USE kinds
2818
2819
2820    IMPLICIT NONE
2821
2822    CHARACTER (LEN=*) ::  grid     !<
2823    CHARACTER (LEN=*) ::  mode     !<
2824    CHARACTER (LEN=*) ::  variable !<
2825
2826    INTEGER(iwp) ::  av !<
2827    INTEGER(iwp) ::  i  !<
2828    INTEGER(iwp) ::  j  !<
2829    INTEGER(iwp) ::  k  !<
2830    INTEGER(iwp) ::  nzb_do  !<
2831    INTEGER(iwp) ::  nzt_do  !<
2832
2833    LOGICAL      ::  found !<
2834    LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
2835
2836    REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !<
2837
2838    found = .TRUE.
2839
2840    SELECT CASE ( TRIM( variable ) )
2841
2842
2843       CASE ( 'c_liq*_xy' )        ! 2d-array
2844          IF ( av == 0 )  THEN
2845             DO  i = nxlg, nxrg
2846                DO  j = nysg, nyng
2847                   local_pf(i,j,nzb+1) = c_liq(j,i) * c_veg(j,i)
2848                ENDDO
2849             ENDDO
2850          ELSE
2851             DO  i = nxlg, nxrg
2852                DO  j = nysg, nyng
2853                   local_pf(i,j,nzb+1) = c_liq_av(j,i)
2854                ENDDO
2855             ENDDO
2856          ENDIF
2857
2858          two_d = .TRUE.
2859          grid = 'zu1'
2860
2861       CASE ( 'c_soil*_xy' )        ! 2d-array
2862          IF ( av == 0 )  THEN
2863             DO  i = nxlg, nxrg
2864                DO  j = nysg, nyng
2865                   local_pf(i,j,nzb+1) = 1.0_wp - c_veg(j,i)
2866                ENDDO
2867             ENDDO
2868          ELSE
2869             DO  i = nxlg, nxrg
2870                DO  j = nysg, nyng
2871                   local_pf(i,j,nzb+1) = c_soil_av(j,i)
2872                ENDDO
2873             ENDDO
2874          ENDIF
2875
2876          two_d = .TRUE.
2877          grid = 'zu1'
2878
2879       CASE ( 'c_veg*_xy' )        ! 2d-array
2880          IF ( av == 0 )  THEN
2881             DO  i = nxlg, nxrg
2882                DO  j = nysg, nyng
2883                   local_pf(i,j,nzb+1) = c_veg(j,i)
2884                ENDDO
2885             ENDDO
2886          ELSE
2887             DO  i = nxlg, nxrg
2888                DO  j = nysg, nyng
2889                   local_pf(i,j,nzb+1) = c_veg_av(j,i)
2890                ENDDO
2891             ENDDO
2892          ENDIF
2893
2894          two_d = .TRUE.
2895          grid = 'zu1'
2896
2897       CASE ( 'ghf_eb*_xy' )        ! 2d-array
2898          IF ( av == 0 )  THEN
2899             DO  i = nxlg, nxrg
2900                DO  j = nysg, nyng
2901                   local_pf(i,j,nzb+1) = ghf_eb(j,i)
2902                ENDDO
2903             ENDDO
2904          ELSE
2905             DO  i = nxlg, nxrg
2906                DO  j = nysg, nyng
2907                   local_pf(i,j,nzb+1) = ghf_eb_av(j,i)
2908                ENDDO
2909             ENDDO
2910          ENDIF
2911
2912          two_d = .TRUE.
2913          grid = 'zu1'
2914
2915       CASE ( 'lai*_xy' )        ! 2d-array
2916          IF ( av == 0 )  THEN
2917             DO  i = nxlg, nxrg
2918                DO  j = nysg, nyng
2919                   local_pf(i,j,nzb+1) = lai(j,i)
2920                ENDDO
2921             ENDDO
2922          ELSE
2923             DO  i = nxlg, nxrg
2924                DO  j = nysg, nyng
2925                   local_pf(i,j,nzb+1) = lai_av(j,i)
2926                ENDDO
2927             ENDDO
2928          ENDIF
2929
2930          two_d = .TRUE.
2931          grid = 'zu1'
2932
2933       CASE ( 'm_liq_eb*_xy' )        ! 2d-array
2934          IF ( av == 0 )  THEN
2935             DO  i = nxlg, nxrg
2936                DO  j = nysg, nyng
2937                   local_pf(i,j,nzb+1) = m_liq_eb(j,i)
2938                ENDDO
2939             ENDDO
2940          ELSE
2941             DO  i = nxlg, nxrg
2942                DO  j = nysg, nyng
2943                   local_pf(i,j,nzb+1) = m_liq_eb_av(j,i)
2944                ENDDO
2945             ENDDO
2946          ENDIF
2947
2948          two_d = .TRUE.
2949          grid = 'zu1'
2950
2951       CASE ( 'm_soil_xy', 'm_soil_xz', 'm_soil_yz' )
2952          IF ( av == 0 )  THEN
2953             DO  i = nxlg, nxrg
2954                DO  j = nysg, nyng
2955                   DO k = nzb_soil, nzt_soil
2956                      local_pf(i,j,k) = m_soil(k,j,i)
2957                   ENDDO
2958                ENDDO
2959             ENDDO
2960          ELSE
2961             DO  i = nxlg, nxrg
2962                DO  j = nysg, nyng
2963                   DO k = nzb_soil, nzt_soil
2964                      local_pf(i,j,k) = m_soil_av(k,j,i)
2965                   ENDDO
2966                ENDDO
2967             ENDDO
2968          ENDIF
2969
2970          nzb_do = nzb_soil
2971          nzt_do = nzt_soil
2972
2973          IF ( mode == 'xy' ) grid = 'zs'
2974
2975       CASE ( 'qsws_eb*_xy' )        ! 2d-array
2976          IF ( av == 0 ) THEN
2977             DO  i = nxlg, nxrg
2978                DO  j = nysg, nyng
2979                   local_pf(i,j,nzb+1) =  qsws_eb(j,i)
2980                ENDDO
2981             ENDDO
2982          ELSE
2983             DO  i = nxlg, nxrg
2984                DO  j = nysg, nyng 
2985                   local_pf(i,j,nzb+1) =  qsws_eb_av(j,i)
2986                ENDDO
2987             ENDDO
2988          ENDIF
2989
2990          two_d = .TRUE.
2991          grid = 'zu1'
2992
2993       CASE ( 'qsws_liq_eb*_xy' )        ! 2d-array
2994          IF ( av == 0 ) THEN
2995             DO  i = nxlg, nxrg
2996                DO  j = nysg, nyng
2997                   local_pf(i,j,nzb+1) =  qsws_liq_eb(j,i)
2998                ENDDO
2999             ENDDO
3000          ELSE
3001             DO  i = nxlg, nxrg
3002                DO  j = nysg, nyng 
3003                   local_pf(i,j,nzb+1) =  qsws_liq_eb_av(j,i)
3004                ENDDO
3005             ENDDO
3006          ENDIF
3007
3008          two_d = .TRUE.
3009          grid = 'zu1'
3010
3011       CASE ( 'qsws_soil_eb*_xy' )        ! 2d-array
3012          IF ( av == 0 ) THEN
3013             DO  i = nxlg, nxrg
3014                DO  j = nysg, nyng
3015                   local_pf(i,j,nzb+1) =  qsws_soil_eb(j,i)
3016                ENDDO
3017             ENDDO
3018          ELSE
3019             DO  i = nxlg, nxrg
3020                DO  j = nysg, nyng 
3021                   local_pf(i,j,nzb+1) =  qsws_soil_eb_av(j,i)
3022                ENDDO
3023             ENDDO
3024          ENDIF
3025
3026          two_d = .TRUE.
3027          grid = 'zu1'
3028
3029       CASE ( 'qsws_veg_eb*_xy' )        ! 2d-array
3030          IF ( av == 0 ) THEN
3031             DO  i = nxlg, nxrg
3032                DO  j = nysg, nyng
3033                   local_pf(i,j,nzb+1) =  qsws_veg_eb(j,i)
3034                ENDDO
3035             ENDDO
3036          ELSE
3037             DO  i = nxlg, nxrg
3038                DO  j = nysg, nyng 
3039                   local_pf(i,j,nzb+1) =  qsws_veg_eb_av(j,i)
3040                ENDDO
3041             ENDDO
3042          ENDIF
3043
3044          two_d = .TRUE.
3045          grid = 'zu1'
3046
3047
3048       CASE ( 'r_a*_xy' )        ! 2d-array
3049          IF ( av == 0 )  THEN
3050             DO  i = nxlg, nxrg
3051                DO  j = nysg, nyng
3052                   local_pf(i,j,nzb+1) = r_a(j,i)
3053                ENDDO
3054             ENDDO
3055          ELSE
3056             DO  i = nxlg, nxrg
3057                DO  j = nysg, nyng
3058                   local_pf(i,j,nzb+1) = r_a_av(j,i)
3059                ENDDO
3060             ENDDO
3061          ENDIF
3062
3063          two_d = .TRUE.
3064          grid = 'zu1'
3065
3066       CASE ( 'r_s*_xy' )        ! 2d-array
3067          IF ( av == 0 )  THEN
3068             DO  i = nxlg, nxrg
3069                DO  j = nysg, nyng
3070                   local_pf(i,j,nzb+1) = r_s(j,i)
3071                ENDDO
3072             ENDDO
3073          ELSE
3074             DO  i = nxlg, nxrg
3075                DO  j = nysg, nyng
3076                   local_pf(i,j,nzb+1) = r_s_av(j,i)
3077                ENDDO
3078             ENDDO
3079          ENDIF
3080
3081          two_d = .TRUE.
3082          grid = 'zu1'
3083
3084       CASE ( 'shf_eb*_xy' )        ! 2d-array
3085          IF ( av == 0 ) THEN
3086             DO  i = nxlg, nxrg
3087                DO  j = nysg, nyng
3088                   local_pf(i,j,nzb+1) =  shf_eb(j,i)
3089                ENDDO
3090             ENDDO
3091          ELSE
3092             DO  i = nxlg, nxrg
3093                DO  j = nysg, nyng
3094                   local_pf(i,j,nzb+1) =  shf_eb_av(j,i)
3095                ENDDO
3096             ENDDO
3097          ENDIF
3098
3099          two_d = .TRUE.
3100          grid = 'zu1'
3101
3102       CASE ( 't_soil_xy', 't_soil_xz', 't_soil_yz' )
3103          IF ( av == 0 )  THEN
3104             DO  i = nxlg, nxrg
3105                DO  j = nysg, nyng
3106                   DO k = nzb_soil, nzt_soil
3107                      local_pf(i,j,k) = t_soil(k,j,i)
3108                   ENDDO
3109                ENDDO
3110             ENDDO
3111          ELSE
3112             DO  i = nxlg, nxrg
3113                DO  j = nysg, nyng
3114                   DO k = nzb_soil, nzt_soil
3115                      local_pf(i,j,k) = t_soil_av(k,j,i)
3116                   ENDDO
3117                ENDDO
3118             ENDDO
3119          ENDIF
3120
3121          nzb_do = nzb_soil
3122          nzt_do = nzt_soil
3123
3124          IF ( mode == 'xy' )  grid = 'zs'
3125
3126       CASE DEFAULT
3127          found = .FALSE.
3128          grid  = 'none'
3129
3130    END SELECT
3131 
3132 END SUBROUTINE lsm_data_output_2d
3133
3134
3135!------------------------------------------------------------------------------!
3136!
3137! Description:
3138! ------------
3139!> Subroutine defining 3D output variables
3140!------------------------------------------------------------------------------!
3141 SUBROUTINE lsm_data_output_3d( av, variable, found, local_pf )
3142 
3143
3144    USE indices
3145
3146    USE kinds
3147
3148
3149    IMPLICIT NONE
3150
3151    CHARACTER (LEN=*) ::  variable !<
3152
3153    INTEGER(iwp) ::  av    !<
3154    INTEGER(iwp) ::  i     !<
3155    INTEGER(iwp) ::  j     !<
3156    INTEGER(iwp) ::  k     !<
3157
3158    LOGICAL      ::  found !<
3159
3160    REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb_soil:nzt_soil) ::  local_pf !<
3161
3162
3163    found = .TRUE.
3164
3165
3166    SELECT CASE ( TRIM( variable ) )
3167
3168
3169      CASE ( 'm_soil' )
3170
3171         IF ( av == 0 )  THEN
3172            DO  i = nxlg, nxrg
3173               DO  j = nysg, nyng
3174                  DO  k = nzb_soil, nzt_soil
3175                     local_pf(i,j,k) = m_soil(k,j,i)
3176                  ENDDO
3177               ENDDO
3178            ENDDO
3179         ELSE
3180            DO  i = nxlg, nxrg
3181               DO  j = nysg, nyng
3182                  DO  k = nzb_soil, nzt_soil
3183                     local_pf(i,j,k) = m_soil_av(k,j,i)
3184                  ENDDO
3185               ENDDO
3186            ENDDO
3187         ENDIF
3188
3189      CASE ( 't_soil' )
3190
3191         IF ( av == 0 )  THEN
3192            DO  i = nxlg, nxrg
3193               DO  j = nysg, nyng
3194                  DO  k = nzb_soil, nzt_soil
3195                     local_pf(i,j,k) = t_soil(k,j,i)
3196                  ENDDO
3197               ENDDO
3198            ENDDO
3199         ELSE
3200            DO  i = nxlg, nxrg
3201               DO  j = nysg, nyng
3202                  DO  k = nzb_soil, nzt_soil
3203                     local_pf(i,j,k) = t_soil_av(k,j,i)
3204                  ENDDO
3205               ENDDO
3206            ENDDO
3207         ENDIF
3208
3209
3210       CASE DEFAULT
3211          found = .FALSE.
3212
3213    END SELECT
3214
3215
3216 END SUBROUTINE lsm_data_output_3d
3217
3218
3219!------------------------------------------------------------------------------!
3220!
3221! Description:
3222! ------------
3223!> Write restart data for land surface model
3224!------------------------------------------------------------------------------!
3225 SUBROUTINE lsm_last_actions
3226 
3227
3228    USE control_parameters
3229       
3230    USE kinds
3231
3232    IMPLICIT NONE
3233
3234    IF ( write_binary(1:4) == 'true' )  THEN
3235       IF ( ALLOCATED( c_liq_av ) )  THEN
3236          WRITE ( 14 )  'c_liq_av            ';  WRITE ( 14 ) c_liq_av
3237       ENDIF
3238       IF ( ALLOCATED( c_soil_av ) )  THEN
3239          WRITE ( 14 )  'c_soil_av           ';  WRITE ( 14 ) c_soil_av
3240       ENDIF
3241       IF ( ALLOCATED( c_veg_av ) )  THEN
3242          WRITE ( 14 )  'c_veg_av            ';  WRITE ( 14 ) c_veg_av
3243       ENDIF
3244       IF ( ALLOCATED( ghf_eb_av ) )  THEN
3245          WRITE ( 14 )  'ghf_eb_av           ';  WRITE ( 14 )  ghf_eb_av
3246       ENDIF
3247       IF ( ALLOCATED( lai_av ) )  THEN
3248          WRITE ( 14 )  'lai_av              ';  WRITE ( 14 )  lai_av
3249       ENDIF
3250       WRITE ( 14 )  'm_liq_eb            ';  WRITE ( 14 )  m_liq_eb
3251       IF ( ALLOCATED( m_liq_eb_av ) )  THEN
3252          WRITE ( 14 )  'm_liq_eb_av         ';  WRITE ( 14 )  m_liq_eb_av
3253       ENDIF
3254       WRITE ( 14 )  'm_soil              ';  WRITE ( 14 )  m_soil
3255       IF ( ALLOCATED( m_soil_av ) )  THEN
3256          WRITE ( 14 )  'm_soil_av           ';  WRITE ( 14 )  m_soil_av
3257       ENDIF
3258       IF ( ALLOCATED( qsws_eb_av ) )  THEN
3259          WRITE ( 14 )  'qsws_eb_av          ';  WRITE ( 14 )  qsws_eb_av
3260       ENDIF   
3261       IF ( ALLOCATED( qsws_liq_eb_av ) )  THEN
3262          WRITE ( 14 )  'qsws_liq_eb_av      ';  WRITE ( 14 )  qsws_liq_eb_av
3263       ENDIF 
3264       IF ( ALLOCATED( qsws_soil_eb_av ) )  THEN
3265          WRITE ( 14 )  'qsws_soil_eb_av     ';  WRITE ( 14 )  qsws_soil_eb_av
3266       ENDIF
3267       IF ( ALLOCATED( qsws_veg_eb_av ) )  THEN
3268          WRITE ( 14 )  'qsws_veg_eb_av      ';  WRITE ( 14 )  qsws_veg_eb_av
3269       ENDIF
3270       IF ( ALLOCATED( shf_eb_av ) )  THEN
3271          WRITE ( 14 )  'shf_eb_av           ';  WRITE ( 14 )  shf_eb_av
3272       ENDIF
3273       WRITE ( 14 )  't_soil              ';  WRITE ( 14 )  t_soil
3274       IF ( ALLOCATED( t_soil_av ) )  THEN
3275          WRITE ( 14 )  't_soil_av           ';  WRITE ( 14 )  t_soil_av
3276       ENDIF
3277
3278       WRITE ( 14 )  '*** end lsm ***     '
3279
3280    ENDIF
3281
3282 END SUBROUTINE lsm_last_actions
3283
3284
3285SUBROUTINE lsm_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file,   &
3286                                     nynfa, nyn_on_file, nysfa, nys_on_file,   &
3287                                     offset_xa, offset_ya, overlap_count,      &
3288                                     tmp_2d )
3289 
3290
3291    USE control_parameters
3292       
3293    USE indices
3294   
3295    USE kinds
3296   
3297    USE pegrid
3298
3299    IMPLICIT NONE
3300
3301    CHARACTER (LEN=20) :: field_char   !<
3302
3303    INTEGER(iwp) ::  i               !<
3304    INTEGER(iwp) ::  k               !<
3305    INTEGER(iwp) ::  nxlc            !<
3306    INTEGER(iwp) ::  nxlf            !<
3307    INTEGER(iwp) ::  nxl_on_file     !<
3308    INTEGER(iwp) ::  nxrc            !<
3309    INTEGER(iwp) ::  nxrf            !<
3310    INTEGER(iwp) ::  nxr_on_file     !<
3311    INTEGER(iwp) ::  nync            !<
3312    INTEGER(iwp) ::  nynf            !<
3313    INTEGER(iwp) ::  nyn_on_file     !<
3314    INTEGER(iwp) ::  nysc            !<
3315    INTEGER(iwp) ::  nysf            !<
3316    INTEGER(iwp) ::  nys_on_file     !<
3317    INTEGER(iwp) ::  overlap_count   !<
3318
3319    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxlfa       !<
3320    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxrfa       !<
3321    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nynfa       !<
3322    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nysfa       !<
3323    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_xa   !<
3324    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !<
3325
3326    REAL(wp),                                                                  &
3327       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3328          tmp_2d   !<
3329
3330    REAL(wp),                                                                  &
3331       DIMENSION(nzb_soil:nzt_soil+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3332          tmp_3d   !<
3333
3334    REAL(wp),                                                                  &
3335       DIMENSION(nzb_soil:nzt_soil,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3336          tmp_3d2   !<
3337
3338
3339   IF ( initializing_actions == 'read_restart_data' )  THEN
3340      READ ( 13 )  field_char
3341
3342      DO  WHILE ( TRIM( field_char ) /= '*** end lsm ***' )
3343
3344         DO  k = 1, overlap_count
3345
3346            nxlf = nxlfa(i,k)
3347            nxlc = nxlfa(i,k) + offset_xa(i,k)
3348            nxrf = nxrfa(i,k)
3349            nxrc = nxrfa(i,k) + offset_xa(i,k)
3350            nysf = nysfa(i,k)
3351            nysc = nysfa(i,k) + offset_ya(i,k)
3352            nynf = nynfa(i,k)
3353            nync = nynfa(i,k) + offset_ya(i,k)
3354
3355
3356            SELECT CASE ( TRIM( field_char ) )
3357
3358                CASE ( 'c_liq_av' )
3359                   IF ( .NOT. ALLOCATED( c_liq_av ) )  THEN
3360                      ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )
3361                   ENDIF
3362                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3363                   c_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3364                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3365
3366                CASE ( 'c_soil_av' )
3367                   IF ( .NOT. ALLOCATED( c_soil_av ) )  THEN
3368                      ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )
3369                   ENDIF
3370                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3371                   c_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3372                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3373
3374                CASE ( 'c_veg_av' )
3375                   IF ( .NOT. ALLOCATED( c_veg_av ) )  THEN
3376                      ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )
3377                   ENDIF
3378                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3379                   c_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3380                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3381
3382                CASE ( 'ghf_eb_av' )
3383                   IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
3384                      ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
3385                   ENDIF
3386                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3387                   ghf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3388                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3389
3390                CASE ( 'm_liq_eb' )
3391                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3392                   m_liq_eb(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =        &
3393                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3394
3395                CASE ( 'lai_av' )
3396                   IF ( .NOT. ALLOCATED( lai_av ) )  THEN
3397                      ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )
3398                   ENDIF
3399                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3400                   lai_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3401                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3402
3403                CASE ( 'm_liq_eb_av' )
3404                   IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
3405                      ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
3406                   ENDIF
3407                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3408                   m_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
3409                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3410
3411                CASE ( 'm_soil' )
3412                   IF ( k == 1 )  READ ( 13 )  tmp_3d2(:,:,:)
3413                   m_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
3414                          tmp_3d2(nzb_soil:nzt_soil,nysf-nbgp:nynf             &
3415                          +nbgp,nxlf-nbgp:nxrf+nbgp)
3416
3417                CASE ( 'm_soil_av' )
3418                   IF ( .NOT. ALLOCATED( m_soil_av ) )  THEN
3419                      ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
3420                   ENDIF
3421                   IF ( k == 1 )  READ ( 13 )  tmp_3d2(:,:,:)
3422                   m_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
3423                                    tmp_3d2(nzb_soil:nzt_soil,nysf             &
3424                                    -nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3425
3426                CASE ( 'qsws_eb_av' )
3427                   IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
3428                      ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
3429                   ENDIF 
3430                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3431                   qsws_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
3432                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3433
3434                CASE ( 'qsws_liq_eb_av' )
3435                   IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
3436                      ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
3437                   ENDIF 
3438                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3439                   qsws_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
3440                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3441                CASE ( 'qsws_soil_eb_av' )
3442                   IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
3443                      ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
3444                   ENDIF 
3445                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3446                   qsws_soil_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
3447                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3448
3449                CASE ( 'qsws_veg_eb_av' )
3450                   IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
3451                      ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
3452                   ENDIF 
3453                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3454                   qsws_veg_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =  &
3455                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3456
3457                CASE ( 'shf_eb_av' )
3458                   IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
3459                      ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
3460                   ENDIF
3461                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3462                   shf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
3463                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3464
3465                CASE ( 't_soil' )
3466                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3467                   t_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
3468                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,               &
3469                                                nxlf-nbgp:nxrf+nbgp)
3470
3471                CASE ( 't_soil_av' )
3472                   IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
3473                      ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
3474                   ENDIF
3475                   IF ( k == 1 )  READ ( 13 )  tmp_3d2(:,:,:)
3476                   t_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
3477                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,             &
3478                                    nxlf-nbgp:nxrf+nbgp)
3479
3480
3481               CASE DEFAULT
3482                  WRITE( message_string, * ) 'unknown variable named "',       &
3483                                        TRIM( field_char ), '" found in',      &
3484                                        '&data from prior run on PE ', myid
3485                  CALL message( 'lsm_read_restart_data', 'PA0441', 1, 2, 0, 6, &
3486                                 0 )
3487
3488            END SELECT
3489
3490         ENDDO
3491
3492         READ ( 13 )  field_char
3493
3494      ENDDO
3495   ENDIF
3496
3497 END SUBROUTINE lsm_read_restart_data
3498
3499!------------------------------------------------------------------------------!
3500! Description:
3501! ------------
3502!> Calculation of roughness length for open water (lakes, ocean). The
3503!> parameterization follows Charnock (1955). Two different implementations
3504!> are available: as in ECMWF-IFS (Beljaars 1994) or as in FLake (Subin et al.
3505!> 2012)
3506!------------------------------------------------------------------------------!
3507    SUBROUTINE calc_z0_water_surface
3508
3509       USE control_parameters,                                                 &
3510           ONLY: g, kappa, molecular_viscosity
3511
3512       IMPLICIT NONE
3513
3514       INTEGER :: i  !< running index
3515       INTEGER :: j  !< running index
3516
3517       REAL(wp), PARAMETER :: alpha_ch  = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF
3518!       REAL(wp), PARAMETER :: pr_number = 0.71_wp !< molecular Prandtl number in the Charnock parameterization (differs from prandtl_number)
3519!       REAL(wp), PARAMETER :: sc_number = 0.66_wp !< molecular Schmidt number in the Charnock parameterization
3520!       REAL(wp) :: re_0 !< near-surface roughness Reynolds number
3521
3522
3523       DO  i = nxlg, nxrg   
3524          DO  j = nysg, nyng
3525             IF ( water_surface(j,i) )  THEN
3526
3527!
3528!--             Disabled: FLake parameterization. Ideally, the Charnock
3529!--             coefficient should depend on the water depth and the fetch
3530!--             length
3531!                re_0 = z0(j,i) * us(j,i) / molecular_viscosity
3532!       
3533!                z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
3534!                              alpha_ch * us(j,i) / g )
3535!
3536!                z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) )
3537!                z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) )
3538
3539!
3540!--              Set minimum roughness length for u* > 0.2
3541!                IF ( us(j,i) > 0.2_wp )  THEN
3542!                   z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) )
3543!                   z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) )
3544!                ENDIF
3545
3546!
3547!--             ECMWF IFS model parameterization after Beljaars (1994). At low
3548!--             wind speed, the sea surface becomes aerodynamically smooth and
3549!--             the roughness scales with the viscosity. At high wind speed, the
3550!--             Charnock relation is used.
3551                z0(j,i) =   ( 0.11_wp * molecular_viscosity / us(j,i) )        &
3552                          + ( alpha_ch * us(j,i)**2 / g )
3553
3554                z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i)
3555                z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i)
3556
3557             ENDIF
3558          ENDDO
3559       ENDDO
3560
3561    END SUBROUTINE calc_z0_water_surface
3562
3563
3564!------------------------------------------------------------------------------!
3565! Description:
3566! ------------
3567!> Calculation of specific humidity of the skin layer (surface). It is assumend
3568!> that the skin is always saturated.
3569!------------------------------------------------------------------------------!
3570    SUBROUTINE calc_q_surface
3571
3572       IMPLICIT NONE
3573
3574       INTEGER :: i              !< running index
3575       INTEGER :: j              !< running index
3576       INTEGER :: k              !< running index
3577
3578       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
3579
3580       DO  i = nxlg, nxrg   
3581          DO  j = nysg, nyng
3582             k = nzb_s_inner(j,i)
3583
3584!
3585!--          Calculate water vapour pressure at saturation
3586             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface_p(j,i)   &
3587                   - 273.16_wp ) / ( t_surface_p(j,i) - 35.86_wp ) )
3588
3589!
3590!--          Calculate specific humidity at saturation
3591             q_s = 0.622_wp * e_s / surface_pressure
3592
3593             resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
3594
3595!
3596!--          Calculate specific humidity at surface
3597             IF ( cloud_physics )  THEN
3598                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
3599                             * ( q(k+1,j,i) - ql(k+1,j,i) )
3600             ELSE
3601                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
3602                             * q(k+1,j,i)
3603             ENDIF
3604
3605!
3606!--          Update virtual potential temperature
3607             vpt(k,j,i) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
3608
3609          ENDDO
3610       ENDDO
3611
3612    END SUBROUTINE calc_q_surface
3613
3614
3615 END MODULE land_surface_model_mod
Note: See TracBrowser for help on using the repository browser.