source: palm/trunk/SOURCE/land_surface_model.f90 @ 1682

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

Code annotations made doxygen readable

  • Property svn:keywords set to Id
File size: 71.5 KB
Line 
1!> @file land_surface_model.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-2014 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21! Code annotations made doxygen readable
22!
23! Former revisions:
24! -----------------
25! $Id: land_surface_model.f90 1682 2015-10-07 23:56:08Z knoop $
26! Bugfix: REAL constants provided with KIND-attribute in call of
27! 1590 2015-05-08 13:56:27Z maronga
28! Bugfix: definition of character strings requires same length for all elements
29!
30! 1585 2015-04-30 07:05:52Z maronga
31! Modifications for RRTMG. Changed tables to PARAMETER type.
32!
33! 1571 2015-03-12 16:12:49Z maronga
34! Removed upper-case variable names. Corrected distribution of precipitation to
35! the liquid water reservoir and the bare soil fractions.
36!
37! 1555 2015-03-04 17:44:27Z maronga
38! Added output of r_a and r_s
39!
40! 1553 2015-03-03 17:33:54Z maronga
41! Improved better treatment of roughness lengths. Added default soil temperature
42! profile
43!
44! 1551 2015-03-03 14:18:16Z maronga
45! Flux calculation is now done in prandtl_fluxes. Added support for data output.
46! Vertical indices have been replaced. Restart runs are now possible. Some
47! variables have beem renamed. Bugfix in the prognostic equation for the surface
48! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
49! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
50! the hydraulic conductivity. Bugfix for root fraction and extraction
51! calculation
52!
53! intrinsic function MAX and MIN
54!
55! 1500 2014-12-03 17:42:41Z maronga
56! Corrected calculation of aerodynamic resistance (r_a).
57! Precipitation is now added to liquid water reservoir using LE_liq.
58! Added support for dry runs.
59!
60! 1496 2014-12-02 17:25:50Z maronga
61! Initial revision
62!
63!
64! Description:
65! ------------
66!> Land surface model, consisting of a solver for the energy balance at the
67!> surface and a four layer soil scheme. The scheme is similar to the TESSEL
68!> scheme implemented in the ECMWF IFS model, with modifications according to
69!> H-TESSEL. The implementation is based on the formulation implemented in the
70!> DALES and UCLA-LES models.
71!>
72!> @todo Check dewfall parametrization for fog simulations.
73!> @todo Consider partial absorption of the net shortwave radiation by the surface layer.
74!> @todo Allow for water surfaces, check performance for bare soils.
75!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, nzt_soil=3)).
76!> @todo Implement surface runoff model (required when performing long-term LES with
77!>       considerable precipitation.
78!>
79!> @note  No time step criterion is required as long as the soil layers do not become
80!>        too thin.
81!------------------------------------------------------------------------------!
82 MODULE land_surface_model_mod
83 
84     USE arrays_3d,                                                            &
85         ONLY:  pt, pt_p, q_p, qsws, rif, shf, ts, us, z0, z0h
86
87     USE cloud_parameters,                                                     &
88         ONLY:  cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
89
90     USE control_parameters,                                                   &
91         ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,   &
92                initializing_actions, intermediate_timestep_count_max,         &
93                max_masks, precipitation, pt_surface, rho_surface,             &
94                roughness_length, surface_pressure, timestep_scheme, tsc,      &
95                z0h_factor, time_since_reference_point
96
97     USE indices,                                                              &
98         ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner 
99
100     USE kinds
101
102    USE netcdf_control,                                                        &
103        ONLY:  dots_label, dots_num, dots_unit
104
105     USE radiation_model_mod,                                                  &
106         ONLY:  radiation_scheme, rad_net, rad_sw_in, sigma_sb
107
108     USE statistics,                                                           &
109         ONLY:  hom, statistic_regions
110
111    IMPLICIT NONE
112
113!
114!-- LSM model constants
115    INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !< bottom of the soil model (to be switched)
116                               nzt_soil = 3, & !< top of the soil model (to be switched)
117                               nzs = 4         !< number of soil layers (fixed for now)
118
119    INTEGER(iwp) :: dots_soil = 0  !< starting index for timeseries output
120
121    INTEGER(iwp), DIMENSION(0:1) :: id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz,  &
122                                    id_dim_zs_3d, id_var_zs_xy,                & 
123                                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d
124                                   
125    INTEGER(iwp), DIMENSION(1:max_masks,0:1) :: id_dim_zs_mask, id_var_zs_mask
126
127    REAL(wp), PARAMETER ::                     &
128              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
129              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil   
130              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
131              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
132              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
133              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil
134              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water
135              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
136
137
138!
139!-- LSM variables
140    INTEGER(iwp) :: veg_type  = 2, & !< vegetation type, 0: user-defined, 1-19: generic (see list)
141                    soil_type = 3    !< soil type, 0: user-defined, 1-6: generic (see list)
142
143    LOGICAL :: conserve_water_content = .TRUE., & !< open or closed bottom surface for the soil model
144               dewfall = .TRUE.,                & !< allow/inhibit dewfall
145               land_surface = .FALSE.             !< flag parameter indicating wheather the lsm is used
146
147!   value 9999999.9_wp -> generic available or user-defined value must be set
148!   otherwise -> no generic variable and user setting is optional
149    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
150                canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
151                c_surface   = 20000.0_wp,               & !< Surface (skin) heat capacity
152                drho_l_lv,                              & !< (rho_l * l_v)**-1
153                exn,                                    & !< value of the Exner function
154                e_s = 0.0_wp,                           & !< saturation water vapour pressure
155                field_capacity = 9999999.9_wp,          & !< NAMELIST m_fc
156                f_shortwave_incoming = 9999999.9_wp,    & !< NAMELIST f_sw_in
157                hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
158                ke = 0.0_wp,                            & !< Kersten number
159                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s
160                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u
161                leaf_area_index = 9999999.9_wp,         & !< NAMELIST lai
162                l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
163                min_canopy_resistance = 9999999.9_wp,   & !< NAMELIST r_canopy_min
164                min_soil_resistance = 50.0_wp,          & !< NAMELIST r_soil_min
165                m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
166                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
167                q_s = 0.0_wp,                           & !< saturation specific humidity
168                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
169                rho_cp,                                 & !< rho_surface * cp
170                rho_lv,                                 & !< rho * l_v
171                rd_d_rv,                                & !< r_d / r_v
172                saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
173                vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
174                wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
175                z0_eb  = 9999999.9_wp,                  & !< NAMELIST z0 (lsm_par)
176                z0h_eb = 9999999.9_wp                    !< NAMELIST z0h (lsm_par)
177
178    REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
179              ddz_soil,                     & !< 1/dz_soil
180              ddz_soil_stag,                & !< 1/dz_soil_stag
181              dz_soil,                      & !< soil grid spacing (center-center)
182              dz_soil_stag,                 & !< soil grid spacing (edge-edge)
183              root_extr = 0.0_wp,           & !< root extraction
184              root_fraction = (/9999999.9_wp, 9999999.9_wp,    &
185                                9999999.9_wp, 9999999.9_wp /), & !< distribution of root surface area to the individual soil layers
186              zs = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),    & !< soil layer depths (m)
187              soil_moisture = 0.0_wp          !< soil moisture content (m3/m3)
188
189    REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) ::   &
190              soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp,  283.0_wp,    &
191                                   283.0_wp /) !< soil temperature (K)
192
193#if defined( __nopointer )
194    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface,   & !< surface temperature (K)
195                                                     t_surface_p, & !< progn. surface temperature (K)
196                                                     m_liq_eb,    & !< liquid water reservoir (m)
197                                                     m_liq_eb_av, & !< liquid water reservoir (m)
198                                                     m_liq_eb_p     !< progn. liquid water reservoir (m)
199#else
200    REAL(wp), DIMENSION(:,:), POINTER :: t_surface,      &
201                                         t_surface_p,    & 
202                                         m_liq_eb,       & 
203                                         m_liq_eb_p
204
205    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, &
206                                                     m_liq_eb_av,              &
207                                                     m_liq_eb_1, m_liq_eb_2
208#endif
209
210!
211!-- Temporal tendencies for time stepping           
212    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m,  & !< surface temperature tendency (K)
213                                             tm_liq_eb_m      !< liquid water reservoir tendency (m)
214
215!
216!-- Energy balance variables               
217    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::                                   &
218              alpha_vg,         & !< coef. of Van Genuchten
219              c_liq,            & !< liquid water coverage (of vegetated area)
220              c_liq_av,         & !< average of c_liq
221              c_soil_av,        & !< average of c_soil
222              c_veg,            & !< vegetation coverage
223              c_veg_av,         & !< average of c_veg
224              f_sw_in,          & !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
225              ghf_eb,           & !< ground heat flux
226              ghf_eb_av,        & !< average of ghf_eb
227              gamma_w_sat,      & !< hydraulic conductivity at saturation
228              g_d,              & !< coefficient for dependence of r_canopy on water vapour pressure deficit
229              lai,              & !< leaf area index
230              lai_av,           & !< average of lai
231              lambda_h_sat,     & !< heat conductivity for dry soil
232              lambda_surface_s,    & !< coupling between surface and soil (depends on vegetation type)
233              lambda_surface_u,    & !< coupling between surface and soil (depends on vegetation type)
234              l_vg,             & !< coef. of Van Genuchten
235              m_fc,             & !< soil moisture at field capacity (m3/m3)
236              m_res,            & !< residual soil moisture
237              m_sat,            & !< saturation soil moisture (m3/m3)
238              m_wilt,           & !< soil moisture at permanent wilting point (m3/m3)
239              n_vg,             & !< coef. Van Genuchten 
240              qsws_eb,          & !< surface flux of latent heat (total)
241              qsws_eb_av,       & !< average of qsws_eb
242              qsws_liq_eb,      & !< surface flux of latent heat (liquid water portion)
243              qsws_liq_eb_av,   & !< average of qsws_liq_eb
244              qsws_soil_eb,     & !< surface flux of latent heat (soil portion)
245              qsws_soil_eb_av,  & !< average of qsws_soil_eb
246              qsws_veg_eb,      & !< surface flux of latent heat (vegetation portion)
247              qsws_veg_eb_av,   & !< average of qsws_veg_eb
248              rad_net_l,        & !< local copy of rad_net (net radiation at surface)
249              r_a,              & !< aerodynamic resistance
250              r_a_av,           & !< avergae of r_a
251              r_canopy,         & !< canopy resistance
252              r_soil,           & !< soil resitance
253              r_soil_min,       & !< minimum soil resistance
254              r_s,              & !< total surface resistance (combination of r_soil and r_canopy)
255              r_s_av,           & !< avergae of r_s
256              r_canopy_min,     & !< minimum canopy (stomatal) resistance
257              shf_eb,           & !< surface flux of sensible heat
258              shf_eb_av           !< average of shf_eb
259
260    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
261              lambda_h, &   !< heat conductivity of soil (?)                           
262              lambda_w, &   !< hydraulic diffusivity of soil (?)
263              gamma_w,  &   !< hydraulic conductivity of soil (?)
264              rho_c_total   !< volumetric heat capacity of the actual soil matrix (?)
265
266#if defined( __nopointer )
267    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
268              t_soil,    & !< Soil temperature (K)
269              t_soil_av, & !< Average of t_soil
270              t_soil_p,  & !< Prog. soil temperature (K)
271              m_soil,    & !< Soil moisture (m3/m3)
272              m_soil_av, & !< Average of m_soil
273              m_soil_p     !< Prog. soil moisture (m3/m3)
274#else
275    REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
276              t_soil, t_soil_p, &
277              m_soil, m_soil_p   
278
279    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
280              t_soil_av, t_soil_1, t_soil_2,                                   &
281              m_soil_av, m_soil_1, m_soil_2
282#endif
283
284
285    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
286              tt_soil_m, & !< t_soil storage array
287              tm_soil_m, & !< m_soil storage array
288              root_fr      !< root fraction (sum=1)
289
290
291!
292!-- Predefined Land surface classes (veg_type)
293    CHARACTER(26), DIMENSION(0:19), PARAMETER :: veg_type_name = (/ &
294                                   'user defined              ',    & ! 0
295                                   'crops, mixed farming      ',    & !  1
296                                   'short grass               ',    & !  2
297                                   'evergreen needleleaf trees',    & !  3
298                                   'deciduous needleleaf trees',    & !  4
299                                   'evergreen broadleaf trees ',    & !  5
300                                   'deciduous broadleaf trees ',    & !  6
301                                   'tall grass                ',    & !  7
302                                   'desert                    ',    & !  8
303                                   'tundra                    ',    & !  9
304                                   'irrigated crops           ',    & ! 10
305                                   'semidesert                ',    & ! 11
306                                   'ice caps and glaciers     ',    & ! 12
307                                   'bogs and marshes          ',    & ! 13
308                                   'inland water              ',    & ! 14
309                                   'ocean                     ',    & ! 15
310                                   'evergreen shrubs          ',    & ! 16
311                                   'deciduous shrubs          ',    & ! 17
312                                   'mixed forest/woodland     ',    & ! 18
313                                   'interrupted forest        '     & ! 19
314                                                                 /)
315
316!
317!-- Soil model classes (soil_type)
318    CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
319                                   'user defined',                  & ! 0
320                                   'coarse      ',                  & ! 1
321                                   'medium      ',                  & ! 2
322                                   'medium-fine ',                  & ! 3
323                                   'fine        ',                  & ! 4
324                                   'very fine   ',                  & ! 5
325                                   'organic     ',                  & ! 6
326                                   'loamy (CH)  '                   & ! 7
327                                                                 /)
328!
329!-- Land surface parameters according to the respective classes (veg_type)
330
331!
332!-- Land surface parameters I
333!--                          r_canopy_min,     lai,   c_veg,     g_d
334    REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: veg_pars = RESHAPE( (/ &
335                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
336                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
337                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  3
338                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  4
339                                 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  5
340                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,  & !  6
341                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,  & !  7
342                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  & !  8
343                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  & !  9
344                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & ! 10
345                                 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,  & ! 11
346                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 12
347                                 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,  & ! 13
348                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 14
349                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 15
350                                 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,  & ! 16
351                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
352                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
353                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp   & ! 19
354                                 /), (/ 4, 19 /) )
355
356!
357!-- Land surface parameters II          z0,         z0h
358    REAL(wp), DIMENSION(0:1,1:19), PARAMETER :: roughness_par = RESHAPE( (/ & 
359                                   0.25_wp,  0.25E-2_wp,                    & !  1
360                                   0.20_wp,  0.20E-2_wp,                    & !  2
361                                   2.00_wp,     2.00_wp,                    & !  3
362                                   2.00_wp,     2.00_wp,                    & !  4
363                                   2.00_wp,     2.00_wp,                    & !  5
364                                   2.00_wp,     2.00_wp,                    & !  6
365                                   0.47_wp,  0.47E-2_wp,                    & !  7
366                                  0.013_wp, 0.013E-2_wp,                    & !  8
367                                  0.034_wp, 0.034E-2_wp,                    & !  9
368                                    0.5_wp,  0.50E-2_wp,                    & ! 10
369                                   0.17_wp,  0.17E-2_wp,                    & ! 11
370                                 1.3E-3_wp,   1.3E-4_wp,                    & ! 12
371                                   0.83_wp,  0.83E-2_wp,                    & ! 13
372                                   0.00_wp,  0.00E-2_wp,                    & ! 14
373                                   0.00_wp,  0.00E-2_wp,                    & ! 15
374                                   0.10_wp,  0.10E-2_wp,                    & ! 16
375                                   0.25_wp,  0.25E-2_wp,                    & ! 17
376                                   2.00_wp,  2.00E-2_wp,                    & ! 18
377                                   1.10_wp,  1.10E-2_wp                     & ! 19
378                                 /), (/ 2, 19 /) )
379
380!
381!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
382    REAL(wp), DIMENSION(0:2,1:19), PARAMETER :: surface_pars = RESHAPE( (/ &
383                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  1
384                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  2
385                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  3
386                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  4
387                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  5
388                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  6
389                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  7
390                                      15.0_wp,       15.0_wp, 0.00_wp,     & !  8
391                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  9
392                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 10
393                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 11
394                                      58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
395                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 13
396                                    1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 14
397                                    1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 15
398                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
399                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
400                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
401                                      20.0_wp,       15.0_wp, 0.03_wp      & ! 19
402                                      /), (/ 3, 19 /) )
403
404!
405!-- Root distribution (sum = 1)  level 1, level 2, level 3, level 4,
406    REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: root_distribution = RESHAPE( (/ &
407                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  1
408                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  2
409                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  3
410                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  4
411                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  5
412                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  6
413                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  7
414                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  8
415                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & !  9
416                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 10
417                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 11
418                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 12
419                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 13
420                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 14
421                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 15
422                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16
423                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 17
424                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 18
425                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 19
426                                 /), (/ 4, 19 /) )
427
428!
429!-- Soil parameters according to the following porosity classes (soil_type)
430
431!
432!-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
433    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/     &
434                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
435                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
436                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, & ! 3
437                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
438                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
439                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, & ! 6
440                                 0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  & ! 7
441                                 /), (/ 4, 7 /) )
442
443!
444!-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res 
445    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
446                                 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
447                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
448                                 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
449                                 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
450                                 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
451                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
452                                 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  & ! 7
453                                 /), (/ 4, 7 /) )
454
455
456    SAVE
457
458
459    PRIVATE
460
461
462!
463!-- Public parameters, constants and initial values
464    PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient,       &
465           conserve_water_content, dewfall, field_capacity,                    &
466           f_shortwave_incoming, hydraulic_conductivity, init_lsm,             &
467           init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable,    &
468           land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model,  &
469           lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance,          &
470           min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp,     &
471           rho_lv, root_fraction, saturation_moisture, soil_moisture,          &
472           soil_temperature, soil_type, soil_type_name, vegetation_coverage,   &
473           veg_type, veg_type_name, wilting_point, z0_eb, z0h_eb
474
475!
476!-- Public grid and NetCDF variables
477    PUBLIC dots_soil, id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz,                &
478           id_dim_zs_3d, id_dim_zs_mask, id_var_zs_xy, id_var_zs_xz,           &
479           id_var_zs_yz, id_var_zs_3d, id_var_zs_mask, nzb_soil, nzs, nzt_soil,&
480           zs
481
482
483!
484!-- Public 2D output variables
485    PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av,     &
486           lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av,      &
487           qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av,         &
488           r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av
489
490
491!
492!-- Public prognostic variables
493    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
494
495    INTERFACE init_lsm
496       MODULE PROCEDURE init_lsm
497    END INTERFACE init_lsm
498
499    INTERFACE lsm_energy_balance
500       MODULE PROCEDURE lsm_energy_balance
501    END INTERFACE lsm_energy_balance
502
503    INTERFACE lsm_soil_model
504       MODULE PROCEDURE lsm_soil_model
505    END INTERFACE lsm_soil_model
506
507    INTERFACE lsm_swap_timelevel
508       MODULE PROCEDURE lsm_swap_timelevel
509    END INTERFACE lsm_swap_timelevel
510
511 CONTAINS
512
513
514!------------------------------------------------------------------------------!
515! Description:
516! ------------
517!> Allocate land surface model arrays and define pointers
518!------------------------------------------------------------------------------!
519    SUBROUTINE init_lsm_arrays
520   
521
522       IMPLICIT NONE
523
524!
525!--    Allocate surface and soil temperature / humidity
526#if defined( __nopointer )
527       ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
528       ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
529       ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
530       ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
531       ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
532       ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
533       ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
534       ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
535#else
536       ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
537       ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
538       ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
539       ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
540       ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
541       ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
542       ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
543       ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
544#endif
545
546!
547!--    Allocate intermediate timestep arrays
548       ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
549       ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
550       ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
551       ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
552
553!
554!--    Allocate 2D vegetation model arrays
555       ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
556       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
557       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
558       ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
559       ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
560       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
561       ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
562       ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
563       ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
564       ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) )
565       ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
566       ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
567       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
568       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
569       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
570       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
571       ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
572       ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
573       ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
574       ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
575       ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
576       ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
577       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
578       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
579       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
580       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
581       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
582       ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
583       ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
584
585#if ! defined( __nopointer )
586!
587!--    Initial assignment of the pointers
588       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
589       t_surface => t_surface_1; t_surface_p => t_surface_2
590       m_soil    => m_soil_1;    m_soil_p    => m_soil_2
591       m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
592#endif
593
594
595    END SUBROUTINE init_lsm_arrays
596
597!------------------------------------------------------------------------------!
598! Description:
599! ------------
600!> Initialization of the land surface model
601!------------------------------------------------------------------------------!
602    SUBROUTINE init_lsm
603   
604
605       IMPLICIT NONE
606
607       INTEGER(iwp) ::  i !< running index
608       INTEGER(iwp) ::  j !< running index
609       INTEGER(iwp) ::  k !< running index
610
611
612!
613!--    Calculate Exner function
614       exn       = ( surface_pressure / 1000.0_wp )**0.286_wp
615
616
617!
618!--    If no cloud physics is used, rho_surface has not been calculated before
619       IF ( .NOT. cloud_physics )  THEN
620          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
621       ENDIF
622
623!
624!--    Calculate frequently used parameters
625       rho_cp    = cp * rho_surface
626       rd_d_rv   = r_d / r_v
627       rho_lv    = rho_surface * l_v
628       drho_l_lv = 1.0_wp / (rho_l * l_v)
629
630!
631!--    Set inital values for prognostic quantities
632       tt_surface_m = 0.0_wp
633       tt_soil_m    = 0.0_wp
634       tm_liq_eb_m  = 0.0_wp
635       c_liq        = 0.0_wp
636
637       ghf_eb = 0.0_wp
638       shf_eb = rho_cp * shf
639
640       IF ( humidity )  THEN
641          qsws_eb = rho_l * l_v * qsws
642       ELSE
643          qsws_eb = 0.0_wp
644       ENDIF
645
646       qsws_liq_eb  = 0.0_wp
647       qsws_soil_eb = qsws_eb
648       qsws_veg_eb  = 0.0_wp
649
650       r_a        = 50.0_wp
651       r_s        = 50.0_wp
652       r_canopy   = 0.0_wp
653       r_soil     = 0.0_wp
654
655!
656!--    Allocate 3D soil model arrays
657       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
658       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
659       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
660
661       lambda_h = 0.0_wp
662!
663!--    If required, allocate humidity-related variables for the soil model
664       IF ( humidity )  THEN
665          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
666          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
667
668          lambda_w = 0.0_wp 
669       ENDIF
670
671!
672!--    Calculate grid spacings. Temperature and moisture are defined at
673!--    the center of the soil layers, whereas gradients/fluxes are defined
674!--    at the edges (_stag)
675       dz_soil_stag(nzb_soil) = zs(nzb_soil)
676
677       DO k = nzb_soil+1, nzt_soil
678          dz_soil_stag(k) = zs(k) - zs(k-1)
679       ENDDO
680
681       DO k = nzb_soil, nzt_soil-1
682          dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k))
683       ENDDO
684       dz_soil(nzt_soil) = dz_soil_stag(nzt_soil)
685
686       ddz_soil      = 1.0_wp / dz_soil
687       ddz_soil_stag = 1.0_wp / dz_soil_stag
688
689!
690!--    Initialize standard soil types. It is possible to overwrite each
691!--    parameter by setting the respecticy NAMELIST variable to a
692!--    value /= 9999999.9.
693       IF ( soil_type /= 0 )  THEN 
694 
695          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
696             alpha_vangenuchten = soil_pars(0,soil_type)
697          ENDIF
698
699          IF ( l_vangenuchten == 9999999.9_wp )  THEN
700             l_vangenuchten = soil_pars(1,soil_type)
701          ENDIF
702
703          IF ( n_vangenuchten == 9999999.9_wp )  THEN
704             n_vangenuchten = soil_pars(2,soil_type)           
705          ENDIF
706
707          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
708             hydraulic_conductivity = soil_pars(3,soil_type)           
709          ENDIF
710
711          IF ( saturation_moisture == 9999999.9_wp )  THEN
712             saturation_moisture = m_soil_pars(0,soil_type)           
713          ENDIF
714
715          IF ( field_capacity == 9999999.9_wp )  THEN
716             field_capacity = m_soil_pars(1,soil_type)           
717          ENDIF
718
719          IF ( wilting_point == 9999999.9_wp )  THEN
720             wilting_point = m_soil_pars(2,soil_type)           
721          ENDIF
722
723          IF ( residual_moisture == 9999999.9_wp )  THEN
724             residual_moisture = m_soil_pars(3,soil_type)       
725          ENDIF
726
727       ENDIF   
728
729       alpha_vg      = alpha_vangenuchten
730       l_vg          = l_vangenuchten
731       n_vg          = n_vangenuchten 
732       gamma_w_sat   = hydraulic_conductivity
733       m_sat         = saturation_moisture
734       m_fc          = field_capacity
735       m_wilt        = wilting_point
736       m_res         = residual_moisture
737       r_soil_min    = min_soil_resistance
738
739!
740!--    Initial run actions
741       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
742
743          t_soil    = 0.0_wp
744          m_liq_eb  = 0.0_wp
745          m_soil    = 0.0_wp
746
747!
748!--       Map user settings of T and q for each soil layer
749!--       (make sure that the soil moisture does not drop below the permanent
750!--       wilting point) -> problems with devision by zero)
751          DO k = nzb_soil, nzt_soil
752             t_soil(k,:,:)    = soil_temperature(k)
753             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
754             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
755          ENDDO
756          t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
757
758!
759!--       Calculate surface temperature
760          t_surface = pt_surface * exn
761
762!
763!--       Set artifical values for ts and us so that r_a has its initial value for
764!--       the first time step
765          DO  i = nxlg, nxrg
766             DO  j = nysg, nyng
767                k = nzb_s_inner(j,i)
768
769!
770!--             Assure that r_a cannot be zero at model start
771                IF ( pt(k+1,j,i) == pt(k,j,i) )  pt(k+1,j,i) = pt(k+1,j,i)     &
772                                                 + 1.0E-10_wp
773
774                us(j,i) = 0.1_wp
775                ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i)
776                shf(j,i) = - us(j,i) * ts(j,i)
777             ENDDO
778          ENDDO
779
780!
781!--    Actions for restart runs
782       ELSE
783
784          DO  i = nxlg, nxrg
785             DO  j = nysg, nyng
786                k = nzb_s_inner(j,i)               
787                t_surface(j,i) = pt(k,j,i) * exn
788             ENDDO
789          ENDDO
790
791       ENDIF
792
793!
794!--    Calculate saturation soil heat conductivity
795       lambda_h_sat(:,:) = lambda_h_sm ** (1.0_wp - m_sat(:,:)) *              &
796                           lambda_h_water ** m_sat(:,:)
797
798
799
800
801       DO k = nzb_soil, nzt_soil
802          root_fr(k,:,:) = root_fraction(k)
803       ENDDO
804
805       IF ( veg_type /= 0 )  THEN
806          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
807             min_canopy_resistance = veg_pars(0,veg_type)
808          ENDIF
809          IF ( leaf_area_index == 9999999.9_wp )  THEN
810             leaf_area_index = veg_pars(1,veg_type)         
811          ENDIF
812          IF ( vegetation_coverage == 9999999.9_wp )  THEN
813             vegetation_coverage = veg_pars(2,veg_type)     
814          ENDIF
815          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
816              canopy_resistance_coefficient= veg_pars(3,veg_type)     
817          ENDIF
818          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
819             lambda_surface_stable = surface_pars(0,veg_type)         
820          ENDIF
821          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
822             lambda_surface_unstable = surface_pars(1,veg_type)       
823          ENDIF
824          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
825             f_shortwave_incoming = surface_pars(2,veg_type)       
826          ENDIF
827          IF ( z0_eb == 9999999.9_wp )  THEN
828             roughness_length = roughness_par(0,veg_type) 
829             z0_eb            = roughness_par(0,veg_type) 
830          ENDIF
831          IF ( z0h_eb == 9999999.9_wp )  THEN
832             z0h_eb = roughness_par(1,veg_type)
833          ENDIF
834          z0h_factor = z0h_eb / z0_eb
835
836          IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
837             DO k = nzb_soil, nzt_soil
838                root_fr(k,:,:) = root_distribution(k,veg_type)
839                root_fraction(k) = root_distribution(k,veg_type)
840             ENDDO
841          ENDIF
842
843       ELSE
844
845          IF ( z0_eb == 9999999.9_wp )  THEN
846             z0_eb = roughness_length
847          ENDIF
848          IF ( z0h_eb == 9999999.9_wp )  THEN
849             z0h_eb = z0_eb * z0h_factor
850          ENDIF
851
852       ENDIF
853
854!
855!--    Initialize vegetation
856       r_canopy_min         = min_canopy_resistance
857       lai                  = leaf_area_index
858       c_veg                = vegetation_coverage
859       g_d                  = canopy_resistance_coefficient
860       lambda_surface_s     = lambda_surface_stable
861       lambda_surface_u     = lambda_surface_unstable
862       f_sw_in              = f_shortwave_incoming
863       z0                   = z0_eb
864       z0h                  = z0h_eb
865
866!
867!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
868       CALL user_init_land_surface
869
870
871       t_soil_p    = t_soil
872       m_soil_p    = m_soil
873       m_liq_eb_p  = m_liq_eb
874
875!--    Store initial profiles of t_soil and m_soil (assuming they are
876!--    horizontally homogeneous on this PE)
877       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
878                                                nysg,nxlg), 2,                 &
879                                                statistic_regions+1 )
880       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
881                                                nysg,nxlg), 2,                 &
882                                                statistic_regions+1 )
883
884!
885!--    Calculate humidity at the surface
886       IF ( humidity )  THEN
887          CALL calc_q_surface
888       ENDIF
889
890
891
892!
893!--    Add timeseries for land surface model
894
895       dots_label(dots_num+1) = "ghf_eb"
896       dots_label(dots_num+2) = "shf_eb"
897       dots_label(dots_num+3) = "qsws_eb"
898       dots_label(dots_num+4) = "qsws_liq_eb"
899       dots_label(dots_num+5) = "qsws_soil_eb"
900       dots_label(dots_num+6) = "qsws_veg_eb"
901       dots_unit(dots_num+1:dots_num+6) = "W/m2"
902
903       dots_label(dots_num+7) = "r_a"
904       dots_label(dots_num+8) = "r_s"
905       dots_unit(dots_num+7:dots_num+8) = "s/m"
906
907       dots_soil = dots_num + 1
908       dots_num  = dots_num + 8
909
910       RETURN
911
912    END SUBROUTINE init_lsm
913
914
915
916!------------------------------------------------------------------------------!
917! Description:
918! ------------
919!> Solver for the energy balance at the surface.
920!> @note Surface fluxes are calculated in the land surface model, but these are
921!>       not used in the atmospheric code. The fluxes are calculated afterwards in
922!>       prandtl_fluxes using the surface values of temperature and humidity as
923!>       provided by the land surface model. In this way, the fluxes in the land
924!>       surface model are not equal to the ones calculated in prandtl_fluxes
925!------------------------------------------------------------------------------!
926    SUBROUTINE lsm_energy_balance
927
928
929       IMPLICIT NONE
930
931       INTEGER(iwp) ::  i         !< running index
932       INTEGER(iwp) ::  j         !< running index
933       INTEGER(iwp) ::  k, ks     !< running index
934
935       REAL(wp) :: f1,          & !< resistance correction term 1
936                   f2,          & !< resistance correction term 2
937                   f3,          & !< resistance correction term 3
938                   m_min,       & !< minimum soil moisture
939                   T_1,         & !< actual temperature at first grid point
940                   e,           & !< water vapour pressure
941                   e_s,         & !< water vapour saturation pressure
942                   e_s_dt,      & !< derivate of e_s with respect to T
943                   tend,        & !< tendency
944                   dq_s_dt,     & !< derivate of q_s with respect to T
945                   coef_1,      & !< coef. for prognostic equation
946                   coef_2,      & !< coef. for prognostic equation
947                   f_qsws,      & !< factor for qsws_eb
948                   f_qsws_veg,  & !< factor for qsws_veg_eb
949                   f_qsws_soil, & !< factor for qsws_soil_eb
950                   f_qsws_liq,  & !< factor for qsws_liq_eb
951                   f_shf,       & !< factor for shf_eb
952                   lambda_surface, & !< Current value of lambda_surface
953                   m_liq_eb_max   !< maxmimum value of the liq. water reservoir
954
955
956!
957!--    Calculate the exner function for the current time step
958       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
959
960       
961
962       DO i = nxlg, nxrg
963          DO j = nysg, nyng
964             k = nzb_s_inner(j,i)
965
966!
967!--          Set lambda_surface according to stratification
968             IF ( rif(j,i) >= 0.0_wp )  THEN
969                lambda_surface = lambda_surface_s(j,i)
970             ELSE
971                lambda_surface = lambda_surface_u(j,i)
972             ENDIF
973
974!
975!--          First step: calculate aerodyamic resistance. As pt, us, ts
976!--          are not available for the prognostic time step, data from the last
977!--          time step is used here. Note that this formulation is the
978!--          equivalent to the ECMWF formulation using drag coefficients
979             r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) +       &
980                         1.0E-20_wp)
981
982!
983!--          Second step: calculate canopy resistance r_canopy
984!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
985 
986!--          f1: correction for incoming shortwave radiation (stomata close at
987!--          night)
988             IF ( radiation_scheme /= 'constant' )  THEN
989                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /&
990                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)        &
991                               + 1.0_wp)) )
992             ELSE
993                f1 = 1.0_wp
994             ENDIF
995
996!
997!--          f2: correction for soil moisture availability to plants (the
998!--          integrated soil moisture must thus be considered here)
999!--          f2 = 0 for very dry soils
1000             m_total = 0.0_wp
1001             DO ks = nzb_soil, nzt_soil
1002                 m_total = m_total + root_fr(ks,j,i)                           &
1003                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
1004             ENDDO 
1005
1006             IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LT. m_fc(j,i) )  THEN
1007                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
1008             ELSEIF ( m_total .GE. m_fc(j,i) )  THEN
1009                f2 = 1.0_wp
1010             ELSE
1011                f2 = 1.0E-20_wp
1012             ENDIF
1013
1014!
1015!--          Calculate water vapour pressure at saturation
1016             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
1017                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
1018
1019!
1020!--          f3: correction for vapour pressure deficit
1021             IF ( g_d(j,i) /= 0.0_wp )  THEN
1022!
1023!--             Calculate vapour pressure
1024                e  = q_p(k+1,j,i) * surface_pressure / 0.622_wp
1025                f3 = EXP ( -g_d(j,i) * (e_s - e) )
1026             ELSE
1027                f3 = 1.0_wp
1028             ENDIF
1029
1030!
1031!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1032!--          this calculation is obsolete, as r_canopy is not used below.
1033!--          To do: check for very dry soil -> r_canopy goes to infinity
1034             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
1035                                             + 1.0E-20_wp)
1036
1037!
1038!--          Third step: calculate bare soil resistance r_soil. The Clapp &
1039!--          Hornberger parametrization does not consider c_veg.
1040             IF ( soil_type /= 7 )  THEN
1041                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
1042                        m_res(j,i)
1043             ELSE
1044                m_min = m_wilt(j,i)
1045             ENDIF
1046
1047             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
1048             f2 = MAX(f2,1.0E-20_wp)
1049             f2 = MIN(f2,1.0_wp)
1050
1051             r_soil(j,i) = r_soil_min(j,i) / f2
1052
1053!
1054!--          Calculate fraction of liquid water reservoir
1055             m_liq_eb_max = m_max_depth * lai(j,i)
1056             c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp))
1057             
1058
1059!
1060!--          Calculate saturation specific humidity
1061             q_s = 0.622_wp * e_s / surface_pressure
1062
1063!
1064!--          In case of dewfall, set evapotranspiration to zero
1065!--          All super-saturated water is then removed from the air
1066             IF ( humidity .AND. dewfall .AND. q_s .LE. q_p(k+1,j,i) )  THEN
1067                r_canopy(j,i) = 0.0_wp
1068                r_soil(j,i)   = 0.0_wp
1069             ENDIF 
1070
1071!
1072!--          Calculate coefficients for the total evapotranspiration
1073             f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
1074                           (r_a(j,i) + r_canopy(j,i))
1075
1076             f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
1077                                                          r_soil(j,i))
1078             f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
1079
1080
1081!
1082!--          If soil moisture is below wilting point, plants do no longer
1083!--          transpirate.
1084!              IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
1085!                 f_qsws_veg = 0.0_wp
1086!              ENDIF
1087
1088!
1089!--          If dewfall is deactivated, vegetation, soil and liquid water
1090!--          reservoir are not allowed to take up water from the super-saturated
1091!--          air.
1092             IF ( humidity )  THEN
1093                IF ( q_s .LE. q_p(k+1,j,i) )  THEN
1094                   IF ( .NOT. dewfall )  THEN
1095                      f_qsws_veg  = 0.0_wp
1096                      f_qsws_soil = 0.0_wp
1097                      f_qsws_liq  = 0.0_wp
1098                   ENDIF
1099                ENDIF
1100             ENDIF
1101
1102             f_shf  = rho_cp / r_a(j,i)
1103             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1104
1105!
1106!--          Calculate derivative of q_s for Taylor series expansion
1107             e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
1108                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
1109                              / (t_surface(j,i) - 35.86_wp)**2 )
1110
1111             dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
1112
1113             T_1 = pt_p(k+1,j,i) * exn
1114
1115!
1116!--          Add LW up so that it can be removed in prognostic equation
1117             rad_net_l(j,i) = rad_net(j,i) + sigma_sb * t_surface(j,i) ** 4
1118
1119             IF ( humidity )  THEN
1120
1121!
1122!--             Numerator of the prognostic equation
1123                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
1124                         + f_shf / exn * T_1 + f_qsws * ( q_p(k+1,j,i) - q_s   &
1125                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
1126                         * t_soil(nzb_soil,j,i)
1127
1128!
1129!--             Denominator of the prognostic equation
1130                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
1131                         * dq_s_dt + lambda_surface + f_shf / exn
1132
1133             ELSE
1134
1135!
1136!--             Numerator of the prognostic equation
1137                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
1138                         + f_shf / exn * T_1 + lambda_surface                  &
1139                         * t_soil(nzb_soil,j,i)
1140
1141!
1142!--             Denominator of the prognostic equation
1143                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
1144                         + lambda_surface + f_shf / exn
1145
1146             ENDIF
1147
1148             tend = 0.0_wp
1149
1150!
1151!--          Implicit solution when the surface layer has no heat capacity,
1152!--          otherwise use RK3 scheme.
1153             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface *        &
1154                                t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d&
1155                                * tsc(2) ) 
1156!
1157!--          Add RK3 term
1158             t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)              &
1159                                * tt_surface_m(j,i)
1160
1161!
1162!--          Calculate true tendency
1163             tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)        &
1164                    * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
1165!
1166!--          Calculate t_surface tendencies for the next Runge-Kutta step
1167             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1168                IF ( intermediate_timestep_count == 1 )  THEN
1169                   tt_surface_m(j,i) = tend
1170                ELSEIF ( intermediate_timestep_count <                         &
1171                         intermediate_timestep_count_max )  THEN
1172                   tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp           &
1173                                       * tt_surface_m(j,i)
1174                ENDIF
1175             ENDIF
1176
1177             pt_p(k,j,i) = t_surface_p(j,i) / exn
1178!
1179!--          Calculate fluxes
1180             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb                 &
1181                              * t_surface(j,i)**4 - 4.0_wp * sigma_sb          &
1182                              * t_surface(j,i)**3 * t_surface_p(j,i)
1183             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
1184                              - t_soil(nzb_soil,j,i))
1185             shf_eb(j,i)    = - f_shf  * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
1186
1187             IF ( humidity )  THEN
1188                qsws_eb(j,i)  = - f_qsws    * ( q_p(k+1,j,i) - q_s + dq_s_dt   &
1189                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
1190
1191                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( q_p(k+1,j,i) - q_s       &
1192                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1193                                    * t_surface_p(j,i) )
1194
1195                qsws_soil_eb(j,i) = - f_qsws_soil * ( q_p(k+1,j,i) - q_s       &
1196                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1197                                    * t_surface_p(j,i) )
1198
1199                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( q_p(k+1,j,i) - q_s       &
1200                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1201                                    * t_surface_p(j,i) )
1202             ENDIF
1203
1204!
1205!--          Calculate the true surface resistance
1206             IF ( qsws_eb(j,i) .EQ. 0.0_wp )  THEN
1207                r_s(j,i) = 1.0E10_wp
1208             ELSE
1209                r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dt           &
1210                           * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
1211                           / qsws_eb(j,i) - r_a(j,i)
1212             ENDIF
1213
1214!
1215!--          Calculate change in liquid water reservoir due to dew fall or
1216!--          evaporation of liquid water
1217             IF ( humidity )  THEN
1218!
1219!--             If precipitation is activated, add rain water to qsws_liq_eb
1220!--             and qsws_soil_eb according the the vegetation coverage.
1221!--             precipitation_rate is given in mm.
1222                IF ( precipitation )  THEN
1223
1224!
1225!--                Add precipitation to liquid water reservoir, if possible.
1226!--                Otherwise, add the water to bare soil fraction.
1227                   IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max )  THEN
1228                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
1229                                       + c_veg(j,i) * precipitation_rate(j,i)  &
1230                                       * 0.001_wp * rho_l * l_v
1231                   ELSE
1232                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1233                                        + c_veg(j,i) * precipitation_rate(j,i) &
1234                                        * 0.001_wp * rho_l * l_v
1235                   ENDIF
1236
1237!--                Add precipitation to bare soil according to the bare soil
1238!--                coverage.
1239                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp             &
1240                                       - c_veg(j,i)) * precipitation_rate(j,i) &
1241                                       * 0.001_wp * rho_l * l_v
1242                ENDIF
1243!
1244!--             If the air is saturated, check the reservoir water level
1245                IF ( q_s .LE. q_p(k+1,j,i))  THEN
1246!
1247!--                Check if reservoir is full (avoid values > m_liq_eb_max)
1248!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this
1249!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
1250!--                so that tend is zero and no further check is needed
1251                   IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max )  THEN
1252                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1253                                           + qsws_liq_eb(j,i)
1254                      qsws_liq_eb(j,i)  = 0.0_wp
1255                   ENDIF
1256
1257!
1258!--                In case qsws_veg_eb becomes negative (unphysical behavior),
1259!--                let the water enter the liquid water reservoir as dew on the
1260!--                plant
1261                   IF ( qsws_veg_eb(j,i) .LT. 0.0_wp )  THEN
1262                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
1263                      qsws_veg_eb(j,i) = 0.0_wp
1264                   ENDIF
1265                ENDIF                   
1266 
1267                tend = - qsws_liq_eb(j,i) * drho_l_lv
1268
1269                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
1270                                                   + tsc(3) * tm_liq_eb_m(j,i) )
1271
1272!
1273!--             Check if reservoir is overfull -> reduce to maximum
1274!--             (conservation of water is violated here)
1275                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
1276
1277!
1278!--             Check if reservoir is empty (avoid values < 0.0)
1279!--             (conservation of water is violated here)
1280                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
1281
1282
1283!
1284!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
1285                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1286                   IF ( intermediate_timestep_count == 1 )  THEN
1287                      tm_liq_eb_m(j,i) = tend
1288                   ELSEIF ( intermediate_timestep_count <                      &
1289                            intermediate_timestep_count_max )  THEN
1290                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
1291                                      * tm_liq_eb_m(j,i)
1292                   ENDIF
1293                ENDIF
1294
1295             ENDIF
1296
1297          ENDDO
1298       ENDDO
1299
1300    END SUBROUTINE lsm_energy_balance
1301
1302
1303!------------------------------------------------------------------------------!
1304! Description:
1305! ------------
1306!> Soil model as part of the land surface model. The model predicts soil
1307!> temperature and water content.
1308!------------------------------------------------------------------------------!
1309    SUBROUTINE lsm_soil_model
1310
1311
1312       IMPLICIT NONE
1313
1314       INTEGER(iwp) ::  i   !< running index
1315       INTEGER(iwp) ::  j   !< running index
1316       INTEGER(iwp) ::  k   !< running index
1317
1318       REAL(wp)     :: h_vg !< Van Genuchten coef. h
1319
1320       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
1321                                               lambda_temp, & !< temp. lambda
1322                                               tend           !< tendency
1323
1324       DO i = nxlg, nxrg   
1325          DO j = nysg, nyng
1326             DO k = nzb_soil, nzt_soil
1327!
1328!--             Calculate volumetric heat capacity of the soil, taking into
1329!--             account water content
1330                rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i))       &
1331                                     + rho_c_water * m_soil(k,j,i))
1332
1333!
1334!--             Calculate soil heat conductivity at the center of the soil
1335!--             layers
1336                ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
1337                lambda_temp(k) = ke * (lambda_h_sat(j,i) + lambda_h_dry) +     &
1338                                 lambda_h_dry
1339
1340             ENDDO
1341
1342!
1343!--          Calculate soil heat conductivity (lambda_h) at the _stag level
1344!--          using linear interpolation
1345             DO k = nzb_soil, nzt_soil-1
1346                 
1347                lambda_h(k,j,i) = lambda_temp(k) +                             &
1348                                  ( lambda_temp(k+1) - lambda_temp(k) )        &
1349                                  * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
1350
1351             ENDDO
1352             lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
1353
1354!
1355!--          Prognostic equation for soil temperature t_soil
1356             tend(:) = 0.0_wp
1357             tend(0) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *                    &
1358                       ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)     &
1359                         - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil)         &
1360                         + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
1361
1362             DO  k = 1, nzt_soil
1363                tend(k) = (1.0_wp/rho_c_total(k,j,i))                          &
1364                          * (   lambda_h(k,j,i)                                &
1365                              * ( t_soil(k+1,j,i) - t_soil(k,j,i) )            &
1366                              * ddz_soil(k)                                    &
1367                              - lambda_h(k-1,j,i)                              &
1368                              * ( t_soil(k,j,i) - t_soil(k-1,j,i) )            &
1369                              * ddz_soil(k-1)                                  &
1370                            ) * ddz_soil_stag(k)
1371             ENDDO
1372
1373             t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)   &
1374                                             + dt_3d * ( tsc(2)                &
1375                                             * tend(:) + tsc(3)                &
1376                                             * tt_soil_m(:,j,i) )   
1377
1378!
1379!--          Calculate t_soil tendencies for the next Runge-Kutta step
1380             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1381                IF ( intermediate_timestep_count == 1 )  THEN
1382                   DO  k = nzb_soil, nzt_soil
1383                      tt_soil_m(k,j,i) = tend(k)
1384                   ENDDO
1385                ELSEIF ( intermediate_timestep_count <                         &
1386                         intermediate_timestep_count_max )  THEN
1387                   DO  k = nzb_soil, nzt_soil
1388                      tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
1389                                         * tt_soil_m(k,j,i)
1390                   ENDDO
1391                ENDIF
1392             ENDIF
1393
1394
1395             DO k = nzb_soil, nzt_soil
1396
1397!
1398!--             Calculate soil diffusivity at the center of the soil layers
1399                lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat          &
1400                                  / m_sat(j,i) ) * ( MAX(m_soil(k,j,i),        &
1401                                  m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp)
1402
1403!
1404!--             Parametrization of Van Genuchten
1405                IF ( soil_type /= 7 )  THEN
1406!
1407!--                Calculate the hydraulic conductivity after Van Genuchten
1408!--                (1980)
1409                   h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -       &
1410                              MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i)   &
1411                              / (n_vg(j,i)-1.0_wp)) - 1.0_wp                   &
1412                          )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i)
1413
1414                   gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +            &
1415                                   (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp   &
1416                                   -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg    &
1417                                   )**(n_vg(j,i)-1.0_wp))**2 )                 &
1418                                   / ( (1.0_wp + (alpha_vg(j,i)*h_vg           &
1419                                   )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) &
1420                                   *(l_vg(j,i) + 2.0_wp)) )
1421
1422!
1423!--             Parametrization of Clapp & Hornberger
1424                ELSE
1425                   gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i)           &
1426                                   / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
1427                ENDIF
1428
1429             ENDDO
1430
1431
1432             IF ( humidity )  THEN
1433!
1434!--             Calculate soil diffusivity (lambda_w) at the _stag level
1435!--             using linear interpolation. To do: replace this with
1436!--             ECMWF-IFS Eq. 8.81
1437                DO k = nzb_soil, nzt_soil-1
1438                     
1439                   lambda_w(k,j,i) = lambda_temp(k) +                          &
1440                                     ( lambda_temp(k+1) - lambda_temp(k) )     &
1441                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)
1442                   gamma_w(k,j,i)  = gamma_temp(k) +                           &
1443                                     ( gamma_temp(k+1) - gamma_temp(k) )       &
1444                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)                 
1445
1446                ENDDO
1447
1448!
1449!
1450!--             In case of a closed bottom (= water content is conserved), set
1451!--             hydraulic conductivity to zero to that no water will be lost
1452!--             in the bottom layer.
1453                IF ( conserve_water_content )  THEN
1454                   gamma_w(nzt_soil,j,i) = 0.0_wp
1455                ELSE
1456                   gamma_w(nzt_soil,j,i) = lambda_temp(nzt_soil)
1457                ENDIF     
1458
1459!--             The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
1460!--             ensures the mass conservation for water. The transpiration of
1461!--             plants equals the cumulative withdrawals by the roots in the
1462!--             soil. The scheme takes into account the availability of water
1463!--             in the soil layers as well as the root fraction in the
1464!--             respective layer. Layer with moisture below wilting point will
1465!--             not contribute, which reflects the preference of plants to
1466!--             take water from moister layers.
1467
1468!
1469!--             Calculate the root extraction (ECMWF 7.69, modified so that the
1470!--             sum of root_extr = 1). The energy balance solver guarantees a
1471!--             positive transpiration, so that there is no need for an
1472!--             additional check.
1473
1474                m_total = 0.0_wp
1475                DO k = nzb_soil, nzt_soil
1476                    IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
1477                       m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
1478                    ENDIF
1479                ENDDO 
1480
1481                DO k = nzb_soil, nzt_soil 
1482                   IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
1483                      root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
1484                   ELSE
1485                      root_extr(k) = 0.0_wp
1486                   ENDIF
1487                ENDDO
1488
1489!
1490!--             Prognostic equation for soil water content m_soil
1491                tend(:) = 0.0_wp
1492                tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (                  &
1493                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
1494                            * ddz_soil(nzb_soil) - gamma_w(nzb_soil,j,i) - (   &
1495                            root_extr(nzb_soil) * qsws_veg_eb(j,i)             &
1496                            + qsws_soil_eb(j,i) ) * drho_l_lv )                &
1497                            * ddz_soil_stag(nzb_soil)
1498
1499                DO  k = nzb_soil+1, nzt_soil-1
1500                   tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)             &
1501                               - m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)&
1502                               - lambda_w(k-1,j,i) * (m_soil(k,j,i) -          &
1503                               m_soil(k-1,j,i)) * ddz_soil(k-1)                &
1504                               + gamma_w(k-1,j,i) - (root_extr(k)              &
1505                               * qsws_veg_eb(j,i) * drho_l_lv)                 &
1506                             ) * ddz_soil_stag(k)
1507
1508                ENDDO
1509                tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                     &
1510                                        - lambda_w(nzt_soil-1,j,i)             &
1511                                        * (m_soil(nzt_soil,j,i)                &
1512                                        - m_soil(nzt_soil-1,j,i))              &
1513                                        * ddz_soil(nzt_soil-1)                 &
1514                                        + gamma_w(nzt_soil-1,j,i) - (          &
1515                                          root_extr(nzt_soil)                  &
1516                                        * qsws_veg_eb(j,i) * drho_l_lv  )      &
1517                                      ) * ddz_soil_stag(nzt_soil)             
1518
1519                m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
1520                                                + dt_3d * ( tsc(2) * tend(:)   &
1521                                                + tsc(3) * tm_soil_m(:,j,i) )   
1522
1523!
1524!--             Account for dry soils (find a better solution here!)
1525                m_soil_p(:,j,i) = MAX(m_soil_p(:,j,i),0.0_wp)
1526
1527!
1528!--             Calculate m_soil tendencies for the next Runge-Kutta step
1529                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1530                   IF ( intermediate_timestep_count == 1 )  THEN
1531                      DO  k = nzb_soil, nzt_soil
1532                         tm_soil_m(k,j,i) = tend(k)
1533                      ENDDO
1534                   ELSEIF ( intermediate_timestep_count <                      &
1535                            intermediate_timestep_count_max )  THEN
1536                      DO  k = nzb_soil, nzt_soil
1537                         tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
1538                                     * tm_soil_m(k,j,i)
1539                      ENDDO
1540                   ENDIF
1541                ENDIF
1542
1543             ENDIF
1544
1545          ENDDO
1546       ENDDO
1547
1548!
1549!--    Calculate surface specific humidity
1550       IF ( humidity )  THEN
1551          CALL calc_q_surface
1552       ENDIF
1553
1554
1555    END SUBROUTINE lsm_soil_model
1556
1557
1558!------------------------------------------------------------------------------!
1559! Description:
1560! ------------
1561!> Calculation of specific humidity of the surface layer (surface)
1562!------------------------------------------------------------------------------!
1563    SUBROUTINE calc_q_surface
1564
1565       IMPLICIT NONE
1566
1567       INTEGER :: i              !< running index
1568       INTEGER :: j              !< running index
1569       INTEGER :: k              !< running index
1570       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
1571
1572       DO i = nxlg, nxrg   
1573          DO j = nysg, nyng
1574             k = nzb_s_inner(j,i)
1575
1576!
1577!--          Calculate water vapour pressure at saturation
1578             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
1579                                              - 273.16_wp ) / ( t_surface(j,i) &
1580                                              - 35.86_wp ) )
1581
1582!
1583!--          Calculate specific humidity at saturation
1584             q_s = 0.622_wp * e_s / surface_pressure
1585
1586
1587             resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
1588
1589!
1590!--          Calculate specific humidity at surface
1591             q_p(k,j,i) = resistance * q_s + (1.0_wp - resistance)             &
1592                          * q_p(k+1,j,i)
1593
1594          ENDDO
1595       ENDDO
1596
1597    END SUBROUTINE calc_q_surface
1598
1599!------------------------------------------------------------------------------!
1600! Description:
1601! ------------
1602!> Swapping of timelevels
1603!------------------------------------------------------------------------------!
1604    SUBROUTINE lsm_swap_timelevel ( mod_count )
1605
1606       IMPLICIT NONE
1607
1608       INTEGER, INTENT(IN) :: mod_count
1609
1610#if defined( __nopointer )
1611
1612       t_surface    = t_surface_p
1613       t_soil       = t_soil_p
1614       IF ( humidity )  THEN
1615          m_soil    = m_soil_p
1616          m_liq_eb  = m_liq_eb_p
1617       ENDIF
1618
1619#else
1620   
1621       SELECT CASE ( mod_count )
1622
1623          CASE ( 0 )
1624
1625             t_surface  => t_surface_1; t_surface_p  => t_surface_2
1626             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
1627             IF ( humidity )  THEN
1628                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
1629                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
1630             ENDIF
1631
1632
1633          CASE ( 1 )
1634
1635             t_surface  => t_surface_2; t_surface_p  => t_surface_1
1636             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
1637             IF ( humidity )  THEN
1638                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
1639                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
1640             ENDIF
1641
1642       END SELECT
1643#endif
1644
1645    END SUBROUTINE lsm_swap_timelevel
1646
1647
1648 END MODULE land_surface_model_mod
Note: See TracBrowser for help on using the repository browser.