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

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

last commit documented

  • Property svn:keywords set to Id
File size: 110.1 KB
Line 
1!> @file land_surface_model_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: land_surface_model_mod.f90 1967 2016-07-18 11:56:58Z maronga $
26!
27! 1966 2016-07-18 11:54:18Z maronga
28! Bugfix: calculation of m_total in soil model was not set to zero at model start
29!
30! 1949 2016-06-17 07:19:16Z maronga
31! Bugfix: calculation of qsws_soil_eb with precipitation = .TRUE. gave
32! qsws_soil_eb = 0 due to a typo
33!
34! 1856 2016-04-13 12:56:17Z maronga
35! Bugfix: for water surfaces, the initial water surface temperature is set equal
36! to the intital skin temperature. Moreover, the minimum value of r_a is now
37! 1.0 to avoid too large fluxes at the first model time step
38!
39! 1849 2016-04-08 11:33:18Z hoffmann
40! prr moved to arrays_3d
41!
42! 1826 2016-04-07 12:01:39Z maronga
43! Cleanup after modularization
44!
45! 1817 2016-04-06 15:44:20Z maronga
46! Added interface for lsm_init_arrays. Added subroutines for check_parameters,
47! header, and parin. Renamed some subroutines.
48!
49! 1788 2016-03-10 11:01:04Z maronga
50! Bugfix: calculate lambda_surface based on temperature gradient between skin
51! layer and soil layer instead of Obukhov length
52! Changed: moved calculation of surface specific humidity to energy balance solver
53! New: water surfaces are available by using a fixed sea surface temperature.
54! The roughness lengths are calculated dynamically using the Charnock
55! parameterization. This involves the new roughness length for moisture z0q.
56! New: modified solution of the energy balance solver and soil model for
57! paved surfaces (i.e. asphalt concrete).
58! Syntax layout improved.
59! Changed: parameter dewfall removed.
60!
61! 1783 2016-03-06 18:36:17Z raasch
62! netcdf variables moved to netcdf module
63!
64! 1757 2016-02-22 15:49:32Z maronga
65! Bugfix: set tm_soil_m to zero after allocation. Added parameter
66! unscheduled_radiation_calls to control calls of the radiation model based on
67! the skin temperature change during one time step (preliminary version). Set
68! qsws_soil_eb to zero at model start (previously set to qsws_eb). Removed MAX
69! function as it cannot be vectorized.
70!
71! 1709 2015-11-04 14:47:01Z maronga
72! Renamed pt_1 and qv_1 to pt1 and qv1.
73! Bugfix: set initial values for t_surface_p in case of restart runs
74! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky'
75! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky'
76! Added todo action
77!
78! 1697 2015-10-28 17:14:10Z raasch
79! bugfix: misplaced cpp-directive
80!
81! 1695 2015-10-27 10:03:11Z maronga
82! Bugfix: REAL constants provided with KIND-attribute in call of
83! Replaced rif with ol
84!
85! 1691 2015-10-26 16:17:44Z maronga
86! Added skip_time_do_lsm to allow for spin-ups without LSM. Various bugfixes:
87! Soil temperatures are now defined at the edges of the layers, calculation of
88! shb_eb corrected, prognostic equation for skin temperature corrected. Surface
89! fluxes are now directly transfered to atmosphere
90!
91! 1682 2015-10-07 23:56:08Z knoop
92! Code annotations made doxygen readable
93!
94! 1590 2015-05-08 13:56:27Z maronga
95! Bugfix: definition of character strings requires same length for all elements
96!
97! 1585 2015-04-30 07:05:52Z maronga
98! Modifications for RRTMG. Changed tables to PARAMETER type.
99!
100! 1571 2015-03-12 16:12:49Z maronga
101! Removed upper-case variable names. Corrected distribution of precipitation to
102! the liquid water reservoir and the bare soil fractions.
103!
104! 1555 2015-03-04 17:44:27Z maronga
105! Added output of r_a and r_s
106!
107! 1553 2015-03-03 17:33:54Z maronga
108! Improved better treatment of roughness lengths. Added default soil temperature
109! profile
110!
111! 1551 2015-03-03 14:18:16Z maronga
112! Flux calculation is now done in prandtl_fluxes. Added support for data output.
113! Vertical indices have been replaced. Restart runs are now possible. Some
114! variables have beem renamed. Bugfix in the prognostic equation for the surface
115! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
116! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
117! the hydraulic conductivity. Bugfix for root fraction and extraction
118! calculation
119!
120! intrinsic function MAX and MIN
121!
122! 1500 2014-12-03 17:42:41Z maronga
123! Corrected calculation of aerodynamic resistance (r_a).
124! Precipitation is now added to liquid water reservoir using LE_liq.
125! Added support for dry runs.
126!
127! 1496 2014-12-02 17:25:50Z maronga
128! Initial revision
129!
130!
131! Description:
132! ------------
133!> Land surface model, consisting of a solver for the energy balance at the
134!> surface and a four layer soil scheme. The scheme is similar to the TESSEL
135!> scheme implemented in the ECMWF IFS model, with modifications according to
136!> H-TESSEL. The implementation is based on the formulation implemented in the
137!> DALES and UCLA-LES models.
138!>
139!> @todo Consider partial absorption of the net shortwave radiation by the
140!>       skin layer.
141!> @todo Improve surface water parameterization
142!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
143!>       nzt_soil=3)).
144!> @todo Implement surface runoff model (required when performing long-term LES
145!>       with considerable precipitation.
146!> @todo Fix crashes with radiation_scheme == 'constant'
147!>
148!> @note No time step criterion is required as long as the soil layers do not
149!>       become too thin.
150!------------------------------------------------------------------------------!
151 MODULE land_surface_model_mod
152 
153    USE arrays_3d,                                                             &
154        ONLY:  hyp, ol, pt, pt_p, prr, q, q_p, ql, qsws, shf, ts, us, vpt, z0, &
155               z0h, z0q
156
157    USE cloud_parameters,                                                      &
158        ONLY:  cp, hyrho, l_d_cp, l_d_r, l_v, pt_d_t, rho_l, r_d, r_v
159
160    USE control_parameters,                                                    &
161        ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,    &
162               initializing_actions, intermediate_timestep_count_max,          &
163               max_masks, precipitation, pt_surface,                           &
164               rho_surface, roughness_length, surface_pressure,                &
165               timestep_scheme, tsc, z0h_factor, time_since_reference_point
166
167    USE indices,                                                               &
168        ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner 
169
170    USE kinds
171
172    USE pegrid
173
174    USE radiation_model_mod,                                                   &
175        ONLY:  force_radiation_call, radiation_scheme, rad_net, rad_sw_in,     &
176               rad_lw_out, rad_lw_out_change_0, sigma_sb,                      &
177               unscheduled_radiation_calls
178       
179#if defined ( __rrtmg )
180    USE radiation_model_mod,                                                   &
181        ONLY:  rrtm_idrv
182#endif
183
184    USE statistics,                                                            &
185        ONLY:  hom, statistic_regions
186
187    IMPLICIT NONE
188
189!
190!-- LSM model constants
191    INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !< bottom of the soil model (to be switched)
192                               nzt_soil = 3, & !< top of the soil model (to be switched)
193                               nzs = 4         !< number of soil layers (fixed for now)
194
195    REAL(wp), PARAMETER ::                     &
196              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
197              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil   
198              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
199              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
200              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
201              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil
202              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water
203              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
204
205
206!
207!-- LSM variables
208    INTEGER(iwp) :: veg_type  = 2, & !< NAMELIST veg_type_2d
209                    soil_type = 3    !< NAMELIST soil_type_2d
210
211    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  soil_type_2d, &  !< soil type, 0: user-defined, 1-7: generic (see list)
212                                                  veg_type_2d      !< vegetation type, 0: user-defined, 1-19: generic (see list)
213
214    LOGICAL, DIMENSION(:,:), ALLOCATABLE :: water_surface,     & !< flag parameter for water surfaces (classes 14+15)
215                                            pave_surface,      & !< flag parameter for pavements (asphalt etc.) (class 20)
216                                            building_surface     !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
217
218    LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
219               force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls
220               land_surface = .FALSE.              !< flag parameter indicating wheather the lsm is used
221
222!   value 9999999.9_wp -> generic available or user-defined value must be set
223!   otherwise -> no generic variable and user setting is optional
224    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
225                canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
226                c_surface   = 20000.0_wp,               & !< Surface (skin) heat capacity
227                drho_l_lv,                              & !< (rho_l * l_v)**-1
228                exn,                                    & !< value of the Exner function
229                e_s = 0.0_wp,                           & !< saturation water vapour pressure
230                field_capacity = 9999999.9_wp,          & !< NAMELIST m_fc
231                f_shortwave_incoming = 9999999.9_wp,    & !< NAMELIST f_sw_in
232                hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
233                ke = 0.0_wp,                            & !< Kersten number
234                lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil
235                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s
236                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u
237                leaf_area_index = 9999999.9_wp,         & !< NAMELIST lai
238                l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
239                min_canopy_resistance = 9999999.9_wp,   & !< NAMELIST r_canopy_min
240                min_soil_resistance = 50.0_wp,          & !< NAMELIST r_soil_min
241                m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
242                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
243                pave_depth = 9999999.9_wp,              & !< depth of the pavement
244                pave_heat_capacity = 1.94E6_wp,         & !< volumetric heat capacity of pavement (e.g. roads)
245                pave_heat_conductivity = 1.00_wp,       & !< heat conductivity for pavements (e.g. roads)
246                q_s = 0.0_wp,                           & !< saturation specific humidity
247                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
248                rho_cp,                                 & !< rho_surface * cp
249                rho_lv,                                 & !< rho * l_v
250                rd_d_rv,                                & !< r_d / r_v
251                saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
252                skip_time_do_lsm = 0.0_wp,              & !< LSM is not called before this time
253                vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
254                wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
255                z0_eb  = 9999999.9_wp,                  & !< NAMELIST z0 (lsm_par)
256                z0h_eb = 9999999.9_wp,                  & !< NAMELIST z0h (lsm_par)
257                z0q_eb = 9999999.9_wp                     !< NAMELIST z0q (lsm_par)
258
259    REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
260              ddz_soil_stag,                  & !< 1/dz_soil_stag
261              dz_soil_stag,                   & !< soil grid spacing (center-center)
262              root_extr = 0.0_wp,             & !< root extraction
263              root_fraction = (/9999999.9_wp, 9999999.9_wp,    &
264                                9999999.9_wp, 9999999.9_wp /), & !< distribution of root surface area to the individual soil layers
265              zs = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),    & !< soil layer depths (m)
266              soil_moisture = 0.0_wp          !< soil moisture content (m3/m3)
267
268    REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) ::   &
269              soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp,  283.0_wp,    & !< soil temperature (K)
270                                   283.0_wp /),                                &                                   
271              ddz_soil,                                                        & !< 1/dz_soil
272              dz_soil                                                            !< soil grid spacing (edge-edge)
273
274#if defined( __nopointer )
275    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface,   & !< surface temperature (K)
276                                                     t_surface_p, & !< progn. surface temperature (K)
277                                                     m_liq_eb,    & !< liquid water reservoir (m)
278                                                     m_liq_eb_av, & !< liquid water reservoir (m)
279                                                     m_liq_eb_p     !< progn. liquid water reservoir (m)
280#else
281    REAL(wp), DIMENSION(:,:), POINTER :: t_surface,      &
282                                         t_surface_p,    & 
283                                         m_liq_eb,       & 
284                                         m_liq_eb_p
285
286    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, &
287                                                     m_liq_eb_av,              &
288                                                     m_liq_eb_1, m_liq_eb_2
289#endif
290
291!
292!-- Temporal tendencies for time stepping           
293    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m,  & !< surface temperature tendency (K)
294                                             tm_liq_eb_m      !< liquid water reservoir tendency (m)
295
296!
297!-- Energy balance variables               
298    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
299              alpha_vg,         & !< coef. of Van Genuchten
300              c_liq,            & !< liquid water coverage (of vegetated area)
301              c_liq_av,         & !< average of c_liq
302              c_soil_av,        & !< average of c_soil
303              c_veg,            & !< vegetation coverage
304              c_veg_av,         & !< average of c_veg
305              f_sw_in,          & !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
306              ghf_eb,           & !< ground heat flux
307              ghf_eb_av,        & !< average of ghf_eb
308              gamma_w_sat,      & !< hydraulic conductivity at saturation
309              g_d,              & !< coefficient for dependence of r_canopy on water vapour pressure deficit
310              lai,              & !< leaf area index
311              lai_av,           & !< average of lai
312              lambda_surface_s, & !< coupling between surface and soil (depends on vegetation type)
313              lambda_surface_u, & !< coupling between surface and soil (depends on vegetation type)
314              l_vg,             & !< coef. of Van Genuchten
315              m_fc,             & !< soil moisture at field capacity (m3/m3)
316              m_res,            & !< residual soil moisture
317              m_sat,            & !< saturation soil moisture (m3/m3)
318              m_wilt,           & !< soil moisture at permanent wilting point (m3/m3)
319              n_vg,             & !< coef. Van Genuchten 
320              qsws_eb,          & !< surface flux of latent heat (total)
321              qsws_eb_av,       & !< average of qsws_eb
322              qsws_liq_eb,      & !< surface flux of latent heat (liquid water portion)
323              qsws_liq_eb_av,   & !< average of qsws_liq_eb
324              qsws_soil_eb,     & !< surface flux of latent heat (soil portion)
325              qsws_soil_eb_av,  & !< average of qsws_soil_eb
326              qsws_veg_eb,      & !< surface flux of latent heat (vegetation portion)
327              qsws_veg_eb_av,   & !< average of qsws_veg_eb
328              rad_net_l,        & !< local copy of rad_net (net radiation at surface)
329              r_a,              & !< aerodynamic resistance
330              r_a_av,           & !< average of r_a
331              r_canopy,         & !< canopy resistance
332              r_soil,           & !< soil resistance
333              r_soil_min,       & !< minimum soil resistance
334              r_s,              & !< total surface resistance (combination of r_soil and r_canopy)
335              r_s_av,           & !< average of r_s
336              r_canopy_min,     & !< minimum canopy (stomatal) resistance
337              shf_eb,           & !< surface flux of sensible heat
338              shf_eb_av           !< average of shf_eb
339
340
341    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
342              lambda_h, &   !< heat conductivity of soil (W/m/K)                           
343              lambda_w, &   !< hydraulic diffusivity of soil (?)
344              gamma_w,  &   !< hydraulic conductivity of soil (W/m/K)
345              rho_c_total   !< volumetric heat capacity of the actual soil matrix (?)
346
347#if defined( __nopointer )
348    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
349              t_soil,    & !< Soil temperature (K)
350              t_soil_av, & !< Average of t_soil
351              t_soil_p,  & !< Prog. soil temperature (K)
352              m_soil,    & !< Soil moisture (m3/m3)
353              m_soil_av, & !< Average of m_soil
354              m_soil_p     !< Prog. soil moisture (m3/m3)
355#else
356    REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
357              t_soil, t_soil_p, &
358              m_soil, m_soil_p   
359
360    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
361              t_soil_av, t_soil_1, t_soil_2,                                   &
362              m_soil_av, m_soil_1, m_soil_2
363#endif
364
365
366    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
367              tt_soil_m, & !< t_soil storage array
368              tm_soil_m, & !< m_soil storage array
369              root_fr      !< root fraction (sum=1)
370
371
372!
373!-- Predefined Land surface classes (veg_type)
374    CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ &
375                                   'user defined              ',    & ! 0
376                                   'crops, mixed farming      ',    & !  1
377                                   'short grass               ',    & !  2
378                                   'evergreen needleleaf trees',    & !  3
379                                   'deciduous needleleaf trees',    & !  4
380                                   'evergreen broadleaf trees ',    & !  5
381                                   'deciduous broadleaf trees ',    & !  6
382                                   'tall grass                ',    & !  7
383                                   'desert                    ',    & !  8
384                                   'tundra                    ',    & !  9
385                                   'irrigated crops           ',    & ! 10
386                                   'semidesert                ',    & ! 11
387                                   'ice caps and glaciers     ',    & ! 12
388                                   'bogs and marshes          ',    & ! 13
389                                   'inland water              ',    & ! 14
390                                   'ocean                     ',    & ! 15
391                                   'evergreen shrubs          ',    & ! 16
392                                   'deciduous shrubs          ',    & ! 17
393                                   'mixed forest/woodland     ',    & ! 18
394                                   'interrupted forest        ',    & ! 19
395                                   'pavements/roads           '     & ! 20
396                                                                 /)
397
398!
399!-- Soil model classes (soil_type)
400    CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
401                                   'user defined',                  & ! 0
402                                   'coarse      ',                  & ! 1
403                                   'medium      ',                  & ! 2
404                                   'medium-fine ',                  & ! 3
405                                   'fine        ',                  & ! 4
406                                   'very fine   ',                  & ! 5
407                                   'organic     ',                  & ! 6
408                                   'loamy (CH)  '                   & ! 7
409                                                                 /)
410!
411!-- Land surface parameters according to the respective classes (veg_type)
412
413!
414!-- Land surface parameters I
415!--                          r_canopy_min,     lai,   c_veg,     g_d
416    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ &
417                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
418                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
419                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  3
420                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  4
421                                 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  5
422                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,  & !  6
423                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,  & !  7
424                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  & !  8
425                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  & !  9
426                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & ! 10
427                                 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,  & ! 11
428                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 12
429                                 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,  & ! 13
430                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 14
431                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 15
432                                 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,  & ! 16
433                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
434                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
435                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,  & ! 19
436                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp   & ! 20
437                                 /), (/ 4, 20 /) )
438
439!
440!-- Land surface parameters II          z0,         z0h,         z0q
441    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: roughness_par = RESHAPE( (/ & 
442                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & !  1
443                                   0.20_wp,  0.20E-2_wp,  0.20E-2_wp,       & !  2
444                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  3
445                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  4
446                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  5
447                                   2.00_wp,     2.00_wp,     2.00_wp,       & !  6
448                                   0.47_wp,  0.47E-2_wp,  0.47E-2_wp,       & !  7
449                                  0.013_wp, 0.013E-2_wp, 0.013E-2_wp,       & !  8
450                                  0.034_wp, 0.034E-2_wp, 0.034E-2_wp,       & !  9
451                                    0.5_wp,  0.50E-2_wp,  0.50E-2_wp,       & ! 10
452                                   0.17_wp,  0.17E-2_wp,  0.17E-2_wp,       & ! 11
453                                 1.3E-3_wp,   1.3E-4_wp,   1.3E-4_wp,       & ! 12
454                                   0.83_wp,  0.83E-2_wp,  0.83E-2_wp,       & ! 13
455                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 14
456                                   0.00_wp,     0.00_wp,     0.00_wp,       & ! 15
457                                   0.10_wp,  0.10E-2_wp,  0.10E-2_wp,       & ! 16
458                                   0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & ! 17
459                                   2.00_wp,  2.00E-2_wp,  2.00E-2_wp,       & ! 18
460                                   1.10_wp,  1.10E-2_wp,  1.10E-2_wp,       & ! 19
461                                 1.0E-4_wp,   1.0E-5_wp,   1.0E-5_wp        & ! 20
462                                 /), (/ 3, 20 /) )
463
464!
465!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
466    REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ &
467                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  1
468                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  2
469                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  3
470                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  4
471                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  5
472                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  6
473                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  7
474                                      15.0_wp,       15.0_wp, 0.00_wp,     & !  8
475                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  9
476                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 10
477                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 11
478                                      58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
479                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 13
480                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 14
481                                    1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 15
482                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
483                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
484                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
485                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 19
486                                       0.0_wp,        0.0_wp, 0.00_wp      & ! 20
487                                      /), (/ 3, 20 /) )
488
489!
490!-- Root distribution (sum = 1)  level 1, level 2, level 3, level 4,
491    REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: root_distribution = RESHAPE( (/ &
492                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  1
493                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  2
494                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  3
495                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  4
496                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  5
497                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  6
498                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  7
499                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  8
500                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & !  9
501                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 10
502                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 11
503                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 12
504                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 13
505                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 14
506                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 15
507                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16
508                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 17
509                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 18
510                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 19
511                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp             & ! 20
512                                 /), (/ 4, 20 /) )
513
514!
515!-- Soil parameters according to the following porosity classes (soil_type)
516
517!
518!-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
519    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/     &
520                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
521                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
522                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, & ! 3
523                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
524                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
525                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, & ! 6
526                                 0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  & ! 7
527                                 /), (/ 4, 7 /) )
528
529!
530!-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res 
531    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
532                                 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
533                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
534                                 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
535                                 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
536                                 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
537                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
538                                 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  & ! 7
539                                 /), (/ 4, 7 /) )
540
541
542    SAVE
543
544
545    PRIVATE
546
547   
548!
549!-- Public functions
550    PUBLIC lsm_check_data_output, lsm_check_data_output_pr,                    &
551           lsm_check_parameters, lsm_energy_balance, lsm_header, lsm_init,     &
552           lsm_init_arrays, lsm_parin, lsm_soil_model, lsm_swap_timelevel 
553!
554!-- Public parameters, constants and initial values
555    PUBLIC land_surface, skip_time_do_lsm
556
557!
558!-- Public grid variables
559    PUBLIC nzb_soil, nzs, nzt_soil, zs
560
561!
562!-- Public 2D output variables
563    PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av,     &
564           lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av,      &
565           qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av,         &
566           r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av
567
568!
569!-- Public prognostic variables
570    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
571
572
573    INTERFACE lsm_check_data_output
574       MODULE PROCEDURE lsm_check_data_output
575    END INTERFACE lsm_check_data_output
576   
577    INTERFACE lsm_check_data_output_pr
578       MODULE PROCEDURE lsm_check_data_output_pr
579    END INTERFACE lsm_check_data_output_pr
580   
581    INTERFACE lsm_check_parameters
582       MODULE PROCEDURE lsm_check_parameters
583    END INTERFACE lsm_check_parameters
584   
585    INTERFACE lsm_energy_balance
586       MODULE PROCEDURE lsm_energy_balance
587    END INTERFACE lsm_energy_balance
588
589    INTERFACE lsm_header
590       MODULE PROCEDURE lsm_header
591    END INTERFACE lsm_header
592   
593    INTERFACE lsm_init
594       MODULE PROCEDURE lsm_init
595    END INTERFACE lsm_init
596
597    INTERFACE lsm_init_arrays
598       MODULE PROCEDURE lsm_init_arrays
599    END INTERFACE lsm_init_arrays
600   
601    INTERFACE lsm_parin
602       MODULE PROCEDURE lsm_parin
603    END INTERFACE lsm_parin
604   
605    INTERFACE lsm_soil_model
606       MODULE PROCEDURE lsm_soil_model
607    END INTERFACE lsm_soil_model
608
609    INTERFACE lsm_swap_timelevel
610       MODULE PROCEDURE lsm_swap_timelevel
611    END INTERFACE lsm_swap_timelevel
612
613 CONTAINS
614
615!------------------------------------------------------------------------------!
616! Description:
617! ------------
618!> Check data output for land surface model
619!------------------------------------------------------------------------------!
620    SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
621 
622 
623       USE control_parameters,                                                 &
624           ONLY: data_output, message_string
625
626       IMPLICIT NONE
627
628       CHARACTER (LEN=*) ::  unit     !<
629       CHARACTER (LEN=*) ::  var      !<
630
631       INTEGER(iwp) :: i
632       INTEGER(iwp) :: ilen   
633       INTEGER(iwp) :: k
634
635       SELECT CASE ( TRIM( var ) )
636
637          CASE ( 'm_soil' )
638             IF (  .NOT.  land_surface )  THEN
639                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
640                         'res land_surface = .TRUE.'
641                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
642             ENDIF
643             unit = 'm3/m3'
644             
645          CASE ( 't_soil' )
646             IF (  .NOT.  land_surface )  THEN
647                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
648                         'res land_surface = .TRUE.'
649                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
650             ENDIF
651             unit = 'K'   
652             
653          CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'm_liq_eb*',&
654                 'qsws_eb*', 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*',  &
655                 'r_a*', 'r_s*', 'shf_eb*' )
656             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
657                message_string = 'illegal value for data_output: "' //         &
658                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
659                                 'cross sections are allowed for this value'
660                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
661             ENDIF
662             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
663                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
664                                 'res land_surface = .TRUE.'
665                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
666             ENDIF
667             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
668                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
669                                 'res land_surface = .TRUE.'
670                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
671             ENDIF
672             IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
673                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
674                                 'res land_surface = .TRUE.'
675                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
676             ENDIF
677             IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
678                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
679                                 'res land_surface = .TRUE.'
680                CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
681             ENDIF
682             IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
683                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
684                                 'res land_surface = .TRUE.'
685                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
686             ENDIF
687             IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
688                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
689                                 'res land_surface = .TRUE.'
690                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
691             ENDIF
692             IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT.  land_surface )  THEN
693                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
694                                 'res land_surface = .TRUE.'
695                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
696             ENDIF
697             IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )   &
698             THEN
699                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
700                                 'res land_surface = .TRUE.'
701                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
702             ENDIF
703             IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
704             THEN
705                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
706                                 'res land_surface = .TRUE.'
707                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
708             ENDIF
709             IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
710             THEN
711                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
712                                 'res land_surface = .TRUE.'
713                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
714             ENDIF
715             IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface ) &
716             THEN
717                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
718                                 'res land_surface = .TRUE.'
719                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
720             ENDIF
721             IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface ) &
722             THEN
723                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
724                                 'res land_surface = .TRUE.'
725                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
726             ENDIF
727
728             IF ( TRIM( var ) == 'lai*'   )  unit = 'none' 
729             IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
730             IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
731             IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
732             IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
733             IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
734             IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
735             IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
736             IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
737             IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
738             IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
739             IF ( TRIM( var ) == 'r_s*')     unit = 's/m' 
740             IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
741             
742          CASE DEFAULT
743             unit = 'illegal'
744
745       END SELECT
746
747
748    END SUBROUTINE lsm_check_data_output
749
750 SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
751 
752       USE control_parameters,                                                 &
753           ONLY: data_output_pr, message_string
754
755       USE indices
756
757       USE profil_parameter
758
759       USE statistics
760
761       IMPLICIT NONE
762   
763       CHARACTER (LEN=*) ::  unit      !<
764       CHARACTER (LEN=*) ::  variable  !<
765       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
766 
767       INTEGER(iwp) ::  user_pr_index !<
768       INTEGER(iwp) ::  var_count     !<
769
770       SELECT CASE ( TRIM( variable ) )
771       
772          CASE ( 't_soil', '#t_soil' )
773             IF (  .NOT.  land_surface )  THEN
774                message_string = 'data_output_pr = ' //                        &
775                                 TRIM( data_output_pr(var_count) ) // ' is' // &
776                                 'not implemented for land_surface = .FALSE.'
777                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
778             ELSE
779                dopr_index(var_count) = 89
780                dopr_unit     = 'K'
781                hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
782                IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
783                   dopr_initial_index(var_count) = 90
784                   hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
785                   data_output_pr(var_count)     = data_output_pr(var_count)(2:)
786                ENDIF
787                unit = dopr_unit
788             ENDIF
789
790          CASE ( 'm_soil', '#m_soil' )
791             IF (  .NOT.  land_surface )  THEN
792                message_string = 'data_output_pr = ' //                        &
793                                 TRIM( data_output_pr(var_count) ) // ' is' // &
794                                 ' not implemented for land_surface = .FALSE.'
795                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
796             ELSE
797                dopr_index(var_count) = 91
798                dopr_unit     = 'm3/m3'
799                hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
800                IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
801                   dopr_initial_index(var_count) = 92
802                   hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
803                   data_output_pr(var_count)     = data_output_pr(var_count)(2:)
804                ENDIF
805                 unit = dopr_unit
806             ENDIF
807
808
809          CASE DEFAULT
810             unit = 'illegal'
811
812       END SELECT
813
814
815    END SUBROUTINE lsm_check_data_output_pr
816 
817 
818!------------------------------------------------------------------------------!
819! Description:
820! ------------
821!> Check parameters routine for land surface model
822!------------------------------------------------------------------------------!
823    SUBROUTINE lsm_check_parameters
824
825       USE control_parameters,                                                 &
826           ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string,         &
827                 most_method, topography
828                 
829       USE radiation_model_mod,                                                &
830           ONLY: radiation
831   
832   
833       IMPLICIT NONE
834
835 
836!
837!--    Dirichlet boundary conditions are required as the surface fluxes are
838!--    calculated from the temperature/humidity gradients in the land surface
839!--    model
840       IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
841          message_string = 'lsm requires setting of'//                         &
842                           'bc_pt_b = "dirichlet" and '//                      &
843                           'bc_q_b  = "dirichlet"'
844          CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
845       ENDIF
846
847       IF (  .NOT.  constant_flux_layer )  THEN
848          message_string = 'lsm requires '//                                   &
849                           'constant_flux_layer = .T.'
850          CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
851       ENDIF
852
853       IF ( topography /= 'flat' )  THEN
854          message_string = 'lsm cannot be used ' //                            & 
855                           'in combination with  topography /= "flat"'
856          CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
857       ENDIF
858
859       IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                     &
860              most_method == 'lookup' )  THEN
861           WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ',     &
862                                      'allowed in combination with ',          &
863                                      'most_method = ', most_method
864          CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
865       ENDIF
866
867       IF ( veg_type == 0 )  THEN
868          IF ( SUM( root_fraction ) /= 1.0_wp )  THEN
869             message_string = 'veg_type = 0 (user_defined)'//                  &
870                              'requires setting of root_fraction(0:3)'//       &
871                              '/= 9999999.9 and SUM(root_fraction) = 1'
872             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
873          ENDIF
874 
875          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
876             message_string = 'veg_type = 0 (user defined)'//                  &
877                              'requires setting of min_canopy_resistance'//    &
878                              '/= 9999999.9'
879             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
880          ENDIF
881
882          IF ( leaf_area_index == 9999999.9_wp )  THEN
883             message_string = 'veg_type = 0 (user_defined)'//                  &
884                              'requires setting of leaf_area_index'//          &
885                              '/= 9999999.9'
886             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
887          ENDIF
888
889          IF ( vegetation_coverage == 9999999.9_wp )  THEN
890             message_string = 'veg_type = 0 (user_defined)'//                  &
891                              'requires setting of vegetation_coverage'//      &
892                              '/= 9999999.9'
893             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
894          ENDIF
895
896          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
897             message_string = 'veg_type = 0 (user_defined)'//                  &
898                              'requires setting of'//                          &
899                              'canopy_resistance_coefficient /= 9999999.9'
900             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
901          ENDIF
902
903          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
904             message_string = 'veg_type = 0 (user_defined)'//                  &
905                              'requires setting of lambda_surface_stable'//    &
906                              '/= 9999999.9'
907             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
908          ENDIF
909
910          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
911             message_string = 'veg_type = 0 (user_defined)'//                  &
912                              'requires setting of lambda_surface_unstable'//  &
913                              '/= 9999999.9'
914             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
915          ENDIF
916
917          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
918             message_string = 'veg_type = 0 (user_defined)'//                  &
919                              'requires setting of f_shortwave_incoming'//     &
920                              '/= 9999999.9'
921             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
922          ENDIF
923
924          IF ( z0_eb == 9999999.9_wp )  THEN
925             message_string = 'veg_type = 0 (user_defined)'//                  &
926                              'requires setting of z0_eb'//                    &
927                              '/= 9999999.9'
928             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
929          ENDIF
930
931          IF ( z0h_eb == 9999999.9_wp )  THEN
932             message_string = 'veg_type = 0 (user_defined)'//                  &
933                              'requires setting of z0h_eb'//                   &
934                              '/= 9999999.9'
935             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
936          ENDIF
937
938
939       ENDIF
940
941       IF ( soil_type == 0 )  THEN
942
943          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
944             message_string = 'soil_type = 0 (user_defined)'//                 &
945                              'requires setting of alpha_vangenuchten'//       &
946                              '/= 9999999.9'
947             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
948          ENDIF
949
950          IF ( l_vangenuchten == 9999999.9_wp )  THEN
951             message_string = 'soil_type = 0 (user_defined)'//                 &
952                              'requires setting of l_vangenuchten'//           &
953                              '/= 9999999.9'
954             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
955          ENDIF
956
957          IF ( n_vangenuchten == 9999999.9_wp )  THEN
958             message_string = 'soil_type = 0 (user_defined)'//                 &
959                              'requires setting of n_vangenuchten'//           &
960                              '/= 9999999.9'
961             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
962          ENDIF
963
964          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
965             message_string = 'soil_type = 0 (user_defined)'//                 &
966                              'requires setting of hydraulic_conductivity'//   &
967                              '/= 9999999.9'
968             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
969          ENDIF
970
971          IF ( saturation_moisture == 9999999.9_wp )  THEN
972             message_string = 'soil_type = 0 (user_defined)'//                 &
973                              'requires setting of saturation_moisture'//      &
974                              '/= 9999999.9'
975             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
976          ENDIF
977
978          IF ( field_capacity == 9999999.9_wp )  THEN
979             message_string = 'soil_type = 0 (user_defined)'//                 &
980                              'requires setting of field_capacity'//           &
981                              '/= 9999999.9'
982             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
983          ENDIF
984
985          IF ( wilting_point == 9999999.9_wp )  THEN
986             message_string = 'soil_type = 0 (user_defined)'//                 &
987                              'requires setting of wilting_point'//            &
988                              '/= 9999999.9'
989             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
990          ENDIF
991
992          IF ( residual_moisture == 9999999.9_wp )  THEN
993             message_string = 'soil_type = 0 (user_defined)'//                 &
994                              'requires setting of residual_moisture'//        &
995                              '/= 9999999.9'
996             CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
997          ENDIF
998
999       ENDIF
1000
1001       IF (  .NOT.  radiation )  THEN
1002          message_string = 'lsm requires '//                                   &
1003                           'radiation = .T.'
1004          CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
1005       ENDIF
1006       
1007       
1008    END SUBROUTINE lsm_check_parameters
1009 
1010!------------------------------------------------------------------------------!
1011! Description:
1012! ------------
1013!> Solver for the energy balance at the surface.
1014!------------------------------------------------------------------------------!
1015    SUBROUTINE lsm_energy_balance
1016
1017
1018       IMPLICIT NONE
1019
1020       INTEGER(iwp) ::  i         !< running index
1021       INTEGER(iwp) ::  j         !< running index
1022       INTEGER(iwp) ::  k, ks     !< running index
1023
1024       REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
1025                   f1,          & !< resistance correction term 1
1026                   f2,          & !< resistance correction term 2
1027                   f3,          & !< resistance correction term 3
1028                   m_min,       & !< minimum soil moisture
1029                   e,           & !< water vapour pressure
1030                   e_s,         & !< water vapour saturation pressure
1031                   e_s_dt,      & !< derivate of e_s with respect to T
1032                   tend,        & !< tendency
1033                   dq_s_dt,     & !< derivate of q_s with respect to T
1034                   coef_1,      & !< coef. for prognostic equation
1035                   coef_2,      & !< coef. for prognostic equation
1036                   f_qsws,      & !< factor for qsws_eb
1037                   f_qsws_veg,  & !< factor for qsws_veg_eb
1038                   f_qsws_soil, & !< factor for qsws_soil_eb
1039                   f_qsws_liq,  & !< factor for qsws_liq_eb
1040                   f_shf,       & !< factor for shf_eb
1041                   lambda_surface, & !< Current value of lambda_surface
1042                   m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
1043                   pt1,         & !< potential temperature at first grid level
1044                   qv1            !< specific humidity at first grid level
1045
1046!
1047!--    Calculate the exner function for the current time step
1048       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
1049
1050       DO  i = nxlg, nxrg
1051          DO  j = nysg, nyng
1052             k = nzb_s_inner(j,i)
1053
1054!
1055!--          Set lambda_surface according to stratification between skin layer and soil
1056             IF (  .NOT.  pave_surface(j,i) )  THEN
1057
1058                c_surface_tmp = c_surface
1059
1060                IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
1061                   lambda_surface = lambda_surface_s(j,i)
1062                ELSE
1063                   lambda_surface = lambda_surface_u(j,i)
1064                ENDIF
1065             ELSE
1066
1067                c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
1068                lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
1069
1070             ENDIF
1071
1072!
1073!--          First step: calculate aerodyamic resistance. As pt, us, ts
1074!--          are not available for the prognostic time step, data from the last
1075!--          time step is used here. Note that this formulation is the
1076!--          equivalent to the ECMWF formulation using drag coefficients
1077             IF ( cloud_physics )  THEN
1078                pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
1079                qv1 = q(k+1,j,i) - ql(k+1,j,i)
1080             ELSE
1081                pt1 = pt(k+1,j,i)
1082                qv1 = q(k+1,j,i)
1083             ENDIF
1084
1085             r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
1086
1087!
1088!--          Make sure that the resistance does not drop to zero for neutral
1089!--          stratification
1090             IF ( ABS(r_a(j,i)) < 1.0_wp )  r_a(j,i) = 1.0_wp
1091
1092!
1093!--          Second step: calculate canopy resistance r_canopy
1094!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
1095 
1096!--          f1: correction for incoming shortwave radiation (stomata close at
1097!--          night)
1098             IF ( radiation_scheme /= 'constant' )  THEN
1099                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /  &
1100                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)          &
1101                               + 1.0_wp)) )
1102             ELSE
1103                f1 = 1.0_wp
1104             ENDIF
1105
1106
1107!
1108!--          f2: correction for soil moisture availability to plants (the
1109!--          integrated soil moisture must thus be considered here)
1110!--          f2 = 0 for very dry soils
1111             m_total = 0.0_wp
1112             DO  ks = nzb_soil, nzt_soil
1113                 m_total = m_total + root_fr(ks,j,i)                           &
1114                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
1115             ENDDO 
1116
1117             IF ( m_total > m_wilt(j,i)  .AND.  m_total < m_fc(j,i) )  THEN
1118                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
1119             ELSEIF ( m_total >= m_fc(j,i) )  THEN
1120                f2 = 1.0_wp
1121             ELSE
1122                f2 = 1.0E-20_wp
1123             ENDIF
1124
1125!
1126!--          Calculate water vapour pressure at saturation
1127             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
1128                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
1129
1130!
1131!--          f3: correction for vapour pressure deficit
1132             IF ( g_d(j,i) /= 0.0_wp )  THEN
1133!
1134!--             Calculate vapour pressure
1135                e  = qv1 * surface_pressure / 0.622_wp
1136                f3 = EXP ( -g_d(j,i) * (e_s - e) )
1137             ELSE
1138                f3 = 1.0_wp
1139             ENDIF
1140
1141!
1142!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1143!--          this calculation is obsolete, as r_canopy is not used below.
1144!--          To do: check for very dry soil -> r_canopy goes to infinity
1145             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
1146                                             + 1.0E-20_wp)
1147
1148!
1149!--          Third step: calculate bare soil resistance r_soil. The Clapp &
1150!--          Hornberger parametrization does not consider c_veg.
1151             IF ( soil_type_2d(j,i) /= 7 )  THEN
1152                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
1153                        m_res(j,i)
1154             ELSE
1155                m_min = m_wilt(j,i)
1156             ENDIF
1157
1158             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
1159             f2 = MAX(f2,1.0E-20_wp)
1160             f2 = MIN(f2,1.0_wp)
1161
1162             r_soil(j,i) = r_soil_min(j,i) / f2
1163
1164!
1165!--          Calculate the maximum possible liquid water amount on plants and
1166!--          bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
1167!--          assumed, while paved surfaces might hold up 1 mm of water. The
1168!--          liquid water fraction for paved surfaces is calculated after
1169!--          Noilhan & Planton (1989), while the ECMWF formulation is used for
1170!--          vegetated surfaces and bare soils.
1171             IF ( pave_surface(j,i) )  THEN
1172                m_liq_eb_max = m_max_depth * 5.0_wp
1173                c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
1174             ELSE
1175                m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)           &
1176                               + (1.0_wp - c_veg(j,i)) )
1177                c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
1178             ENDIF
1179
1180!
1181!--          Calculate saturation specific humidity
1182             q_s = 0.622_wp * e_s / surface_pressure
1183
1184!
1185!--          In case of dewfall, set evapotranspiration to zero
1186!--          All super-saturated water is then removed from the air
1187             IF ( humidity  .AND.  q_s <= qv1 )  THEN
1188                r_canopy(j,i) = 0.0_wp
1189                r_soil(j,i)   = 0.0_wp
1190             ENDIF
1191
1192!
1193!--          Calculate coefficients for the total evapotranspiration
1194!--          In case of water surface, set vegetation and soil fluxes to zero.
1195!--          For pavements, only evaporation of liquid water is possible.
1196             IF ( water_surface(j,i) )  THEN
1197                f_qsws_veg  = 0.0_wp
1198                f_qsws_soil = 0.0_wp
1199                f_qsws_liq  = rho_lv / r_a(j,i)
1200             ELSEIF ( pave_surface (j,i) )  THEN
1201                f_qsws_veg  = 0.0_wp
1202                f_qsws_soil = 0.0_wp
1203                f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
1204             ELSE
1205                f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/     &
1206                              (r_a(j,i) + r_canopy(j,i))
1207                f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +     &
1208                                                          r_soil(j,i))
1209                f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
1210             ENDIF
1211!
1212!--          If soil moisture is below wilting point, plants do no longer
1213!--          transpirate.
1214!              IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
1215!                 f_qsws_veg = 0.0_wp
1216!              ENDIF
1217
1218             f_shf  = rho_cp / r_a(j,i)
1219             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1220
1221!
1222!--          Calculate derivative of q_s for Taylor series expansion
1223             e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
1224                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
1225                              / (t_surface(j,i) - 35.86_wp)**2 )
1226
1227             dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
1228
1229!
1230!--          Add LW up so that it can be removed in prognostic equation
1231             rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
1232
1233!
1234!--          Calculate new skin temperature
1235             IF ( humidity )  THEN
1236#if defined ( __rrtmg )
1237!
1238!--             Numerator of the prognostic equation
1239                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
1240                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
1241                         + f_shf * pt1 + f_qsws * ( qv1 - q_s                  &
1242                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
1243                         * t_soil(nzb_soil,j,i)
1244
1245!
1246!--             Denominator of the prognostic equation
1247                coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt           &
1248                         + lambda_surface + f_shf / exn
1249#else
1250
1251!
1252!--             Numerator of the prognostic equation
1253                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
1254                         * t_surface(j,i) ** 4                                 &
1255                         + f_shf * pt1 + f_qsws * ( qv1                        &
1256                         - q_s + dq_s_dt * t_surface(j,i) )                    &
1257                         + lambda_surface * t_soil(nzb_soil,j,i)
1258
1259!
1260!--             Denominator of the prognostic equation
1261                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
1262                         * dq_s_dt + lambda_surface + f_shf / exn
1263 
1264#endif
1265             ELSE
1266
1267#if defined ( __rrtmg )
1268!
1269!--             Numerator of the prognostic equation
1270                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
1271                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
1272                         + f_shf * pt1  + lambda_surface                       &
1273                         * t_soil(nzb_soil,j,i)
1274
1275!
1276!--             Denominator of the prognostic equation
1277                coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
1278#else
1279
1280!
1281!--             Numerator of the prognostic equation
1282                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
1283                          * t_surface(j,i) ** 4 + f_shf * pt1                  &
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               &
1289                         + lambda_surface + f_shf / exn
1290#endif
1291             ENDIF
1292
1293             tend = 0.0_wp
1294
1295!
1296!--          Implicit solution when the surface layer has no heat capacity,
1297!--          otherwise use RK3 scheme.
1298             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *    &
1299                                t_surface(j,i) ) / ( c_surface_tmp + coef_2    &
1300                                * dt_3d * tsc(2) ) 
1301
1302!
1303!--          Add RK3 term
1304             IF ( c_surface_tmp /= 0.0_wp )  THEN
1305
1306                t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
1307                                   * tt_surface_m(j,i)
1308
1309!
1310!--             Calculate true tendency
1311                tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)     &
1312                       * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
1313!
1314!--             Calculate t_surface tendencies for the next Runge-Kutta step
1315                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1316                   IF ( intermediate_timestep_count == 1 )  THEN
1317                      tt_surface_m(j,i) = tend
1318                   ELSEIF ( intermediate_timestep_count <                      &
1319                            intermediate_timestep_count_max )  THEN
1320                      tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp        &
1321                                          * tt_surface_m(j,i)
1322                   ENDIF
1323                ENDIF
1324             ENDIF
1325
1326!
1327!--          In case of fast changes in the skin temperature, it is possible to
1328!--          update the radiative fluxes independently from the prescribed
1329!--          radiation call frequency. This effectively prevents oscillations,
1330!--          especially when setting skip_time_do_radiation /= 0. The threshold
1331!--          value of 0.2 used here is just a first guess. This method should be
1332!--          revised in the future as tests have shown that the threshold is
1333!--          often reached, when no oscillations would occur (causes immense
1334!--          computing time for the radiation code).
1335             IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.     &
1336                  unscheduled_radiation_calls )  THEN
1337                force_radiation_call_l = .TRUE.
1338             ENDIF
1339
1340             pt(k,j,i) = t_surface_p(j,i) / exn
1341
1342!
1343!--          Calculate fluxes
1344#if defined ( __rrtmg )
1345             rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)      &
1346                                * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
1347                                - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
1348
1349             IF ( rrtm_idrv == 1 )  THEN
1350                rad_net(j,i) = rad_net_l(j,i)
1351                rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
1352                                      + rad_lw_out_change_0(j,i)               &
1353                                      * ( t_surface_p(j,i) - t_surface(j,i) )
1354             ENDIF
1355#else
1356             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb             &
1357                                * t_surface(j,i)**4 - 4.0_wp * sigma_sb        &
1358                                * t_surface(j,i)**3 * t_surface_p(j,i)
1359#endif
1360
1361             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
1362                              - t_soil(nzb_soil,j,i))
1363
1364             shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
1365
1366             shf(j,i) = shf_eb(j,i) / rho_cp
1367
1368             IF ( humidity )  THEN
1369                qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt            &
1370                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
1371
1372                qsws(j,i) = qsws_eb(j,i) / rho_lv
1373
1374                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                &
1375                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1376                                    * t_surface_p(j,i) )
1377
1378                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                &
1379                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1380                                    * t_surface_p(j,i) )
1381
1382                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                &
1383                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1384                                    * t_surface_p(j,i) )
1385             ENDIF
1386
1387!
1388!--          Calculate the true surface resistance
1389             IF ( qsws_eb(j,i) == 0.0_wp )  THEN
1390                r_s(j,i) = 1.0E10_wp
1391             ELSE
1392                r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                    &
1393                           * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
1394                           / qsws_eb(j,i) - r_a(j,i)
1395             ENDIF
1396
1397!
1398!--          Calculate change in liquid water reservoir due to dew fall or
1399!--          evaporation of liquid water
1400             IF ( humidity )  THEN
1401!
1402!--             If precipitation is activated, add rain water to qsws_liq_eb
1403!--             and qsws_soil_eb according the the vegetation coverage.
1404!--             precipitation_rate is given in mm.
1405                IF ( precipitation )  THEN
1406
1407!
1408!--                Add precipitation to liquid water reservoir, if possible.
1409!--                Otherwise, add the water to soil. In case of
1410!--                pavements, the exceeding water amount is implicitely removed
1411!--                as runoff as qsws_soil_eb is then not used in the soil model
1412                   IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
1413                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
1414                                       + c_veg(j,i) * prr(k,j,i) * hyrho(k)    &
1415                                       * 0.001_wp * rho_l * l_v
1416                   ELSE
1417                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1418                                        + c_veg(j,i) * prr(k,j,i) * hyrho(k)   &
1419                                        * 0.001_wp * rho_l * l_v
1420                   ENDIF
1421
1422!--                Add precipitation to bare soil according to the bare soil
1423!--                coverage.
1424                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i) + (1.0_wp             &
1425                                       - c_veg(j,i)) * prr(k,j,i) * hyrho(k)   &
1426                                       * 0.001_wp * rho_l * l_v
1427                ENDIF
1428
1429!
1430!--             If the air is saturated, check the reservoir water level
1431                IF ( qsws_eb(j,i) < 0.0_wp )  THEN
1432
1433!
1434!--                Check if reservoir is full (avoid values > m_liq_eb_max)
1435!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this
1436!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
1437!--                so that tend is zero and no further check is needed
1438                   IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
1439                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1440                                           + qsws_liq_eb(j,i)
1441
1442                      qsws_liq_eb(j,i)  = 0.0_wp
1443                   ENDIF
1444
1445!
1446!--                In case qsws_veg_eb becomes negative (unphysical behavior),
1447!--                let the water enter the liquid water reservoir as dew on the
1448!--                plant
1449                   IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
1450                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
1451                      qsws_veg_eb(j,i) = 0.0_wp
1452                   ENDIF
1453                ENDIF                   
1454 
1455                tend = - qsws_liq_eb(j,i) * drho_l_lv
1456
1457                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
1458                                                   + tsc(3) * tm_liq_eb_m(j,i) )
1459
1460!
1461!--             Check if reservoir is overfull -> reduce to maximum
1462!--             (conservation of water is violated here)
1463                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
1464
1465!
1466!--             Check if reservoir is empty (avoid values < 0.0)
1467!--             (conservation of water is violated here)
1468                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
1469
1470
1471!
1472!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
1473                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1474                   IF ( intermediate_timestep_count == 1 )  THEN
1475                      tm_liq_eb_m(j,i) = tend
1476                   ELSEIF ( intermediate_timestep_count <                      &
1477                            intermediate_timestep_count_max )  THEN
1478                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
1479                                      * tm_liq_eb_m(j,i)
1480                   ENDIF
1481                ENDIF
1482
1483             ENDIF
1484
1485          ENDDO
1486       ENDDO
1487
1488!
1489!--    Make a logical OR for all processes. Force radiation call if at
1490!--    least one processor reached the threshold change in skin temperature
1491       IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count    &
1492            == intermediate_timestep_count_max-1 )  THEN
1493#if defined( __parallel )
1494          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1495          CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
1496                              1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
1497#else
1498          force_radiation_call = force_radiation_call_l
1499#endif
1500          force_radiation_call_l = .FALSE.
1501       ENDIF
1502
1503!
1504!--    Calculate surface specific humidity
1505       IF ( humidity )  THEN
1506          CALL calc_q_surface
1507       ENDIF
1508
1509!
1510!--    Calculate new roughness lengths (for water surfaces only)
1511       CALL calc_z0_water_surface
1512
1513
1514    END SUBROUTINE lsm_energy_balance
1515
1516!------------------------------------------------------------------------------!
1517! Description:
1518! ------------
1519!> Header output for land surface model
1520!------------------------------------------------------------------------------!
1521    SUBROUTINE lsm_header ( io )
1522
1523
1524       IMPLICIT NONE
1525
1526       CHARACTER (LEN=86) ::  t_soil_chr          !< String for soil temperature profile
1527       CHARACTER (LEN=86) ::  roots_chr           !< String for root profile
1528       CHARACTER (LEN=86) ::  vertical_index_chr  !< String for the vertical index
1529       CHARACTER (LEN=86) ::  m_soil_chr          !< String for soil moisture
1530       CHARACTER (LEN=86) ::  soil_depth_chr      !< String for soil depth
1531       CHARACTER (LEN=10) ::  coor_chr            !< Temporary string
1532   
1533       INTEGER(iwp) ::  i                         !< Loop index over soil layers
1534 
1535       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1536 
1537       t_soil_chr = ''
1538       m_soil_chr    = ''
1539       soil_depth_chr  = '' 
1540       roots_chr        = '' 
1541       vertical_index_chr   = ''
1542
1543       i = 1
1544       DO i = nzb_soil, nzt_soil
1545          WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i)
1546          t_soil_chr = TRIM( t_soil_chr ) // ' ' // TRIM( coor_chr )
1547
1548          WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i)
1549          m_soil_chr = TRIM( m_soil_chr ) // ' ' // TRIM( coor_chr )
1550
1551          WRITE (coor_chr,'(F10.2,7X)')  - zs(i)
1552          soil_depth_chr = TRIM( soil_depth_chr ) // ' '  // TRIM( coor_chr )
1553
1554          WRITE (coor_chr,'(F10.2,7X)')  root_fraction(i)
1555          roots_chr = TRIM( roots_chr ) // ' '  // TRIM( coor_chr )
1556
1557          WRITE (coor_chr,'(I10,7X)')  i
1558          vertical_index_chr = TRIM( vertical_index_chr ) // ' '  //           &
1559                               TRIM( coor_chr )
1560       ENDDO
1561
1562!
1563!--    Write land surface model header
1564       WRITE( io,  1 )
1565       IF ( conserve_water_content )  THEN
1566          WRITE( io, 2 )
1567       ELSE
1568          WRITE( io, 3 )
1569       ENDIF
1570
1571       WRITE( io, 4 ) TRIM( veg_type_name(veg_type) ),                         &
1572                        TRIM (soil_type_name(soil_type) )
1573       WRITE( io, 5 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
1574                        TRIM( m_soil_chr ), TRIM( roots_chr ),                 &
1575                        TRIM( vertical_index_chr )
1576
15771   FORMAT (//' Land surface model information:'/                              &
1578              ' ------------------------------'/)
15792   FORMAT ('    --> Soil bottom is closed (water content is conserved',       &
1580            ', default)')
15813   FORMAT ('    --> Soil bottom is open (water content is not conserved)')         
15824   FORMAT ('    --> Land surface type  : ',A,/                                &
1583            '    --> Soil porosity type : ',A)
15845   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
1585            '       Height:        ',A,'  m'/                                  &
1586            '       Temperature:   ',A,'  K'/                                  &
1587            '       Moisture:      ',A,'  m**3/m**3'/                          &
1588            '       Root fraction: ',A,'  '/                                   &
1589            '       Grid point:    ',A)
1590
1591    END SUBROUTINE lsm_header
1592
1593
1594!------------------------------------------------------------------------------!
1595! Description:
1596! ------------
1597!> Initialization of the land surface model
1598!------------------------------------------------------------------------------!
1599    SUBROUTINE lsm_init
1600   
1601
1602       IMPLICIT NONE
1603
1604       INTEGER(iwp) ::  i !< running index
1605       INTEGER(iwp) ::  j !< running index
1606       INTEGER(iwp) ::  k !< running index
1607
1608       REAL(wp) :: pt1   !< potential temperature at first grid level
1609
1610
1611!
1612!--    Calculate Exner function
1613       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
1614
1615
1616!
1617!--    If no cloud physics is used, rho_surface has not been calculated before
1618       IF (  .NOT.  cloud_physics )  THEN
1619          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
1620       ENDIF
1621
1622!
1623!--    Calculate frequently used parameters
1624       rho_cp    = cp * rho_surface
1625       rd_d_rv   = r_d / r_v
1626       rho_lv    = rho_surface * l_v
1627       drho_l_lv = 1.0_wp / (rho_l * l_v)
1628
1629!
1630!--    Set inital values for prognostic quantities
1631       tt_surface_m = 0.0_wp
1632       tt_soil_m    = 0.0_wp
1633       tm_soil_m    = 0.0_wp
1634       tm_liq_eb_m  = 0.0_wp
1635       c_liq        = 0.0_wp
1636
1637       ghf_eb = 0.0_wp
1638       shf_eb = rho_cp * shf
1639
1640       IF ( humidity )  THEN
1641          qsws_eb = rho_lv * qsws
1642       ELSE
1643          qsws_eb = 0.0_wp
1644       ENDIF
1645
1646       qsws_liq_eb  = 0.0_wp
1647       qsws_soil_eb = 0.0_wp
1648       qsws_veg_eb  = 0.0_wp
1649
1650       r_a        = 50.0_wp
1651       r_s        = 50.0_wp
1652       r_canopy   = 0.0_wp
1653       r_soil     = 0.0_wp
1654
1655!
1656!--    Allocate 3D soil model arrays
1657       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1658       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1659       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1660
1661       lambda_h = 0.0_wp
1662!
1663!--    If required, allocate humidity-related variables for the soil model
1664       IF ( humidity )  THEN
1665          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1666          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
1667
1668          lambda_w = 0.0_wp 
1669       ENDIF
1670
1671!
1672!--    Calculate grid spacings. Temperature and moisture are defined at
1673!--    the edges of the soil layers (_stag), whereas gradients/fluxes are defined
1674!--    at the centers
1675       dz_soil(nzb_soil) = zs(nzb_soil)
1676
1677       DO  k = nzb_soil+1, nzt_soil
1678          dz_soil(k) = zs(k) - zs(k-1)
1679       ENDDO
1680       dz_soil(nzt_soil+1) = dz_soil(nzt_soil)
1681
1682       DO  k = nzb_soil, nzt_soil-1
1683          dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k))
1684       ENDDO
1685       dz_soil_stag(nzt_soil) = dz_soil(nzt_soil)
1686
1687       ddz_soil      = 1.0_wp / dz_soil
1688       ddz_soil_stag = 1.0_wp / dz_soil_stag
1689
1690!
1691!--    Initialize standard soil types. It is possible to overwrite each
1692!--    parameter by setting the respecticy NAMELIST variable to a
1693!--    value /= 9999999.9.
1694       IF ( soil_type /= 0 )  THEN 
1695 
1696          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
1697             alpha_vangenuchten = soil_pars(0,soil_type)
1698          ENDIF
1699
1700          IF ( l_vangenuchten == 9999999.9_wp )  THEN
1701             l_vangenuchten = soil_pars(1,soil_type)
1702          ENDIF
1703
1704          IF ( n_vangenuchten == 9999999.9_wp )  THEN
1705             n_vangenuchten = soil_pars(2,soil_type)           
1706          ENDIF
1707
1708          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
1709             hydraulic_conductivity = soil_pars(3,soil_type)           
1710          ENDIF
1711
1712          IF ( saturation_moisture == 9999999.9_wp )  THEN
1713             saturation_moisture = m_soil_pars(0,soil_type)           
1714          ENDIF
1715
1716          IF ( field_capacity == 9999999.9_wp )  THEN
1717             field_capacity = m_soil_pars(1,soil_type)           
1718          ENDIF
1719
1720          IF ( wilting_point == 9999999.9_wp )  THEN
1721             wilting_point = m_soil_pars(2,soil_type)           
1722          ENDIF
1723
1724          IF ( residual_moisture == 9999999.9_wp )  THEN
1725             residual_moisture = m_soil_pars(3,soil_type)       
1726          ENDIF
1727
1728       ENDIF   
1729
1730!
1731!--    Map values to the respective 2D arrays
1732       alpha_vg      = alpha_vangenuchten
1733       l_vg          = l_vangenuchten
1734       n_vg          = n_vangenuchten 
1735       gamma_w_sat   = hydraulic_conductivity
1736       m_sat         = saturation_moisture
1737       m_fc          = field_capacity
1738       m_wilt        = wilting_point
1739       m_res         = residual_moisture
1740       r_soil_min    = min_soil_resistance
1741
1742!
1743!--    Initial run actions
1744       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1745
1746          t_soil    = 0.0_wp
1747          m_liq_eb  = 0.0_wp
1748          m_soil    = 0.0_wp
1749
1750!
1751!--       Map user settings of T and q for each soil layer
1752!--       (make sure that the soil moisture does not drop below the permanent
1753!--       wilting point) -> problems with devision by zero)
1754          DO  k = nzb_soil, nzt_soil
1755             t_soil(k,:,:)    = soil_temperature(k)
1756             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
1757             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
1758          ENDDO
1759          t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
1760
1761!
1762!--       Calculate surface temperature
1763          t_surface   = pt_surface * exn
1764
1765!
1766!--       Set artifical values for ts and us so that r_a has its initial value
1767!--       for the first time step
1768          DO  i = nxlg, nxrg
1769             DO  j = nysg, nyng
1770                k = nzb_s_inner(j,i)
1771
1772                IF ( cloud_physics )  THEN
1773                   pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
1774                ELSE
1775                   pt1 = pt(k+1,j,i)
1776                ENDIF
1777
1778!
1779!--             Assure that r_a cannot be zero at model start
1780                IF ( pt1 == pt(k,j,i) )  pt1 = pt1 + 1.0E-10_wp
1781
1782                us(j,i)  = 0.1_wp
1783                ts(j,i)  = (pt1 - pt(k,j,i)) / r_a(j,i)
1784                shf(j,i) = - us(j,i) * ts(j,i)
1785             ENDDO
1786          ENDDO
1787
1788!
1789!--    Actions for restart runs
1790       ELSE
1791
1792          DO  i = nxlg, nxrg
1793             DO  j = nysg, nyng
1794                k = nzb_s_inner(j,i)               
1795                t_surface(j,i) = pt(k,j,i) * exn
1796             ENDDO
1797          ENDDO
1798
1799       ENDIF
1800
1801       DO  k = nzb_soil, nzt_soil
1802          root_fr(k,:,:) = root_fraction(k)
1803       ENDDO
1804
1805       IF ( veg_type /= 0 )  THEN
1806          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
1807             min_canopy_resistance = veg_pars(0,veg_type)
1808          ENDIF
1809          IF ( leaf_area_index == 9999999.9_wp )  THEN
1810             leaf_area_index = veg_pars(1,veg_type)         
1811          ENDIF
1812          IF ( vegetation_coverage == 9999999.9_wp )  THEN
1813             vegetation_coverage = veg_pars(2,veg_type)     
1814          ENDIF
1815          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
1816              canopy_resistance_coefficient= veg_pars(3,veg_type)     
1817          ENDIF
1818          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
1819             lambda_surface_stable = surface_pars(0,veg_type)         
1820          ENDIF
1821          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
1822             lambda_surface_unstable = surface_pars(1,veg_type)       
1823          ENDIF
1824          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
1825             f_shortwave_incoming = surface_pars(2,veg_type)       
1826          ENDIF
1827          IF ( z0_eb == 9999999.9_wp )  THEN
1828             roughness_length = roughness_par(0,veg_type) 
1829             z0_eb            = roughness_par(0,veg_type) 
1830          ENDIF
1831          IF ( z0h_eb == 9999999.9_wp )  THEN
1832             z0h_eb = roughness_par(1,veg_type)
1833          ENDIF
1834          IF ( z0q_eb == 9999999.9_wp )  THEN
1835             z0q_eb = roughness_par(2,veg_type)
1836          ENDIF
1837          z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp )
1838
1839          IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
1840             DO  k = nzb_soil, nzt_soil
1841                root_fr(k,:,:) = root_distribution(k,veg_type)
1842                root_fraction(k) = root_distribution(k,veg_type)
1843             ENDDO
1844          ENDIF
1845
1846       ELSE
1847
1848          IF ( z0_eb == 9999999.9_wp )  THEN
1849             z0_eb = roughness_length
1850          ENDIF
1851          IF ( z0h_eb == 9999999.9_wp )  THEN
1852             z0h_eb = z0_eb * z0h_factor
1853          ENDIF
1854          IF ( z0q_eb == 9999999.9_wp )  THEN
1855             z0q_eb = z0_eb * z0h_factor
1856          ENDIF
1857
1858       ENDIF
1859
1860!
1861!--    For surfaces covered with pavement, set depth of the pavement (with dry
1862!--    soil below). The depth must be greater than the first soil layer depth
1863       IF ( veg_type == 20 )  THEN
1864          IF ( pave_depth == 9999999.9_wp )  THEN
1865             pave_depth = zs(nzb_soil) 
1866          ELSE
1867             pave_depth = MAX( zs(nzb_soil), pave_depth )
1868          ENDIF
1869       ENDIF
1870
1871!
1872!--    Map vegetation and soil types to 2D array to allow for heterogeneous
1873!--    surfaces via user interface see below
1874       veg_type_2d = veg_type
1875       soil_type_2d = soil_type
1876
1877!
1878!--    Map vegetation parameters to the respective 2D arrays
1879       r_canopy_min         = min_canopy_resistance
1880       lai                  = leaf_area_index
1881       c_veg                = vegetation_coverage
1882       g_d                  = canopy_resistance_coefficient
1883       lambda_surface_s     = lambda_surface_stable
1884       lambda_surface_u     = lambda_surface_unstable
1885       f_sw_in              = f_shortwave_incoming
1886       z0                   = z0_eb
1887       z0h                  = z0h_eb
1888       z0q                  = z0q_eb
1889
1890!
1891!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
1892       CALL user_init_land_surface
1893
1894!
1895!--    Set flag parameter if vegetation type was set to a water surface. Also
1896!--    set temperature to a constant value in all "soil" layers.
1897       DO  i = nxlg, nxrg
1898          DO  j = nysg, nyng
1899             IF ( veg_type_2d(j,i) == 14  .OR.  veg_type_2d(j,i) == 15 )  THEN
1900                water_surface(j,i) = .TRUE.
1901                t_soil(:,j,i) = t_surface(j,i)
1902             ELSEIF ( veg_type_2d(j,i) == 20 )  THEN
1903                pave_surface(j,i) = .TRUE.
1904                m_soil(:,j,i) = 0.0_wp
1905             ENDIF
1906
1907          ENDDO
1908       ENDDO
1909
1910!
1911!--    Calculate new roughness lengths (for water surfaces only)
1912       CALL calc_z0_water_surface
1913
1914       t_soil_p    = t_soil
1915       m_soil_p    = m_soil
1916       m_liq_eb_p  = m_liq_eb
1917       t_surface_p = t_surface
1918
1919
1920
1921!--    Store initial profiles of t_soil and m_soil (assuming they are
1922!--    horizontally homogeneous on this PE)
1923       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
1924                                                nysg,nxlg), 2,                 &
1925                                                statistic_regions+1 )
1926       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
1927                                                nysg,nxlg), 2,                 &
1928                                                statistic_regions+1 )
1929
1930    END SUBROUTINE lsm_init
1931
1932
1933!------------------------------------------------------------------------------!
1934! Description:
1935! ------------
1936!> Allocate land surface model arrays and define pointers
1937!------------------------------------------------------------------------------!
1938    SUBROUTINE lsm_init_arrays
1939   
1940
1941       IMPLICIT NONE
1942
1943!
1944!--    Allocate surface and soil temperature / humidity
1945#if defined( __nopointer )
1946       ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
1947       ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
1948       ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1949       ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1950       ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
1951       ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
1952       ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
1953       ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
1954#else
1955       ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
1956       ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
1957       ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1958       ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1959       ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
1960       ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
1961       ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
1962       ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
1963#endif
1964
1965!
1966!--    Allocate intermediate timestep arrays
1967       ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
1968       ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1969       ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
1970       ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
1971
1972!
1973!--    Allocate 2D vegetation model arrays
1974       ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
1975       ALLOCATE ( building_surface(nysg:nyng,nxlg:nxrg) )
1976       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
1977       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
1978       ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
1979       ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
1980       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
1981       ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
1982       ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
1983       ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
1984       ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
1985       ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
1986       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
1987       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
1988       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
1989       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
1990       ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
1991       ALLOCATE ( pave_surface(nysg:nyng,nxlg:nxrg) )
1992       ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
1993       ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
1994       ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
1995       ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
1996       ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
1997       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
1998       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
1999       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
2000       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
2001       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
2002       ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
2003       ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
2004       ALLOCATE ( soil_type_2d(nysg:nyng,nxlg:nxrg) )
2005       ALLOCATE ( veg_type_2d(nysg:nyng,nxlg:nxrg) )
2006       ALLOCATE ( water_surface(nysg:nyng,nxlg:nxrg) )
2007
2008#if ! defined( __nopointer )
2009!
2010!--    Initial assignment of the pointers
2011       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
2012       t_surface => t_surface_1; t_surface_p => t_surface_2
2013       m_soil    => m_soil_1;    m_soil_p    => m_soil_2
2014       m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
2015#endif
2016
2017
2018    END SUBROUTINE lsm_init_arrays
2019
2020
2021!------------------------------------------------------------------------------!
2022! Description:
2023! ------------
2024!> Parin for &lsmpar for land surface model
2025!------------------------------------------------------------------------------!
2026    SUBROUTINE lsm_parin
2027
2028
2029       IMPLICIT NONE
2030
2031       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
2032       
2033       NAMELIST /lsm_par/         alpha_vangenuchten, c_surface,               &
2034                                  canopy_resistance_coefficient,               &
2035                                  conserve_water_content,                      &
2036                                  f_shortwave_incoming, field_capacity,        & 
2037                                  hydraulic_conductivity,                      &
2038                                  lambda_surface_stable,                       &
2039                                  lambda_surface_unstable, leaf_area_index,    &
2040                                  l_vangenuchten, min_canopy_resistance,       &
2041                                  min_soil_resistance, n_vangenuchten,         &
2042                                  pave_depth, pave_heat_capacity,              &
2043                                  pave_heat_conductivity,                      &
2044                                  residual_moisture, root_fraction,            &
2045                                  saturation_moisture, skip_time_do_lsm,       &
2046                                  soil_moisture, soil_temperature, soil_type,  &
2047                                  vegetation_coverage, veg_type, wilting_point,& 
2048                                  zs, z0_eb, z0h_eb, z0q_eb
2049       
2050       line = ' '
2051       
2052!
2053!--    Try to find land surface model package
2054       REWIND ( 11 )
2055       line = ' '
2056       DO   WHILE ( INDEX( line, '&lsm_par' ) == 0 )
2057          READ ( 11, '(A)', END=10 )  line
2058       ENDDO
2059       BACKSPACE ( 11 )
2060
2061!
2062!--    Read user-defined namelist
2063       READ ( 11, lsm_par )
2064
2065!
2066!--    Set flag that indicates that the land surface model is switched on
2067       land_surface = .TRUE.
2068
2069 10    CONTINUE
2070       
2071
2072    END SUBROUTINE lsm_parin
2073
2074
2075!------------------------------------------------------------------------------!
2076! Description:
2077! ------------
2078!> Soil model as part of the land surface model. The model predicts soil
2079!> temperature and water content.
2080!------------------------------------------------------------------------------!
2081    SUBROUTINE lsm_soil_model
2082
2083
2084       IMPLICIT NONE
2085
2086       INTEGER(iwp) ::  i   !< running index
2087       INTEGER(iwp) ::  j   !< running index
2088       INTEGER(iwp) ::  k   !< running index
2089
2090       REAL(wp)     :: h_vg !< Van Genuchten coef. h
2091
2092       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
2093                                                 lambda_temp, & !< temp. lambda
2094                                                 tend           !< tendency
2095
2096       DO  i = nxlg, nxrg
2097          DO  j = nysg, nyng
2098
2099             IF ( pave_surface(j,i) )  THEN
2100                rho_c_total(nzb_soil,j,i) = pave_heat_capacity
2101                lambda_temp(nzb_soil)     = pave_heat_conductivity
2102             ENDIF
2103
2104             IF (  .NOT.  water_surface(j,i) )  THEN
2105                DO  k = nzb_soil, nzt_soil
2106
2107
2108                   IF ( pave_surface(j,i)  .AND.  zs(k) <= pave_depth )  THEN
2109                   
2110                      rho_c_total(k,j,i) = pave_heat_capacity
2111                      lambda_temp(k)     = pave_heat_conductivity   
2112
2113                   ELSE           
2114!
2115!--                   Calculate volumetric heat capacity of the soil, taking
2116!--                   into account water content
2117                      rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) &
2118                                           + rho_c_water * m_soil(k,j,i))
2119
2120!
2121!--                   Calculate soil heat conductivity at the center of the soil
2122!--                   layers
2123                      lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *    &
2124                                     lambda_h_water ** m_soil(k,j,i)
2125
2126                      ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i)             &
2127                                                     / m_sat(j,i)))
2128
2129                      lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +    &
2130                                       lambda_h_dry
2131                   ENDIF
2132
2133                ENDDO
2134
2135!
2136!--             Calculate soil heat conductivity (lambda_h) at the _stag level
2137!--             using linear interpolation. For pavement surface, the
2138!--             true pavement depth is considered
2139                DO  k = nzb_soil, nzt_soil-1
2140                   IF ( pave_surface(j,i)  .AND.  zs(k)   < pave_depth         &
2141                                           .AND.  zs(k+1) > pave_depth )  THEN
2142                      lambda_h(k,j,i) = ( pave_depth - zs(k) ) / dz_soil(k+1)  &
2143                                        * lambda_temp(k)                       &
2144                                        + ( 1.0_wp - ( pave_depth - zs(k) )    &
2145                                        / dz_soil(k+1) ) * lambda_temp(k+1)
2146                   ELSE
2147                      lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
2148                                        * 0.5_wp
2149                   ENDIF
2150                ENDDO
2151                lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
2152
2153
2154
2155
2156!
2157!--             Prognostic equation for soil temperature t_soil
2158                tend(:) = 0.0_wp
2159
2160                tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *          &
2161                          ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)  &
2162                            - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)    &
2163                            + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
2164
2165                DO  k = nzb_soil+1, nzt_soil
2166                   tend(k) = (1.0_wp/rho_c_total(k,j,i))                       &
2167                             * (   lambda_h(k,j,i)                             &
2168                                 * ( t_soil(k+1,j,i) - t_soil(k,j,i) )         &
2169                                 * ddz_soil(k+1)                               &
2170                                 - lambda_h(k-1,j,i)                           &
2171                                 * ( t_soil(k,j,i) - t_soil(k-1,j,i) )         &
2172                                 * ddz_soil(k)                                 &
2173                               ) * ddz_soil_stag(k)
2174
2175                ENDDO
2176
2177                t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)&
2178                                                  + dt_3d * ( tsc(2)           &
2179                                                  * tend(nzb_soil:nzt_soil)    & 
2180                                                  + tsc(3)                     &
2181                                                  * tt_soil_m(:,j,i) )   
2182
2183!
2184!--             Calculate t_soil tendencies for the next Runge-Kutta step
2185                IF ( timestep_scheme(1:5) == 'runge' )  THEN
2186                   IF ( intermediate_timestep_count == 1 )  THEN
2187                      DO  k = nzb_soil, nzt_soil
2188                         tt_soil_m(k,j,i) = tend(k)
2189                      ENDDO
2190                   ELSEIF ( intermediate_timestep_count <                      &
2191                            intermediate_timestep_count_max )  THEN
2192                      DO  k = nzb_soil, nzt_soil
2193                         tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
2194                                         * tt_soil_m(k,j,i)
2195                      ENDDO
2196                   ENDIF
2197                ENDIF
2198
2199
2200                DO  k = nzb_soil, nzt_soil
2201
2202!
2203!--                Calculate soil diffusivity at the center of the soil layers
2204                   lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat       &
2205                                     / m_sat(j,i) ) * ( MAX( m_soil(k,j,i),    &
2206                                     m_wilt(j,i) ) / m_sat(j,i) )**(           &
2207                                     b_ch + 2.0_wp )
2208
2209!
2210!--                Parametrization of Van Genuchten
2211                   IF ( soil_type /= 7 )  THEN
2212!
2213!--                   Calculate the hydraulic conductivity after Van Genuchten
2214!--                   (1980)
2215                      h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -    &
2216                                 MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**(      &
2217                                 n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp  &
2218                             )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i)
2219
2220
2221                      gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +         &
2222                                      ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**(  &
2223                                      1.0_wp - 1.0_wp / n_vg(j,i) ) - (        &
2224                                      alpha_vg(j,i) * h_vg )**( n_vg(j,i)      &
2225                                      - 1.0_wp) )**2 )                         &
2226                                      / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg    &
2227                                      )**n_vg(j,i) )**( ( 1.0_wp  - 1.0_wp     &
2228                                      / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) )
2229
2230!
2231!--                Parametrization of Clapp & Hornberger
2232                   ELSE
2233                      gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i)       &
2234                                      / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
2235                   ENDIF
2236
2237                ENDDO
2238
2239!
2240!--             Prognostic equation for soil moisture content. Only performed,
2241!--             when humidity is enabled in the atmosphere and the surface type
2242!--             is not pavement (implies dry soil below).
2243                IF ( humidity  .AND.  .NOT.  pave_surface(j,i) )  THEN
2244!
2245!--                Calculate soil diffusivity (lambda_w) at the _stag level
2246!--                using linear interpolation. To do: replace this with
2247!--                ECMWF-IFS Eq. 8.81
2248                   DO  k = nzb_soil, nzt_soil-1
2249                     
2250                      lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
2251                                        * 0.5_wp
2252                      gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )    &
2253                                        * 0.5_wp
2254
2255                   ENDDO
2256
2257!
2258!
2259!--                In case of a closed bottom (= water content is conserved),
2260!--                set hydraulic conductivity to zero to that no water will be
2261!--                lost in the bottom layer.
2262                   IF ( conserve_water_content )  THEN
2263                      gamma_w(nzt_soil,j,i) = 0.0_wp
2264                   ELSE
2265                      gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
2266                   ENDIF     
2267
2268!--                The root extraction (= root_extr * qsws_veg_eb / (rho_l     
2269!--                * l_v)) ensures the mass conservation for water. The         
2270!--                transpiration of plants equals the cumulative withdrawals by
2271!--                the roots in the soil. The scheme takes into account the
2272!--                availability of water in the soil layers as well as the root
2273!--                fraction in the respective layer. Layer with moisture below
2274!--                wilting point will not contribute, which reflects the
2275!--                preference of plants to take water from moister layers.
2276
2277!
2278!--                Calculate the root extraction (ECMWF 7.69, the sum of
2279!--                root_extr = 1). The energy balance solver guarantees a
2280!--                positive transpiration, so that there is no need for an
2281!--                additional check.
2282                   m_total = 0.0_wp
2283                   DO  k = nzb_soil, nzt_soil
2284                       IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
2285                          m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
2286                       ENDIF
2287                   ENDDO 
2288
2289                   IF ( m_total > 0.0_wp )  THEN
2290                      DO  k = nzb_soil, nzt_soil
2291                         IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
2292                            root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i)      &
2293                                                            / m_total
2294                         ELSE
2295                            root_extr(k) = 0.0_wp
2296                         ENDIF
2297                      ENDDO
2298                   ENDIF
2299
2300!
2301!--                Prognostic equation for soil water content m_soil.
2302                   tend(:) = 0.0_wp
2303
2304                   tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (               &
2305                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
2306                            * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
2307                               root_extr(nzb_soil) * qsws_veg_eb(j,i)          &
2308                               + qsws_soil_eb(j,i) ) * drho_l_lv )             &
2309                               * ddz_soil_stag(nzb_soil)
2310
2311                   DO  k = nzb_soil+1, nzt_soil-1
2312                      tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)          &
2313                                - m_soil(k,j,i) ) * ddz_soil(k+1)              &
2314                                - gamma_w(k,j,i)                               &
2315                                - lambda_w(k-1,j,i) * (m_soil(k,j,i) -         &
2316                                m_soil(k-1,j,i)) * ddz_soil(k)                 &
2317                                + gamma_w(k-1,j,i) - (root_extr(k)             &
2318                                * qsws_veg_eb(j,i) * drho_l_lv)                &
2319                                ) * ddz_soil_stag(k)
2320
2321                   ENDDO
2322                   tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                  &
2323                                           - lambda_w(nzt_soil-1,j,i)          &
2324                                           * (m_soil(nzt_soil,j,i)             &
2325                                           - m_soil(nzt_soil-1,j,i))           &
2326                                           * ddz_soil(nzt_soil)                &
2327                                           + gamma_w(nzt_soil-1,j,i) - (       &
2328                                             root_extr(nzt_soil)               &
2329                                           * qsws_veg_eb(j,i) * drho_l_lv  )   &
2330                                     ) * ddz_soil_stag(nzt_soil)             
2331
2332                   m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
2333                                                   + dt_3d * ( tsc(2) * tend(:)   &
2334                                                   + tsc(3) * tm_soil_m(:,j,i) )   
2335   
2336!
2337!--                Account for dry soils (find a better solution here!)
2338                   DO  k = nzb_soil, nzt_soil
2339                      IF ( m_soil_p(k,j,i) < 0.0_wp )  m_soil_p(k,j,i) = 0.0_wp
2340                   ENDDO
2341
2342!
2343!--                Calculate m_soil tendencies for the next Runge-Kutta step
2344                   IF ( timestep_scheme(1:5) == 'runge' )  THEN
2345                      IF ( intermediate_timestep_count == 1 )  THEN
2346                         DO  k = nzb_soil, nzt_soil
2347                            tm_soil_m(k,j,i) = tend(k)
2348                         ENDDO
2349                      ELSEIF ( intermediate_timestep_count <                   &
2350                               intermediate_timestep_count_max )  THEN
2351                         DO  k = nzb_soil, nzt_soil
2352                            tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp&
2353                                     * tm_soil_m(k,j,i)
2354                         ENDDO
2355                      ENDIF
2356                   ENDIF
2357
2358                ENDIF
2359
2360             ENDIF
2361
2362          ENDDO
2363       ENDDO
2364
2365    END SUBROUTINE lsm_soil_model
2366
2367 
2368!------------------------------------------------------------------------------!
2369! Description:
2370! ------------
2371!> Swapping of timelevels
2372!------------------------------------------------------------------------------!
2373    SUBROUTINE lsm_swap_timelevel ( mod_count )
2374
2375       IMPLICIT NONE
2376
2377       INTEGER, INTENT(IN) :: mod_count
2378
2379#if defined( __nopointer )
2380
2381       t_surface    = t_surface_p
2382       t_soil       = t_soil_p
2383       IF ( humidity )  THEN
2384          m_soil    = m_soil_p
2385          m_liq_eb  = m_liq_eb_p
2386       ENDIF
2387
2388#else
2389   
2390       SELECT CASE ( mod_count )
2391
2392          CASE ( 0 )
2393
2394             t_surface  => t_surface_1; t_surface_p  => t_surface_2
2395             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
2396             IF ( humidity )  THEN
2397                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
2398                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
2399             ENDIF
2400
2401
2402          CASE ( 1 )
2403
2404             t_surface  => t_surface_2; t_surface_p  => t_surface_1
2405             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
2406             IF ( humidity )  THEN
2407                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
2408                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
2409             ENDIF
2410
2411       END SELECT
2412#endif
2413
2414    END SUBROUTINE lsm_swap_timelevel
2415
2416
2417!------------------------------------------------------------------------------!
2418! Description:
2419! ------------
2420!> Calculation of roughness length for open water (lakes, ocean). The
2421!> parameterization follows Charnock (1955). Two different implementations
2422!> are available: as in ECMWF-IFS (Beljaars 1994) or as in FLake (Subin et al.
2423!> 2012)
2424!------------------------------------------------------------------------------!
2425    SUBROUTINE calc_z0_water_surface
2426
2427       USE control_parameters,                                                 &
2428           ONLY: g, kappa, molecular_viscosity
2429
2430       IMPLICIT NONE
2431
2432       INTEGER :: i  !< running index
2433       INTEGER :: j  !< running index
2434
2435       REAL(wp), PARAMETER :: alpha_ch  = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF
2436!       REAL(wp), PARAMETER :: pr_number = 0.71_wp !< molecular Prandtl number in the Charnock parameterization (differs from prandtl_number)
2437!       REAL(wp), PARAMETER :: sc_number = 0.66_wp !< molecular Schmidt number in the Charnock parameterization
2438!       REAL(wp) :: re_0 !< near-surface roughness Reynolds number
2439
2440
2441       DO  i = nxlg, nxrg   
2442          DO  j = nysg, nyng
2443             IF ( water_surface(j,i) )  THEN
2444
2445!
2446!--             Disabled: FLake parameterization. Ideally, the Charnock
2447!--             coefficient should depend on the water depth and the fetch
2448!--             length
2449!                re_0 = z0(j,i) * us(j,i) / molecular_viscosity
2450!       
2451!                z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
2452!                              alpha_ch * us(j,i) / g )
2453!
2454!                z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) )
2455!                z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) )
2456
2457!
2458!--              Set minimum roughness length for u* > 0.2
2459!                IF ( us(j,i) > 0.2_wp )  THEN
2460!                   z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) )
2461!                   z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) )
2462!                ENDIF
2463
2464!
2465!--             ECMWF IFS model parameterization after Beljaars (1994). At low
2466!--             wind speed, the sea surface becomes aerodynamically smooth and
2467!--             the roughness scales with the viscosity. At high wind speed, the
2468!--             Charnock relation is used.
2469                z0(j,i) =   ( 0.11_wp * molecular_viscosity / us(j,i) )        &
2470                          + ( alpha_ch * us(j,i)**2 / g )
2471
2472                z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i)
2473                z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i)
2474
2475             ENDIF
2476          ENDDO
2477       ENDDO
2478
2479    END SUBROUTINE calc_z0_water_surface
2480
2481
2482!------------------------------------------------------------------------------!
2483! Description:
2484! ------------
2485!> Calculation of specific humidity of the skin layer (surface). It is assumend
2486!> that the skin is always saturated.
2487!------------------------------------------------------------------------------!
2488    SUBROUTINE calc_q_surface
2489
2490       IMPLICIT NONE
2491
2492       INTEGER :: i              !< running index
2493       INTEGER :: j              !< running index
2494       INTEGER :: k              !< running index
2495
2496       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
2497
2498       DO  i = nxlg, nxrg   
2499          DO  j = nysg, nyng
2500             k = nzb_s_inner(j,i)
2501
2502!
2503!--          Calculate water vapour pressure at saturation
2504             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface_p(j,i)   &
2505                   - 273.16_wp ) / ( t_surface_p(j,i) - 35.86_wp ) )
2506
2507!
2508!--          Calculate specific humidity at saturation
2509             q_s = 0.622_wp * e_s / surface_pressure
2510
2511             resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
2512
2513!
2514!--          Calculate specific humidity at surface
2515             IF ( cloud_physics )  THEN
2516                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
2517                             * ( q(k+1,j,i) - ql(k+1,j,i) )
2518             ELSE
2519                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
2520                             * q(k+1,j,i)
2521             ENDIF
2522
2523!
2524!--          Update virtual potential temperature
2525             vpt(k,j,i) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
2526
2527          ENDDO
2528       ENDDO
2529
2530    END SUBROUTINE calc_q_surface
2531
2532
2533 END MODULE land_surface_model_mod
Note: See TracBrowser for help on using the repository browser.