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

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

last commit documented

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