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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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