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

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

further modularization of land surface model (2D/3D output and restart data)

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