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

Last change on this file since 1590 was 1590, checked in by maronga, 7 years ago

bugfix in definition of character arrays within the LSM and radiation code

  • 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! Bugfix: definition of character strings requires same length for all elements
23!
24! Former revisions:
25! -----------------
26! $Id: land_surface_model.f90 1590 2015-05-08 13:56:27Z maronga $
27! Bugfix: REAL constants provided with KIND-attribute in call of
28! 1585 2015-04-30 07:05:52Z maronga
29! Modifications for RRTMG. Changed tables to PARAMETER type.
30!
31! 1571 2015-03-12 16:12:49Z maronga
32! Removed upper-case variable names. Corrected distribution of precipitation to
33! the liquid water reservoir and the bare soil fractions.
34!
35! 1555 2015-03-04 17:44:27Z maronga
36! Added output of r_a and r_s
37!
38! 1553 2015-03-03 17:33:54Z maronga
39! Improved better treatment of roughness lengths. Added default soil temperature
40! profile
41!
42! 1551 2015-03-03 14:18:16Z maronga
43! Flux calculation is now done in prandtl_fluxes. Added support for data output.
44! Vertical indices have been replaced. Restart runs are now possible. Some
45! variables have beem renamed. Bugfix in the prognostic equation for the surface
46! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
47! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
48! the hydraulic conductivity. Bugfix for root fraction and extraction
49! calculation
50!
51! intrinsic function MAX and MIN
52!
53! 1500 2014-12-03 17:42:41Z maronga
54! Corrected calculation of aerodynamic resistance (r_a).
55! Precipitation is now added to liquid water reservoir using LE_liq.
56! Added support for dry runs.
57!
58! 1496 2014-12-02 17:25:50Z maronga
59! Initial revision
60!
61!
62! Description:
63! ------------
64! Land surface model, consisting of a solver for the energy balance at the
65! surface and a four layer soil scheme. The scheme is similar to the TESSEL
66! scheme implemented in the ECMWF IFS model, with modifications according to
67! H-TESSEL. The implementation is based on the formulation implemented in the
68! DALES and UCLA-LES models.
69!
70! To do list:
71! -----------
72! - Check dewfall parametrization for fog simulations
73! - Consider partial absorption of the net shortwave radiation by the surface layer
74! - Allow for water surfaces, check performance for bare soils
75! - Invert indices (running from -3 to 0. Currently: nzb_soil=0, nzt_soil=3))
76! - Implement surface runoff model (required when performing long-term LES with
77!   considerable precipitation
78!
79! Notes:
80! ------
81! - No time step criterion is required as long as the soil layers do not become
82!   too thin
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!
921! Note: surface fluxes are calculated in the land surface model, but these are
922! not used in the atmospheric code. The fluxes are calculated afterwards in
923! prandtl_fluxes using the surface values of temperature and humidity as
924! provided by the land surface model. In this way, the fluxes in the land
925! surface model are not equal to the ones calculated in prandtl_fluxes
926!------------------------------------------------------------------------------!
927    SUBROUTINE lsm_energy_balance
928
929
930       IMPLICIT NONE
931
932       INTEGER(iwp) ::  i         !: running index
933       INTEGER(iwp) ::  j         !: running index
934       INTEGER(iwp) ::  k, ks     !: running index
935
936       REAL(wp) :: f1,          & !: resistance correction term 1
937                   f2,          & !: resistance correction term 2
938                   f3,          & !: resistance correction term 3
939                   m_min,       & !: minimum soil moisture
940                   T_1,         & !: actual temperature at first grid point
941                   e,           & !: water vapour pressure
942                   e_s,         & !: water vapour saturation pressure
943                   e_s_dt,      & !: derivate of e_s with respect to T
944                   tend,        & !: tendency
945                   dq_s_dt,     & !: derivate of q_s with respect to T
946                   coef_1,      & !: coef. for prognostic equation
947                   coef_2,      & !: coef. for prognostic equation
948                   f_qsws,      & !: factor for qsws_eb
949                   f_qsws_veg,  & !: factor for qsws_veg_eb
950                   f_qsws_soil, & !: factor for qsws_soil_eb
951                   f_qsws_liq,  & !: factor for qsws_liq_eb
952                   f_shf,       & !: factor for shf_eb
953                   lambda_surface, & !: Current value of lambda_surface
954                   m_liq_eb_max   !: maxmimum value of the liq. water reservoir
955
956
957!
958!--    Calculate the exner function for the current time step
959       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
960
961       
962
963       DO i = nxlg, nxrg
964          DO j = nysg, nyng
965             k = nzb_s_inner(j,i)
966
967!
968!--          Set lambda_surface according to stratification
969             IF ( rif(j,i) >= 0.0_wp )  THEN
970                lambda_surface = lambda_surface_s(j,i)
971             ELSE
972                lambda_surface = lambda_surface_u(j,i)
973             ENDIF
974
975!
976!--          First step: calculate aerodyamic resistance. As pt, us, ts
977!--          are not available for the prognostic time step, data from the last
978!--          time step is used here. Note that this formulation is the
979!--          equivalent to the ECMWF formulation using drag coefficients
980             r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) +       &
981                         1.0E-20_wp)
982
983!
984!--          Second step: calculate canopy resistance r_canopy
985!--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
986 
987!--          f1: correction for incoming shortwave radiation (stomata close at
988!--          night)
989             IF ( radiation_scheme /= 'constant' )  THEN
990                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /&
991                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)        &
992                               + 1.0_wp)) )
993             ELSE
994                f1 = 1.0_wp
995             ENDIF
996
997!
998!--          f2: correction for soil moisture availability to plants (the
999!--          integrated soil moisture must thus be considered here)
1000!--          f2 = 0 for very dry soils
1001             m_total = 0.0_wp
1002             DO ks = nzb_soil, nzt_soil
1003                 m_total = m_total + root_fr(ks,j,i)                           &
1004                           * MAX(m_soil(ks,j,i),m_wilt(j,i))
1005             ENDDO
1006
1007             IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LT. m_fc(j,i) )  THEN
1008                f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
1009             ELSEIF ( m_total .GE. m_fc(j,i) )  THEN
1010                f2 = 1.0_wp
1011             ELSE
1012                f2 = 1.0E-20_wp
1013             ENDIF
1014
1015!
1016!--          Calculate water vapour pressure at saturation
1017             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
1018                           - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
1019
1020!
1021!--          f3: correction for vapour pressure deficit
1022             IF ( g_d(j,i) /= 0.0_wp )  THEN
1023!
1024!--             Calculate vapour pressure
1025                e  = q_p(k+1,j,i) * surface_pressure / 0.622_wp
1026                f3 = EXP ( -g_d(j,i) * (e_s - e) )
1027             ELSE
1028                f3 = 1.0_wp
1029             ENDIF
1030
1031!
1032!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
1033!--          this calculation is obsolete, as r_canopy is not used below.
1034!--          To do: check for very dry soil -> r_canopy goes to infinity
1035             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
1036                                             + 1.0E-20_wp)
1037
1038!
1039!--          Third step: calculate bare soil resistance r_soil. The Clapp &
1040!--          Hornberger parametrization does not consider c_veg.
1041             IF ( soil_type /= 7 )  THEN
1042                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
1043                        m_res(j,i)
1044             ELSE
1045                m_min = m_wilt(j,i)
1046             ENDIF
1047
1048             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
1049             f2 = MAX(f2,1.0E-20_wp)
1050             f2 = MIN(f2,1.0_wp)
1051
1052             r_soil(j,i) = r_soil_min(j,i) / f2
1053
1054!
1055!--          Calculate fraction of liquid water reservoir
1056             m_liq_eb_max = m_max_depth * lai(j,i)
1057             c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp))
1058             
1059
1060!
1061!--          Calculate saturation specific humidity
1062             q_s = 0.622_wp * e_s / surface_pressure
1063
1064!
1065!--          In case of dewfall, set evapotranspiration to zero
1066!--          All super-saturated water is then removed from the air
1067             IF ( humidity .AND. dewfall .AND. q_s .LE. q_p(k+1,j,i) )  THEN
1068                r_canopy(j,i) = 0.0_wp
1069                r_soil(j,i)   = 0.0_wp
1070             ENDIF
1071
1072!
1073!--          Calculate coefficients for the total evapotranspiration
1074             f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
1075                           (r_a(j,i) + r_canopy(j,i))
1076
1077             f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
1078                                                          r_soil(j,i))
1079             f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
1080
1081
1082!
1083!--          If soil moisture is below wilting point, plants do no longer
1084!--          transpirate.
1085!              IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
1086!                 f_qsws_veg = 0.0_wp
1087!              ENDIF
1088
1089!
1090!--          If dewfall is deactivated, vegetation, soil and liquid water
1091!--          reservoir are not allowed to take up water from the super-saturated
1092!--          air.
1093             IF ( humidity )  THEN
1094                IF ( q_s .LE. q_p(k+1,j,i) )  THEN
1095                   IF ( .NOT. dewfall )  THEN
1096                      f_qsws_veg  = 0.0_wp
1097                      f_qsws_soil = 0.0_wp
1098                      f_qsws_liq  = 0.0_wp
1099                   ENDIF
1100                ENDIF
1101             ENDIF
1102
1103             f_shf  = rho_cp / r_a(j,i)
1104             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
1105
1106!
1107!--          Calculate derivative of q_s for Taylor series expansion
1108             e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
1109                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
1110                              / (t_surface(j,i) - 35.86_wp)**2 )
1111
1112             dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
1113
1114             T_1 = pt_p(k+1,j,i) * exn
1115
1116!
1117!--          Add LW up so that it can be removed in prognostic equation
1118             rad_net_l(j,i) = rad_net(j,i) + sigma_sb * t_surface(j,i) ** 4
1119
1120             IF ( humidity )  THEN
1121
1122!
1123!--             Numerator of the prognostic equation
1124                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
1125                         + f_shf / exn * T_1 + f_qsws * ( q_p(k+1,j,i) - q_s   &
1126                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
1127                         * t_soil(nzb_soil,j,i)
1128
1129!
1130!--             Denominator of the prognostic equation
1131                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
1132                         * dq_s_dt + lambda_surface + f_shf / exn
1133
1134             ELSE
1135
1136!
1137!--             Numerator of the prognostic equation
1138                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
1139                         + f_shf / exn * T_1 + lambda_surface                  &
1140                         * t_soil(nzb_soil,j,i)
1141
1142!
1143!--             Denominator of the prognostic equation
1144                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
1145                         + lambda_surface + f_shf / exn
1146
1147             ENDIF
1148
1149             tend = 0.0_wp
1150
1151!
1152!--          Implicit solution when the surface layer has no heat capacity,
1153!--          otherwise use RK3 scheme.
1154             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface *        &
1155                                t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d&
1156                                * tsc(2) ) 
1157!
1158!--          Add RK3 term
1159             t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)              &
1160                                * tt_surface_m(j,i)
1161
1162!
1163!--          Calculate true tendency
1164             tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)        &
1165                    * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
1166!
1167!--          Calculate t_surface tendencies for the next Runge-Kutta step
1168             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1169                IF ( intermediate_timestep_count == 1 )  THEN
1170                   tt_surface_m(j,i) = tend
1171                ELSEIF ( intermediate_timestep_count <                         &
1172                         intermediate_timestep_count_max )  THEN
1173                   tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp           &
1174                                       * tt_surface_m(j,i)
1175                ENDIF
1176             ENDIF
1177
1178             pt_p(k,j,i) = t_surface_p(j,i) / exn
1179!
1180!--          Calculate fluxes
1181             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb                 &
1182                              * t_surface(j,i)**4 - 4.0_wp * sigma_sb          &
1183                              * t_surface(j,i)**3 * t_surface_p(j,i)
1184             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
1185                              - t_soil(nzb_soil,j,i))
1186             shf_eb(j,i)    = - f_shf  * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
1187
1188             IF ( humidity )  THEN
1189                qsws_eb(j,i)  = - f_qsws    * ( q_p(k+1,j,i) - q_s + dq_s_dt   &
1190                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
1191
1192                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( q_p(k+1,j,i) - q_s       &
1193                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1194                                    * t_surface_p(j,i) )
1195
1196                qsws_soil_eb(j,i) = - f_qsws_soil * ( q_p(k+1,j,i) - q_s       &
1197                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1198                                    * t_surface_p(j,i) )
1199
1200                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( q_p(k+1,j,i) - q_s       &
1201                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
1202                                    * t_surface_p(j,i) )
1203             ENDIF
1204
1205!
1206!--          Calculate the true surface resistance
1207             IF ( qsws_eb(j,i) .EQ. 0.0_wp )  THEN
1208                r_s(j,i) = 1.0E10_wp
1209             ELSE
1210                r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dt           &
1211                           * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
1212                           / qsws_eb(j,i) - r_a(j,i)
1213             ENDIF
1214
1215!
1216!--          Calculate change in liquid water reservoir due to dew fall or
1217!--          evaporation of liquid water
1218             IF ( humidity )  THEN
1219!
1220!--             If precipitation is activated, add rain water to qsws_liq_eb
1221!--             and qsws_soil_eb according the the vegetation coverage.
1222!--             precipitation_rate is given in mm.
1223                IF ( precipitation )  THEN
1224
1225!
1226!--                Add precipitation to liquid water reservoir, if possible.
1227!--                Otherwise, add the water to bare soil fraction.
1228                   IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max )  THEN
1229                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
1230                                       + c_veg(j,i) * precipitation_rate(j,i)  &
1231                                       * 0.001_wp * rho_l * l_v
1232                   ELSE
1233                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1234                                        + c_veg(j,i) * precipitation_rate(j,i) &
1235                                        * 0.001_wp * rho_l * l_v
1236                   ENDIF
1237
1238!--                Add precipitation to bare soil according to the bare soil
1239!--                coverage.
1240                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp             &
1241                                       - c_veg(j,i)) * precipitation_rate(j,i) &
1242                                       * 0.001_wp * rho_l * l_v
1243                ENDIF
1244!
1245!--             If the air is saturated, check the reservoir water level
1246                IF ( q_s .LE. q_p(k+1,j,i))  THEN
1247!
1248!--                Check if reservoir is full (avoid values > m_liq_eb_max)
1249!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this
1250!--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
1251!--                so that tend is zero and no further check is needed
1252                   IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max )  THEN
1253                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
1254                                           + qsws_liq_eb(j,i)
1255                      qsws_liq_eb(j,i)  = 0.0_wp
1256                   ENDIF
1257
1258!
1259!--                In case qsws_veg_eb becomes negative (unphysical behavior),
1260!--                let the water enter the liquid water reservoir as dew on the
1261!--                plant
1262                   IF ( qsws_veg_eb(j,i) .LT. 0.0_wp )  THEN
1263                      qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
1264                      qsws_veg_eb(j,i) = 0.0_wp
1265                   ENDIF
1266                ENDIF                   
1267 
1268                tend = - qsws_liq_eb(j,i) * drho_l_lv
1269
1270                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
1271                                                   + tsc(3) * tm_liq_eb_m(j,i) )
1272
1273!
1274!--             Check if reservoir is overfull -> reduce to maximum
1275!--             (conservation of water is violated here)
1276                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
1277
1278!
1279!--             Check if reservoir is empty (avoid values < 0.0)
1280!--             (conservation of water is violated here)
1281                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
1282
1283
1284!
1285!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
1286                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1287                   IF ( intermediate_timestep_count == 1 )  THEN
1288                      tm_liq_eb_m(j,i) = tend
1289                   ELSEIF ( intermediate_timestep_count <                      &
1290                            intermediate_timestep_count_max )  THEN
1291                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
1292                                      * tm_liq_eb_m(j,i)
1293                   ENDIF
1294                ENDIF
1295
1296             ENDIF
1297
1298          ENDDO
1299       ENDDO
1300
1301    END SUBROUTINE lsm_energy_balance
1302
1303
1304!------------------------------------------------------------------------------!
1305! Description:
1306! ------------
1307!
1308! Soil model as part of the land surface model. The model predicts soil
1309! temperature and water content.
1310!------------------------------------------------------------------------------!
1311    SUBROUTINE lsm_soil_model
1312
1313
1314       IMPLICIT NONE
1315
1316       INTEGER(iwp) ::  i   !: running index
1317       INTEGER(iwp) ::  j   !: running index
1318       INTEGER(iwp) ::  k   !: running index
1319
1320       REAL(wp)     :: h_vg !: Van Genuchten coef. h
1321
1322       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !: temp. gamma
1323                                               lambda_temp, & !: temp. lambda
1324                                               tend           !: tendency
1325
1326       DO i = nxlg, nxrg   
1327          DO j = nysg, nyng
1328             DO k = nzb_soil, nzt_soil
1329!
1330!--             Calculate volumetric heat capacity of the soil, taking into
1331!--             account water content
1332                rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i))       &
1333                                     + rho_c_water * m_soil(k,j,i))
1334
1335!
1336!--             Calculate soil heat conductivity at the center of the soil
1337!--             layers
1338                ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
1339                lambda_temp(k) = ke * (lambda_h_sat(j,i) + lambda_h_dry) +     &
1340                                 lambda_h_dry
1341
1342             ENDDO
1343
1344!
1345!--          Calculate soil heat conductivity (lambda_h) at the _stag level
1346!--          using linear interpolation
1347             DO k = nzb_soil, nzt_soil-1
1348                 
1349                lambda_h(k,j,i) = lambda_temp(k) +                             &
1350                                  ( lambda_temp(k+1) - lambda_temp(k) )        &
1351                                  * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
1352
1353             ENDDO
1354             lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
1355
1356!
1357!--          Prognostic equation for soil temperature t_soil
1358             tend(:) = 0.0_wp
1359             tend(0) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *                    &
1360                       ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)     &
1361                         - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil)         &
1362                         + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
1363
1364             DO  k = 1, nzt_soil
1365                tend(k) = (1.0_wp/rho_c_total(k,j,i))                          &
1366                          * (   lambda_h(k,j,i)                                &
1367                              * ( t_soil(k+1,j,i) - t_soil(k,j,i) )            &
1368                              * ddz_soil(k)                                    &
1369                              - lambda_h(k-1,j,i)                              &
1370                              * ( t_soil(k,j,i) - t_soil(k-1,j,i) )            &
1371                              * ddz_soil(k-1)                                  &
1372                            ) * ddz_soil_stag(k)
1373             ENDDO
1374
1375             t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)   &
1376                                             + dt_3d * ( tsc(2)                &
1377                                             * tend(:) + tsc(3)                &
1378                                             * tt_soil_m(:,j,i) )   
1379
1380!
1381!--          Calculate t_soil tendencies for the next Runge-Kutta step
1382             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1383                IF ( intermediate_timestep_count == 1 )  THEN
1384                   DO  k = nzb_soil, nzt_soil
1385                      tt_soil_m(k,j,i) = tend(k)
1386                   ENDDO
1387                ELSEIF ( intermediate_timestep_count <                         &
1388                         intermediate_timestep_count_max )  THEN
1389                   DO  k = nzb_soil, nzt_soil
1390                      tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
1391                                         * tt_soil_m(k,j,i)
1392                   ENDDO
1393                ENDIF
1394             ENDIF
1395
1396
1397             DO k = nzb_soil, nzt_soil
1398
1399!
1400!--             Calculate soil diffusivity at the center of the soil layers
1401                lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat          &
1402                                  / m_sat(j,i) ) * ( MAX(m_soil(k,j,i),        &
1403                                  m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp)
1404
1405!
1406!--             Parametrization of Van Genuchten
1407                IF ( soil_type /= 7 )  THEN
1408!
1409!--                Calculate the hydraulic conductivity after Van Genuchten
1410!--                (1980)
1411                   h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -       &
1412                              MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i)   &
1413                              / (n_vg(j,i)-1.0_wp)) - 1.0_wp                   &
1414                          )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i)
1415
1416                   gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +            &
1417                                   (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp   &
1418                                   -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg    &
1419                                   )**(n_vg(j,i)-1.0_wp))**2 )                 &
1420                                   / ( (1.0_wp + (alpha_vg(j,i)*h_vg           &
1421                                   )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) &
1422                                   *(l_vg(j,i) + 2.0_wp)) )
1423
1424!
1425!--             Parametrization of Clapp & Hornberger
1426                ELSE
1427                   gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i)           &
1428                                   / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
1429                ENDIF
1430
1431             ENDDO
1432
1433
1434             IF ( humidity )  THEN
1435!
1436!--             Calculate soil diffusivity (lambda_w) at the _stag level
1437!--             using linear interpolation. To do: replace this with
1438!--             ECMWF-IFS Eq. 8.81
1439                DO k = nzb_soil, nzt_soil-1
1440                     
1441                   lambda_w(k,j,i) = lambda_temp(k) +                          &
1442                                     ( lambda_temp(k+1) - lambda_temp(k) )     &
1443                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)
1444                   gamma_w(k,j,i)  = gamma_temp(k) +                           &
1445                                     ( gamma_temp(k+1) - gamma_temp(k) )       &
1446                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)                 
1447
1448                ENDDO
1449
1450!
1451!
1452!--             In case of a closed bottom (= water content is conserved), set
1453!--             hydraulic conductivity to zero to that no water will be lost
1454!--             in the bottom layer.
1455                IF ( conserve_water_content )  THEN
1456                   gamma_w(nzt_soil,j,i) = 0.0_wp
1457                ELSE
1458                   gamma_w(nzt_soil,j,i) = lambda_temp(nzt_soil)
1459                ENDIF     
1460
1461!--             The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
1462!--             ensures the mass conservation for water. The transpiration of
1463!--             plants equals the cumulative withdrawals by the roots in the
1464!--             soil. The scheme takes into account the availability of water
1465!--             in the soil layers as well as the root fraction in the
1466!--             respective layer. Layer with moisture below wilting point will
1467!--             not contribute, which reflects the preference of plants to
1468!--             take water from moister layers.
1469
1470!
1471!--             Calculate the root extraction (ECMWF 7.69, modified so that the
1472!--             sum of root_extr = 1). The energy balance solver guarantees a
1473!--             positive transpiration, so that there is no need for an
1474!--             additional check.
1475
1476                m_total = 0.0_wp
1477                DO k = nzb_soil, nzt_soil
1478                    IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
1479                       m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
1480                    ENDIF
1481                ENDDO 
1482
1483                DO k = nzb_soil, nzt_soil
1484                   IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
1485                      root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
1486                   ELSE
1487                      root_extr(k) = 0.0_wp
1488                   ENDIF
1489                ENDDO
1490
1491!
1492!--             Prognostic equation for soil water content m_soil
1493                tend(:) = 0.0_wp
1494                tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (                  &
1495                            m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
1496                            * ddz_soil(nzb_soil) - gamma_w(nzb_soil,j,i) - (   &
1497                            root_extr(nzb_soil) * qsws_veg_eb(j,i)             &
1498                            + qsws_soil_eb(j,i) ) * drho_l_lv )                &
1499                            * ddz_soil_stag(nzb_soil)
1500
1501                DO  k = nzb_soil+1, nzt_soil-1
1502                   tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)             &
1503                               - m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)&
1504                               - lambda_w(k-1,j,i) * (m_soil(k,j,i) -          &
1505                               m_soil(k-1,j,i)) * ddz_soil(k-1)                &
1506                               + gamma_w(k-1,j,i) - (root_extr(k)              &
1507                               * qsws_veg_eb(j,i) * drho_l_lv)                 &
1508                             ) * ddz_soil_stag(k)
1509
1510                ENDDO
1511                tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                     &
1512                                        - lambda_w(nzt_soil-1,j,i)             &
1513                                        * (m_soil(nzt_soil,j,i)                &
1514                                        - m_soil(nzt_soil-1,j,i))              &
1515                                        * ddz_soil(nzt_soil-1)                 &
1516                                        + gamma_w(nzt_soil-1,j,i) - (          &
1517                                          root_extr(nzt_soil)                  &
1518                                        * qsws_veg_eb(j,i) * drho_l_lv  )      &
1519                                      ) * ddz_soil_stag(nzt_soil)             
1520
1521                m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
1522                                                + dt_3d * ( tsc(2) * tend(:)   &
1523                                                + tsc(3) * tm_soil_m(:,j,i) )   
1524
1525!
1526!--             Account for dry soils (find a better solution here!)
1527                m_soil_p(:,j,i) = MAX(m_soil_p(:,j,i),0.0_wp)
1528
1529!
1530!--             Calculate m_soil tendencies for the next Runge-Kutta step
1531                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1532                   IF ( intermediate_timestep_count == 1 )  THEN
1533                      DO  k = nzb_soil, nzt_soil
1534                         tm_soil_m(k,j,i) = tend(k)
1535                      ENDDO
1536                   ELSEIF ( intermediate_timestep_count <                      &
1537                            intermediate_timestep_count_max )  THEN
1538                      DO  k = nzb_soil, nzt_soil
1539                         tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
1540                                     * tm_soil_m(k,j,i)
1541                      ENDDO
1542                   ENDIF
1543                ENDIF
1544
1545             ENDIF
1546
1547          ENDDO
1548       ENDDO
1549
1550!
1551!--    Calculate surface specific humidity
1552       IF ( humidity )  THEN
1553          CALL calc_q_surface
1554       ENDIF
1555
1556
1557    END SUBROUTINE lsm_soil_model
1558
1559
1560!------------------------------------------------------------------------------!
1561! Description:
1562! ------------
1563! Calculation of specific humidity of the surface layer (surface)
1564!------------------------------------------------------------------------------!
1565    SUBROUTINE calc_q_surface
1566
1567       IMPLICIT NONE
1568
1569       INTEGER :: i              !: running index
1570       INTEGER :: j              !: running index
1571       INTEGER :: k              !: running index
1572       REAL(wp) :: resistance    !: aerodynamic and soil resistance term
1573
1574       DO i = nxlg, nxrg   
1575          DO j = nysg, nyng
1576             k = nzb_s_inner(j,i)
1577
1578!
1579!--          Calculate water vapour pressure at saturation
1580             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
1581                                              - 273.16_wp ) / ( t_surface(j,i) &
1582                                              - 35.86_wp ) )
1583
1584!
1585!--          Calculate specific humidity at saturation
1586             q_s = 0.622_wp * e_s / surface_pressure
1587
1588
1589             resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
1590
1591!
1592!--          Calculate specific humidity at surface
1593             q_p(k,j,i) = resistance * q_s + (1.0_wp - resistance)             &
1594                          * q_p(k+1,j,i)
1595
1596          ENDDO
1597       ENDDO
1598
1599    END SUBROUTINE calc_q_surface
1600
1601!------------------------------------------------------------------------------!
1602! Description:
1603! ------------
1604! Swapping of timelevels
1605!------------------------------------------------------------------------------!
1606    SUBROUTINE lsm_swap_timelevel ( mod_count )
1607
1608       IMPLICIT NONE
1609
1610       INTEGER, INTENT(IN) :: mod_count
1611
1612#if defined( __nopointer )
1613
1614       t_surface    = t_surface_p
1615       t_soil       = t_soil_p
1616       IF ( humidity )  THEN
1617          m_soil    = m_soil_p
1618          m_liq_eb  = m_liq_eb_p
1619       ENDIF
1620
1621#else
1622   
1623       SELECT CASE ( mod_count )
1624
1625          CASE ( 0 )
1626
1627             t_surface  => t_surface_1; t_surface_p  => t_surface_2
1628             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
1629             IF ( humidity )  THEN
1630                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
1631                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
1632             ENDIF
1633
1634
1635          CASE ( 1 )
1636
1637             t_surface  => t_surface_2; t_surface_p  => t_surface_1
1638             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
1639             IF ( humidity )  THEN
1640                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
1641                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
1642             ENDIF
1643
1644       END SELECT
1645#endif
1646
1647    END SUBROUTINE lsm_swap_timelevel
1648
1649
1650 END MODULE land_surface_model_mod
Note: See TracBrowser for help on using the repository browser.