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

Last change on this file since 2068 was 2032, checked in by knoop, 8 years ago

last commit documented

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