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

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

minor bugfixes in lasm/radiation

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