Line data Source code
1 : !> \file mo_soil_moisture.f90
2 : !> \brief \copybrief mo_soil_moisture
3 : !> \details \copydetails mo_soil_moisture
4 :
5 : !> \brief Soil moisture of the different layers
6 : !> \details Soil moisture in the different layers is calculated with
7 : !! infiltration as \f$ (\theta / \theta_{sat})^\beta \f$
8 : !! Then evapotranspiration is calculated from PET with a soil water stress factor \f$ f_{SM} \f$
9 : !! either using the Feddes equation - processCase(3) = 1:
10 : !! \f[ f_{SM} = \frac{\theta - \theta_\mathit{pwp}}{\theta_\mathit{fc} - \theta_\mathit{pwp}} \f]
11 : !! or using the Jarvis equation - processCase(3) = 2:
12 : !! \f[ f_{SM} = \frac{1}{\theta_\mathit{stress-index-C1}}
13 : !! \frac{\theta - \theta_\mathit{pwp}}{\theta_\mathit{sat} - \theta_\mathit{pwp}} \f]
14 : !> \authors Matthias Cuntz, Luis Samaniego
15 : !> \date Dec 2012
16 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
17 : !! mHM is released under the LGPLv3+ license \license_note
18 : !> \ingroup f_mhm
19 : MODULE mo_soil_moisture
20 :
21 : USE mo_kind, ONLY : i4, dp
22 :
23 : IMPLICIT NONE
24 :
25 : PUBLIC :: feddes_et_reduction
26 : PUBLIC :: jarvis_et_reduction
27 :
28 : PUBLIC :: soil_moisture ! Soil moisture in different soil horizons
29 :
30 : ! ------------------------------------------------------------------
31 :
32 : CONTAINS
33 :
34 : ! ------------------------------------------------------------------
35 :
36 : ! NAME
37 : ! soil_moisture
38 :
39 : ! PURPOSE
40 : !> \brief Soil moisture in different soil horizons
41 :
42 : !> \details Infiltration \f$I\f$ from one layer \f$k-1\f$ to the next \f$k\f$ on
43 : !> pervious areas is calculated as (omit \f$t\f$)
44 : !> \f[ I[k] = I[k-1] (\theta[k] / \theta_{sat}[k])^{\beta[k]} \f]
45 : !> Then soil moisture can be calculated as (omit \f$k\f$)
46 : !> \f[ \theta[t] = \theta[t-1] + I[t] - \mathit{ET}[t] \f]
47 : !> with \f$ \mathit{ET} \f$ (omit \f$[k,t]\f$) being
48 : !> \f[ \mathit{ET} = f_\mathrm{roots} \cdot f_{SM} \cdot \mathit{PET} \f].
49 :
50 : ! INTENT(IN)
51 : !> \param[in] "integer(i4) :: processCase" 1 - Feddes equation for PET reduction2 - Jarvis
52 : !> equation for PET reduction3 - Jarvis equation for PET reduction and FC dependency on root fraction coefficient
53 : !> \param[in] "real(dp) :: frac_sealed" Fraction of sealed area
54 : !> \param[in] "real(dp) :: water_thresh_sealed" Threshhold water depth in impervious areas [mm TS-1]
55 : !> \param[in] "real(dp) :: pet" Reference evapotranspiration [mm TS-1]
56 : !> \param[in] "real(dp) :: evap_coeff" Evaporation coefficent for free-water surface of
57 : !> that current month
58 : !> \param[in] "real(dp), dimension(:) :: soil_moist_sat" Saturation soil moisture for each horizon [mm]
59 : !> \param[in] "real(dp), dimension(:) :: frac_roots" Fraction of Roots in soil horizon
60 : !> \param[in] "real(dp), dimension(:) :: soil_moist_FC" Soil moisture below which actual ET is reduced [mm]
61 : !> \param[in] "real(dp), dimension(:) :: wilting_point" Permanent wilting point for each horizon [mm]
62 : !> \param[in] "real(dp), dimension(:) :: soil_moist_exponen" Exponential parameter to how non-linear is the soil
63 : !> water retention
64 : !> \param[in] "real(dp) :: jarvis_thresh_c1" Jarvis critical value for normalized soil water
65 : !> content
66 : !> \param[in] "real(dp) :: aet_canopy" Actual ET from canopy [mm TS-1]
67 :
68 : ! INTENT(INOUT)
69 : !> \param[inout] "real(dp) :: prec_effec" Effective precipitation (rain +
70 : !> snow melt) [mm]
71 : !> \param[inout] "real(dp) :: runoff_sealed" Direct runoff from impervious
72 : !> areas
73 : !> \param[inout] "real(dp) :: storage_sealed" Retention storage of impervious
74 : !> areas
75 : !> \param[inout] "real(dp), dimension(size(soil_moist_sat, 1)) :: infiltration" Recharge, infiltration intensity
76 : !> oreffective precipitation of each horizon [mm TS-1]
77 : !> \param[inout] "real(dp), dimension(size(soil_moist_sat, 1)) :: soil_moist" Soil moisture of each horizon
78 : !> [mm]
79 :
80 : ! INTENT(OUT)
81 : !> \param[out] "real(dp), dimension(size(soil_moist_sat, 1)) :: aet" actual ET [mm TS-1]
82 : !> \param[out] "real(dp) :: aet_sealed" actual ET from free-water surfaces,i.e
83 : !> impervious cover [mm TS-1]
84 :
85 : ! HISTORY
86 : !> \authors Matthias Cuntz
87 :
88 : !> \date Dec 2012
89 :
90 : ! Modifications:
91 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
92 : ! M. Cuneyd Demirel, Simon Stisen Jun 2020 - added Feddes and FC dependency on root fraction coefficient processCase(3) = 4
93 :
94 76647048 : subroutine soil_moisture(processCase, frac_sealed, water_thresh_sealed, pet, evap_coeff, soil_moist_sat, frac_roots, &
95 76647048 : soil_moist_FC, wilting_point, soil_moist_exponen, jarvis_thresh_c1, aet_canopy, prec_effec, &
96 76647048 : runoff_sealed, storage_sealed, infiltration, soil_moist, aet, aet_sealed)
97 :
98 : use mo_common_constants, only : eps_dp
99 :
100 : implicit none
101 :
102 : ! 1 - Feddes equation for PET reduction2 - Jarvis equation for PET reduction3 - Jarvis equation for PET reduction
103 : ! and FC dependency on root fraction coefficient
104 : integer(i4), intent(in) :: processCase
105 :
106 : ! Fraction of sealed area
107 : real(dp), intent(in) :: frac_sealed
108 :
109 : ! Threshhold water depth in impervious areas [mm TS-1]
110 : real(dp), intent(in) :: water_thresh_sealed
111 :
112 : ! Reference evapotranspiration [mm TS-1]
113 : real(dp), intent(in) :: pet
114 :
115 : ! Evaporation coefficent for free-water surface of that current month
116 : real(dp), intent(in) :: evap_coeff
117 :
118 : ! Saturation soil moisture for each horizon [mm]
119 : real(dp), dimension(:), intent(in) :: soil_moist_sat
120 :
121 : ! Fraction of Roots in soil horizon
122 : real(dp), dimension(:), intent(in) :: frac_roots
123 :
124 : ! Soil moisture below which actual ET is reduced [mm]
125 : real(dp), dimension(:), intent(in) :: soil_moist_FC
126 :
127 : ! Permanent wilting point for each horizon [mm]
128 : real(dp), dimension(:), intent(in) :: wilting_point
129 :
130 : ! Exponential parameter to how non-linear is the soil water retention
131 : real(dp), dimension(:), intent(in) :: soil_moist_exponen
132 :
133 : ! Jarvis critical value for normalized soil water content
134 : real(dp), intent(in) :: jarvis_thresh_c1
135 :
136 : ! Actual ET from canopy [mm TS-1]
137 : real(dp), intent(in) :: aet_canopy
138 :
139 : ! Effective precipitation (rain + snow melt) [mm]
140 : real(dp), intent(in) :: prec_effec
141 :
142 : ! Direct runoff from impervious areas
143 : real(dp), intent(inout) :: runoff_sealed
144 :
145 : ! Retention storage of impervious areas
146 : real(dp), intent(inout) :: storage_sealed
147 :
148 : ! Recharge, infiltration intensity oreffective precipitation of each horizon [mm TS-1]
149 : real(dp), dimension(size(soil_moist_sat, 1)), intent(inout) :: infiltration
150 :
151 : ! Soil moisture of each horizon [mm]
152 : real(dp), dimension(size(soil_moist_sat, 1)), intent(inout) :: soil_moist
153 :
154 : ! actual ET [mm TS-1]
155 : real(dp), dimension(size(soil_moist_sat, 1)), intent(out) :: aet
156 :
157 : ! actual ET from free-water surfaces,i.e impervious cover [mm TS-1]
158 : real(dp), intent(out) :: aet_sealed
159 :
160 : ! counter
161 : integer(i4) :: hh
162 :
163 : ! Effective Prec or infiltration from above
164 76647048 : real(dp) :: prec_effec_soil
165 :
166 : ! Runoof fraction
167 76647048 : real(dp) :: frac_runoff
168 :
169 : ! PET reduction factor according to actual soil moisture
170 76647048 : real(dp) :: soil_stress_factor
171 :
172 : ! temporary variable for misc use
173 76647048 : real(dp) :: tmp
174 :
175 :
176 : ! ----------------------------------------------------------------
177 : ! IMPERVIOUS COVER PROCESS
178 : ! ----------------------------------------------------------------
179 76647048 : runoff_sealed = 0.0_dp
180 76647048 : aet_sealed = 0.0_dp
181 :
182 76647048 : if (frac_sealed > 0.0_dp) then
183 71921832 : tmp = storage_sealed + prec_effec
184 :
185 71921832 : if (tmp > water_thresh_sealed) then
186 17562477 : runoff_sealed = tmp - water_thresh_sealed
187 17562477 : storage_sealed = water_thresh_sealed
188 : else
189 : runoff_sealed = 0.0_dp
190 54359355 : storage_sealed = tmp
191 : end if
192 :
193 : ! aET from sealed area is propotional to the available water content
194 71921832 : if(water_thresh_sealed .gt. eps_dp) then
195 71921832 : aet_sealed = (pet / evap_coeff - aet_canopy) * (storage_sealed / water_thresh_sealed)
196 : ! numerical problem
197 71921832 : if (aet_sealed .lt. 0.0_dp) aet_sealed = 0.0_dp
198 : else
199 0 : aet_sealed = huge(1.0_dp)
200 : end if
201 :
202 : ! sealed storage updata
203 71921832 : if (storage_sealed .gt. aet_sealed) then
204 68320422 : storage_sealed = storage_sealed - aet_sealed
205 : else
206 3601410 : aet_sealed = storage_sealed
207 3601410 : storage_sealed = 0.0_dp
208 : end if
209 :
210 : end if
211 : ! ----------------------------------------------------------------
212 : ! N-LAYER SOIL MODULE
213 : ! ----------------------------------------------------------------
214 229941144 : aet(:) = 0.0_dp
215 229941144 : infiltration(:) = 0.0_dp
216 :
217 : ! for 1st layer input is prec_effec
218 76647048 : prec_effec_soil = prec_effec
219 :
220 229941144 : do hh = 1, size(soil_moist_sat, 1) ! nHorizons
221 : ! input for other layers is the infiltration from its immediate upper layer will be input
222 153294096 : if (hh .NE. 1) prec_effec_soil = infiltration(hh - 1)
223 :
224 : ! start processing for soil moisture process
225 : ! BASED ON SMs as its upper LIMIT
226 :
227 153294096 : if (soil_moist(hh) > soil_moist_sat(hh)) then
228 0 : infiltration(hh) = prec_effec_soil
229 : else
230 : ! to avoid underflow -- or numerical errors
231 153294096 : if(soil_moist(hh) > eps_dp) then
232 : !frac_runoff = (soil_moist(hh) / soil_moist_sat(hh))**soil_moist_exponen(hh)
233 153294096 : frac_runoff = exp(soil_moist_exponen(hh) * log(soil_moist(hh) / soil_moist_sat(hh)))
234 : else
235 : frac_runoff = 0.0_dp
236 : end if
237 153294096 : tmp = prec_effec_soil * (1.0_dp - frac_runoff)
238 :
239 153294096 : if ((soil_moist(hh) + tmp) > soil_moist_sat(hh)) then
240 0 : infiltration(hh) = prec_effec_soil + (soil_moist(hh) - soil_moist_sat(hh))
241 0 : soil_moist(hh) = soil_moist_sat(hh)
242 : else
243 153294096 : infiltration(hh) = prec_effec_soil - tmp
244 153294096 : soil_moist(hh) = soil_moist(hh) + tmp
245 : end if
246 : end if
247 :
248 : ! aET calculations
249 :
250 : ! Satisfying ET demand sequentially from top to the bottom layer
251 : ! Note that the potential ET for the first soil layer is reduced after
252 : ! satisfying ET demands of the canopy surface
253 :
254 153294096 : aet(hh) = pet - aet_canopy ! First layer
255 229941144 : if (hh /= 1) aet(hh) = aet(hh) - sum(aet(1 : hh - 1), mask = (aet(1 : hh - 1) > 0.0_dp)) ! remaining layers
256 :
257 : ! estimate fraction of ET demand based on root fraction and SM status
258 153294096 : select case(processCase)
259 : ! FEDDES EQUATION: https://doi.org/10.1016/0022-1694(76)90017-2
260 : case(1 , 4)
261 151511952 : soil_stress_factor = feddes_et_reduction(soil_moist(hh), soil_moist_FC(hh), wilting_point(hh), &
262 151511952 : frac_roots(hh))
263 : ! JARVIS EQUATION: https://doi.org/10.1016/0022-1694(89)90050-4
264 : case(2 , 3)
265 : !!!!!!!!! INTRODUCING STRESS FACTOR FOR SOIL MOISTURE ET REDUCTION !!!!!!!!!!!!!!!!!
266 1782144 : soil_stress_factor = jarvis_et_reduction(soil_moist(hh), soil_moist_sat(hh), wilting_point(hh), &
267 155076240 : frac_roots(hh), jarvis_thresh_c1)
268 : end select
269 :
270 153294096 : aet(hh) = aet(hh) * soil_stress_factor
271 :
272 : ! avoid numerical error
273 153294096 : if(aet(hh) < 0.0_dp) aet(hh) = 0.0_dp
274 :
275 : ! reduce SM state
276 153294096 : if(soil_moist(hh) > aet(hh)) then
277 153294096 : soil_moist(hh) = soil_moist(hh) - aet(hh)
278 : else
279 0 : aet(hh) = soil_moist(hh) - eps_dp
280 0 : soil_moist(hh) = eps_dp
281 : end if
282 :
283 : ! avoid numerical error of underflow
284 229941144 : if(soil_moist(hh) < eps_dp) soil_moist(hh) = eps_dp
285 :
286 : end do ! hh
287 :
288 :
289 :
290 76647048 : end subroutine soil_moisture
291 :
292 :
293 : ! ------------------------------------------------------------------
294 :
295 : ! NAME
296 : ! feddes_et_reduction
297 :
298 : ! PURPOSE
299 : !> \brief stress factor for reducing evapotranspiration based on actual soil moisture
300 :
301 : !> \details Potential evapotranspiration is reduced to 0 if SM is lower PWP. PET is equal
302 : !> fraction of roots if soil moisture is exceeding field capacity. If soil moisture is
303 : !> in between PWP and FC PET is reduced by fraction of roots times a stress factor.
304 :
305 : !> The ET reduction factor \f$ f \f$ is estimated as
306 : !> \f[ f = \left\{
307 : !> \begin{array}{lr}
308 : !> f_{roots} & if \theta \ge \theta_{fc}\\
309 : !> f_{roots} \cdot \frac{\theta - \theta_\mathit{pwp}}{\theta_\mathit{fc} - \theta_\mathit{pwp}} &
310 : !> if \theta < \theta_{fc} \\
311 : !> 0 & if \theta < \theta_{pwp}
312 : !> \end{array}
313 : !> \right. \f]
314 :
315 : ! INTENT(IN)
316 : !> \param[in] "real(dp) :: soil_moist" Soil moisture of each horizon [mm]
317 : !> \param[in] "real(dp) :: soil_moist_FC" Soil moisture below which actual ET is reduced [mm]
318 : !> \param[in] "real(dp) :: wilting_point" Permanent wilting point
319 : !> \param[in] "real(dp) :: frac_roots" Fraction of Roots in soil horizon is reduced [mm]
320 :
321 : ! RETURN
322 : !> \return real(dp) :: feddes_et_reduction; et reduction factor
323 :
324 : ! HISTORY
325 : !> \authors Matthias Cuntz, Cueneyd Demirel, Matthias Zink
326 :
327 : !> \date March 2017
328 :
329 : ! Modifications:
330 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
331 : ! M. Cuneyd Demirel, Simon Stisen Jun 2020 - added Feddes and FC dependency on root fraction coefficient processCase(3) = 4
332 :
333 151511952 : elemental pure FUNCTION feddes_et_reduction(soil_moist, soil_moist_FC, wilting_point, frac_roots)
334 : implicit none
335 :
336 : ! Soil moisture of each horizon [mm]
337 : real(dp), intent(in) :: soil_moist
338 :
339 : ! Soil moisture below which actual ET is reduced [mm]
340 : real(dp), intent(in) :: soil_moist_FC
341 :
342 : ! Permanent wilting point
343 : real(dp), intent(in) :: wilting_point
344 :
345 : ! Fraction of Roots in soil horizon is reduced [mm]
346 : real(dp), intent(in) :: frac_roots
347 :
348 : ! reference evapotranspiration in [mm s-1]
349 : real(dp) :: feddes_et_reduction
350 :
351 :
352 : ! SM >= FC
353 151511952 : if (soil_moist >= soil_moist_FC) then
354 82717277 : feddes_et_reduction = frac_roots
355 : ! PW < SM < FC
356 68794675 : else if (soil_moist > wilting_point) then
357 66702131 : feddes_et_reduction = frac_roots * (soil_moist - wilting_point) / (soil_moist_FC - wilting_point)
358 : ! SM <= PW
359 : else
360 : feddes_et_reduction = 0.0_dp
361 : end if
362 :
363 76647048 : END FUNCTION feddes_et_reduction
364 :
365 : ! ------------------------------------------------------------------
366 :
367 : ! NAME
368 : ! jarvis_et_reduction
369 :
370 : ! PURPOSE
371 : !> \brief stress factor for reducing evapotranspiration based on actual soil moisture
372 :
373 : !> \details The soil moisture stress factor is estimated based on the normalized soil water
374 : !> content. The normalized soil water content \f$ \theta_{norm} \f$ is estimated as:
375 : !> \f[ \theta_{norm} = \frac{\theta - \theta_\mathit{pwp}}
376 : !> {\theta_{sat} - \theta_{pwp}} \f]
377 : !> The ET reduction factor \f$ f \f$ is estimated as
378 : !> \f[ f = \left\{
379 : !> \begin{array}{lr}
380 : !> f_{roots} & if \theta_{norm} \ge jarvis\_sm\_threshold\_c1 \\
381 : !> f_{roots}\frac{\theta_{norm}}{jarvis\_sm\_threshold\_c1} &
382 : !> if \theta_{norm} < jarvis\_sm\_threshold\_c1 \\
383 : !> \end{array}
384 : !> \right. \f]
385 :
386 : ! INTENT(IN)
387 : !> \param[in] "real(dp) :: soil_moist" Soil moisture of each horizon [mm]
388 : !> \param[in] "real(dp) :: soil_moist_sat" saturated Soil moisture content [mm]
389 : !> \param[in] "real(dp) :: wilting_point" Permanent wilting point
390 : !> \param[in] "real(dp) :: frac_roots" Fraction of Roots in soil horizon is reduced [mm]
391 : !> \param[in] "real(dp) :: jarvis_thresh_c1" parameter C1 from Jarvis formulation
392 :
393 : ! RETURN
394 : !> \return real(dp) :: jarvis_et_reduction; et reduction factor
395 :
396 : ! HISTORY
397 : !> \authors Cueneyd Demirel, Matthias Zink
398 :
399 : !> \date March 2017
400 :
401 : ! Modifications:
402 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
403 : ! M. Cuneyd Demirel, Simon Stisen Jun 2020 - added Feddes and FC dependency on root fraction coefficient processCase(3) = 4
404 :
405 1782144 : elemental pure FUNCTION jarvis_et_reduction(soil_moist, soil_moist_sat, wilting_point, frac_roots, jarvis_thresh_c1)
406 : implicit none
407 :
408 : ! Soil moisture of each horizon [mm]
409 : real(dp), intent(in) :: soil_moist
410 :
411 : ! saturated Soil moisture content [mm]
412 : real(dp), intent(in) :: soil_moist_sat
413 :
414 : ! Permanent wilting point
415 : real(dp), intent(in) :: wilting_point
416 :
417 : ! Fraction of Roots in soil horizon is reduced [mm]
418 : real(dp), intent(in) :: frac_roots
419 :
420 : ! parameter C1 from Jarvis formulation
421 : real(dp), intent(in) :: jarvis_thresh_c1
422 :
423 : ! reference evapotranspiration in [mm]
424 : real(dp) :: jarvis_et_reduction
425 :
426 : ! normalized soil water content
427 1782144 : real(dp) :: theta_inorm
428 :
429 :
430 : ! Calculating normalized Soil Water Content
431 1782144 : theta_inorm = (soil_moist - wilting_point) / (soil_moist_sat - wilting_point)
432 :
433 : ! correct for numerical unaccuracies
434 1782144 : if (theta_inorm .LT. 0.0_dp) theta_inorm = 0.0_dp
435 1782144 : if (theta_inorm .GT. 1.0_dp) theta_inorm = 1.0_dp
436 :
437 : ! estimate fraction of ET demand based on root fraction and SM status using theta_inorm
438 : ! theta_inorm >= jarvis_thresh_c1
439 1782144 : if (theta_inorm .GE. jarvis_thresh_c1) then
440 815201 : jarvis_et_reduction = frac_roots
441 : ! 0 < theta_inorm < jarvis_thresh_c1
442 966943 : else if (theta_inorm .LT. jarvis_thresh_c1) then
443 966943 : jarvis_et_reduction = frac_roots * (theta_inorm / jarvis_thresh_c1)
444 : end if
445 :
446 151511952 : END FUNCTION jarvis_et_reduction
447 :
448 : END MODULE mo_soil_moisture
|
