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

Last change on this file since 1612 was 1591, checked in by maronga, 9 years ago

last commit documented

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