Line data Source code
1 : !> \file mo_mpr_soilmoist.f90
2 : !> \brief \copybrief mo_mpr_soilmoist
3 : !> \details \copydetails mo_mpr_soilmoist
4 :
5 : !> \brief Multiscale parameter regionalization (MPR) for soil moisture
6 : !> \details This module contains all routines required for parametrizing soil moisture processes.
7 : !> \authors Stephan Thober, Rohini Kumar
8 : !> \date Dec 2012
9 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
10 : !! mHM is released under the LGPLv3+ license \license_note
11 : !> \ingroup f_mpr
12 : module mo_mpr_soilmoist
13 :
14 : use mo_kind, only : i4, dp
15 : use mo_message, only : message, error_message
16 :
17 : implicit none
18 :
19 : public :: mpr_sm
20 :
21 : private
22 :
23 : contains
24 : ! ----------------------------------------------------------------------------
25 :
26 : ! NAME
27 : ! mpr_sm
28 :
29 : ! PURPOSE
30 : !> \brief multiscale parameter regionalization for soil moisture
31 :
32 : !> \details This subroutine is a wrapper around all soil moisture
33 : !> parameter routines. This subroutine requires 13 parameters. These
34 : !> parameters have to correspond to the parameters in the original
35 : !> parameter array at the following locations: 10-12, 13-18, 27-30.
36 : !> Global parameters needed (see mhm_parameter.nml):
37 : !> - param( 1) = orgMatterContent_forest
38 : !> - param( 2) = orgMatterContent_impervious
39 : !> - param( 3) = orgMatterContent_pervious
40 : !> - param( 4) = PTF_lower66_5_constant
41 : !> - param( 5) = PTF_lower66_5_clay
42 : !> - param( 6) = PTF_lower66_5_Db
43 : !> - param( 7) = PTF_higher66_5_constant
44 : !> - param( 8) = PTF_higher66_5_clay
45 : !> - param( 9) = PTF_higher66_5_Db
46 : !> - param(10) = PTF_Ks_constant
47 : !> - param(11) = PTF_Ks_sand
48 : !> - param(12) = PTF_Ks_clay
49 : !> - param(13) = PTF_Ks_curveSlope
50 :
51 : ! INTENT(IN)
52 : !> \param[in] "real(dp), dimension(13) :: param" global parameters
53 : !> \param[in] "integer(i4), dimension(:) :: is_present" indicates whether soiltype is present
54 : !> \param[in] "integer(i4), dimension(:) :: nHorizons" Number of Horizons per soiltype
55 : !> \param[in] "integer(i4), dimension(:) :: nTillHorizons" Number of Tillage Horizons
56 : !> \param[in] "real(dp), dimension(:, :) :: sand" sand content
57 : !> \param[in] "real(dp), dimension(:, :) :: clay" clay content
58 : !> \param[in] "real(dp), dimension(:, :) :: DbM" mineral Bulk density
59 : !> \param[in] "integer(i4), dimension(:) :: ID0" cell ids at level 0
60 : !> \param[in] "integer(i4), dimension(:, :) :: soilId0" soil ids at level 0
61 : !> \param[in] "integer(i4), dimension(:) :: LCOVER0" land cover ids at level 0
62 :
63 : ! INTENT(OUT)
64 : !> \param[out] "real(dp), dimension(:, :, :) :: thetaS_till" saturated soil moisture tillage layer
65 : !> \param[out] "real(dp), dimension(:, :, :) :: thetaFC_till" field capacity tillage layer
66 : !> \param[out] "real(dp), dimension(:, :, :) :: thetaPW_till" permanent wilting point tillage layer
67 : !> \param[out] "real(dp), dimension(:, :) :: thetaS" saturated soil moisture
68 : !> \param[out] "real(dp), dimension(:, :) :: thetaFC" field capacity
69 : !> \param[out] "real(dp), dimension(:, :) :: thetaPW" permanent wilting point
70 : !> \param[out] "real(dp), dimension(:, :, :) :: Ks" saturated hydraulic conductivity
71 : !> \param[out] "real(dp), dimension(:, :, :) :: Db" Bulk density
72 : !> \param[out] "real(dp), dimension(:) :: KsVar_H0" rel. var. of Ks for horizontal flow
73 : !> \param[out] "real(dp), dimension(:) :: KsVar_V0" rel. var. of Ks for vertical flow
74 : !> \param[out] "real(dp), dimension(:) :: SMs_FC0" soil mositure deficit from
75 : !> field cap. w.r.t to saturation
76 :
77 : ! HISTORY
78 : !> \authors Stephan Thober, Rohini Kumar
79 :
80 : !> \date Dec 2012
81 :
82 : ! Modifications:
83 : ! Juliane Mai Oct 2013 - OLD parametrization --> param(1) = orgMatterContent_forest
84 : ! --> param(2) = orgMatterContent_impervious
85 : ! --> param(3) = orgMatterContent_pervious
86 : ! --> param(4:13) = ... orgMatterContent_forest = orgMatterContent_perv
87 : ! + delta_1 NEW parametrization
88 : ! --> param(1) = delta_1 --> param(2) = orgMatterContent_impervious
89 : ! --> param(3) = orgMatterContent_pervious --> param(4:13) = ...
90 : ! Matthias Zink Nov 2013 - documentation, inouts --> out moved constants to mhm_constants
91 : ! Stephan Thober Mar 2014 - separated cell loop from soil loop for better scaling in parallelisation
92 : ! David Schaefer Mar 2015 - Added dummy variable to avoid redundant computations
93 : ! -> Total number of instruction is reduced by ~25%
94 : ! (tested on packaged example/gnu48/{release,debug})
95 : ! Rohini Kumar Mar 2016 - changes for handling multiple soil database options
96 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
97 : ! M. Cuneyd Demirel, Simon Stisen Jun 2020 - added Feddes and FC dependency on root fraction coefficient processCase(3) = 4
98 : ! Rohini Kumar Oct 2021 - Neutron count module to mHM integrate into develop branch (5.11.2)
99 :
100 1242 : subroutine mpr_sm(param, processMatrix, is_present, nHorizons, nTillHorizons, sand, clay, DbM, ID0, soilId0, LCover0, &
101 414 : thetaS_till, thetaFC_till, thetaPW_till, thetaS, thetaFC, thetaPW, Ks, Db, KsVar_H0, KsVar_V0, SMs_FC0)
102 :
103 : use mo_common_constants, only : nodata_dp, nodata_i4
104 : use mo_mpr_constants, only : BulkDens_OrgMatter
105 : use mo_mpr_global_variables, only : iFlag_soilDB
106 : !$ use omp_lib
107 :
108 : implicit none
109 :
110 : ! global parameters
111 : real(dp), dimension(13), intent(in) :: param
112 :
113 : ! - matrix specifying user defined processes
114 : integer(i4), dimension(:, :), intent(in) :: processMatrix
115 :
116 :
117 : ! indicates whether soiltype is present
118 : integer(i4), dimension(:), intent(in) :: is_present
119 :
120 : ! Number of Horizons per soiltype
121 : integer(i4), dimension(:), intent(in) :: nHorizons
122 :
123 : ! Number of Tillage Horizons
124 : integer(i4), dimension(:), intent(in) :: nTillHorizons
125 :
126 : ! sand content
127 : real(dp), dimension(:, :), intent(in) :: sand
128 :
129 : ! clay content
130 : real(dp), dimension(:, :), intent(in) :: clay
131 :
132 : ! mineral Bulk density
133 : real(dp), dimension(:, :), intent(in) :: DbM
134 :
135 : ! cell ids at level 0
136 : integer(i4), dimension(:), intent(in) :: ID0
137 :
138 : ! soil ids at level 0
139 : integer(i4), dimension(:, :), intent(in) :: soilId0
140 :
141 : ! land cover ids at level 0
142 : integer(i4), dimension(:), intent(in) :: LCOVER0
143 :
144 : ! saturated soil moisture tillage layer
145 : real(dp), dimension(:, :, :), intent(out) :: thetaS_till
146 :
147 : ! field capacity tillage layer
148 : real(dp), dimension(:, :, :), intent(out) :: thetaFC_till
149 :
150 : ! permanent wilting point tillage layer
151 : real(dp), dimension(:, :, :), intent(out) :: thetaPW_till
152 :
153 : ! saturated soil moisture
154 : real(dp), dimension(:, :), intent(out) :: thetaS
155 :
156 : ! field capacity
157 : real(dp), dimension(:, :), intent(out) :: thetaFC
158 :
159 : ! permanent wilting point
160 : real(dp), dimension(:, :), intent(out) :: thetaPW
161 :
162 : ! saturated hydraulic conductivity
163 : real(dp), dimension(:, :, :), intent(out) :: Ks
164 :
165 : ! Bulk density
166 : real(dp), dimension(:, :, :), intent(out) :: Db
167 :
168 : ! rel. var. of Ks for horizontal flow
169 : real(dp), dimension(:), intent(out) :: KsVar_H0
170 :
171 : ! rel. var. of Ks for vertical flow
172 : real(dp), dimension(:), intent(out) :: KsVar_V0
173 :
174 : ! soil mositure deficit from
175 : ! field cap. w.r.t to saturation
176 : real(dp), dimension(:), intent(out) :: SMs_FC0
177 :
178 : ! loop index
179 : integer(i4) :: i
180 :
181 : ! loop index
182 : integer(i4) :: j
183 :
184 : ! loop index
185 : integer(i4) :: l
186 :
187 : ! dummy variable for storing soil class
188 : integer(i4) :: s
189 :
190 : integer(i4) :: tmp_minSoilHorizon
191 :
192 414 : real(dp) :: pM
193 :
194 414 : real(dp) :: pOM
195 :
196 : ! temporal saturated hydr. cond
197 414 : real(dp) :: Ks_tmp
198 :
199 : ! van Genuchten shape param
200 414 : real(dp) :: Genu_Mual_n
201 :
202 : ! van Genuchten shape param
203 414 : real(dp) :: Genu_Mual_alpha
204 :
205 414 : real(dp) :: tmp_orgMatterContent_forest
206 :
207 414 : real(dp) :: tmp_orgMatterContent_pervious
208 :
209 414 : real(dp) :: tmp_orgMatterContent_impervious
210 :
211 : ! total saturated soil mositure content
212 414 : real(dp), dimension(:), allocatable :: SMs_tot0
213 :
214 : ! maximum LCover class in L0
215 : integer(i4) :: max_LCover
216 :
217 : ! saturated hydraulic conductivity
218 414 : real(dp), dimension(:, :), allocatable :: Ks_non_till
219 :
220 : ! Case 1 and 4 is based on Jarvis https://doi.org/10.1016/0022-1694(89)90050-4
221 : ! Case 2 and 3 is based on Feddes https://doi.org/10.1016/0022-1694(76)90017-2
222 414 : select case (processMatrix(3, 1))
223 : case(1,2)
224 410 : tmp_orgMatterContent_forest = param(3) + param(1)
225 : case(3,4)
226 414 : tmp_orgMatterContent_forest = param(1)
227 : end select
228 :
229 414 : tmp_orgMatterContent_impervious = param(2)
230 414 : tmp_orgMatterContent_pervious = param(3)
231 : !write(*,*) 'tmp_orgMatterContent_forest = ', tmp_orgMatterContent_forest
232 : !write(*,*) 'tmp_orgMatterContent_impervious = ', tmp_orgMatterContent_impervious
233 : !write(*,*) 'tmp_orgMatterContent_pervious = ', tmp_orgMatterContent_pervious
234 :
235 611064 : tmp_minSoilHorizon = minval(nTillHorizons(:))
236 :
237 : ! allocatable local variables
238 : ! total saturated soil moisture content
239 1242 : allocate(SMs_tot0(size(ID0, 1)))
240 :
241 : ! some additional variables for iFlag_soil = 1
242 414 : if(iFlag_soilDB .eq. 1) then
243 0 : s = size(is_present, 1)
244 : ! note: although second dimension is not required
245 : ! we assign value 1 to be comparable to other assigned variables
246 0 : allocate(Ks_non_till(s, 1))
247 : end if
248 :
249 : ! initializing soil hydraulic properties related params
250 19168378 : KsVar_H0 = merge(0.0_dp, nodata_dp, ID0 .ne. nodata_i4)
251 19168378 : KsVar_V0 = merge(0.0_dp, nodata_dp, ID0 .ne. nodata_i4)
252 19168378 : SMs_tot0 = merge(0.0_dp, nodata_dp, ID0 .ne. nodata_i4)
253 19168378 : SMs_FC0 = merge(0.0_dp, nodata_dp, ID0 .ne. nodata_i4)
254 :
255 : ! initialization
256 1834848 : thetaS_till = 0.0_dp
257 1834848 : thetaFC_till = 0.0_dp
258 1834848 : thetaPW_till = 0.0_dp
259 1833606 : thetaS = 0.0_dp
260 1833606 : thetaFC = 0.0_dp
261 1833606 : thetaPW = 0.0_dp
262 5501232 : Ks = 0.0_dp
263 5501232 : Db = 0.0_dp
264 414 : if(allocated(Ks_non_till)) Ks_non_till = 0.0_dp
265 19167964 : max_LCover = maxval(LCOVER0)
266 : ! select case according to a given soil database flag
267 414 : SELECT CASE(iFlag_soilDB)
268 : ! classical mHM soil database format
269 : CASE(0)
270 : !$OMP PARALLEL default(shared)
271 : !$OMP DO &
272 : !$OMP PRIVATE( i, j, L, pOM, pM, Ks_tmp, Genu_Mual_alpha, Genu_Mual_n ) &
273 : !$OMP SCHEDULE( STATIC )
274 611064 : do i = 1, size(is_present)
275 610650 : if (is_present(i) .lt. 1) cycle
276 57318 : horizon : do j = 1, nHorizons(i)
277 : ! calculating vertical hydraulic conductivity
278 42678 : call hydro_cond(Ks_tmp, param(10 : 13), sand(i, j), clay(i, j))
279 170712 : Ks(i, j, :) = Ks_tmp
280 : ! calculating other soil hydraulic properties
281 : ! tillage horizons
282 653328 : if (j .le. nTillHorizons(i)) then
283 : ! LC class
284 56904 : do L = 1, max_LCover
285 : select case (L)
286 : case(1) ! forest
287 14226 : pOM = tmp_orgMatterContent_forest
288 : case(2) ! impervious
289 14226 : pOM = tmp_orgMatterContent_impervious !param(2)
290 : case(3) ! permeable
291 14226 : pOM = tmp_orgMatterContent_pervious
292 : case default
293 42678 : call error_message('Error mpr_sm: pOM used uninitialized.')
294 : end select
295 42678 : pM = 100.0_dp - pOM
296 : ! bulk density acording to Rawl's (1982) paper
297 42678 : Db(i, j, L) = 100.0_dp / ((pOM / BulkDens_OrgMatter) + (pM / DbM(i, j)))
298 : ! Effect of organic matter content
299 : ! This is taken into account in a simplified form by using
300 : ! the ratio of(Bd / BdOM)
301 42678 : Ks(i, j, L) = Ks(i, j, L) * (DbM(i, j) / Db(i, j, L))
302 : ! estimated SMs_till & van Genuchten's shape parameter (n)
303 42678 : call Genuchten(thetaS_till(i, j, L), Genu_Mual_n, Genu_Mual_alpha, &
304 42678 : param(4 : 9), sand(i, j), clay(i, j), Db(i, j, L))
305 : ! estimating field capacity
306 42678 : call field_cap(thetaFC_till(i, j, L), Ks(i, j, L), thetaS_till(i, j, L), Genu_Mual_n)
307 : ! estimating permanent wilting point
308 56904 : call PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS_till(i, j, L), thetaPW_till(i, j, L))
309 : end do
310 : ! deeper layers
311 : else
312 : ! estimate SMs & van Genuchten's shape parameter (n)
313 28452 : call Genuchten(thetaS(i, j - tmp_minSoilHorizon), Genu_Mual_n, Genu_Mual_alpha, &
314 56904 : param(4 : 9), sand(i, j), clay(i, j), DbM(i, j))
315 : ! estimate field capacity
316 28452 : call field_cap(thetaFC(i, j - tmp_minSoilHorizon), &
317 28452 : Ks(i, j, 1), thetaS(i, j - tmp_minSoilHorizon), Genu_Mual_n)
318 : ! estimate permanent wilting point
319 : call PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS(i, j - tmp_minSoilHorizon), &
320 28452 : thetaPW(i, j - tmp_minSoilHorizon))
321 : end if
322 : end do horizon
323 : end do
324 : !$OMP END DO
325 :
326 : ! calculate other soil properties at each location [L0] for regionalising model parameters
327 : !$OMP DO PRIVATE( s, j ) SCHEDULE( STATIC )
328 19168378 : cellloop : do i = 1, size(soilId0, 1) ! >> here = ncells0
329 19167550 : s = soilId0(i, 1) ! >> in this case the second dimension of soilId0 = 1
330 76670200 : do j = 1, nHorizons(s)
331 76670200 : if (j .le. nTillHorizons(s)) then
332 : ! Soil properties over the whole soil coloum depth
333 19167550 : KsVar_H0(i) = KsVar_H0(i) + thetaS_till(s, j, LCover0(i)) * Ks(s, j, LCover0(i))
334 19167550 : KsVar_V0(i) = KsVar_V0(i) + thetaS_till(s, j, LCover0(i)) / Ks(s, j, LCover0(i))
335 19167550 : SMs_FC0(i) = SMs_FC0(i) + thetaFC_till(s, j, LCover0(i))
336 19167550 : SMs_tot0(i) = SMs_tot0(i) + thetaS_till (s, j, LCover0(i))
337 : else
338 : ! soil_properties over the whole soil column
339 38335100 : KsVar_H0(i) = KsVar_H0(i) + thetaS(s, j - tmp_minSoilHorizon) * Ks(s, j, 1)
340 38335100 : KsVar_V0(i) = KsVar_V0(i) + thetaS(s, j - tmp_minSoilHorizon) / Ks(s, j, 1)
341 38335100 : SMs_FC0(i) = SMs_FC0(i) + thetaFC(s, j - tmp_minSoilHorizon)
342 38335100 : SMs_tot0(i) = SMs_tot0(i) + thetaS (s, j - tmp_minSoilHorizon)
343 : end if
344 : end do
345 : ! ------------------------------------------------------------------
346 : ! DETERMINE RELATIVE VARIABILITIES OF
347 : ! Ks FOR HORIZONTAL FLOW (KsVar_H)
348 : ! &
349 : ! Ks FOR VERTICAL FLOW (KsVar_V)
350 : ! ------------------------------------------------------------------
351 : ! soil moisture saturation deficit relative to the field capacity soil moisture
352 19167550 : SMs_FC0(i) = (SMs_tot0(i) - SMs_FC0(i)) / SMs_tot0(i)
353 : ! Ks variability over the whole soil coloum depth for
354 : ! both horizontal and vertical flows including relative variabilities
355 19167550 : KsVar_H0(i) = KsVar_H0(i) / SMs_tot0(i) / param(13)
356 19167964 : KsVar_V0(i) = SMs_tot0(i) / KsVar_V0(i) / param(13)
357 : end do cellloop
358 : !$OMP END DO
359 : !$OMP END PARALLEL
360 :
361 : ! to handle multiple soil horizons with unique soil class
362 : CASE(1)
363 0 : do i = 1, size(is_present)
364 0 : if (is_present(i) .lt. 1) cycle
365 : ! **** FOR THE TILLAGE TYPE OF SOIL *****
366 : ! there is actually no soil horizons/soil type in this case
367 : ! but we assign of j = 1 to use variables as defined in the classical option (iFlag_soil = 0)
368 0 : do j = 1, 1
369 : ! calculating vertical hydraulic conductivity
370 0 : call hydro_cond(Ks_tmp, param(10 : 13), sand(i, j), clay(i, j))
371 0 : Ks_non_till(i, j) = Ks_tmp ! >> non-till
372 0 : Ks(i, j, :) = Ks_tmp ! >> till layers
373 : ! calculating other soil hydraulic properties
374 : ! tillage horizons properties depending on the LC class
375 0 : do L = 1, max_LCover
376 : select case (L)
377 : case(1) ! forest
378 0 : pOM = tmp_orgMatterContent_forest
379 : case(2) ! impervious
380 0 : pOM = tmp_orgMatterContent_impervious !param(2)
381 : case(3) ! permeable
382 0 : pOM = tmp_orgMatterContent_pervious
383 : case default
384 0 : call error_message('Error mpr_sm: pOM used is not initialized.')
385 : end select
386 0 : pM = 100.0_dp - pOM
387 : ! bulk density acording to Rawl's (1982) paper
388 0 : Db(i, j, L) = 100.0_dp / ((pOM / BulkDens_OrgMatter) + (pM / DbM(i, j)))
389 : ! Effect of organic matter content on Ks estimates
390 : ! This is taken into account in a simplified form by using
391 : ! the ratio of (Bd/BdOM)
392 0 : Ks_tmp = Ks_tmp * (DbM(i, j) / Db(i, j, L))
393 0 : Ks(i, j, L) = Ks_tmp
394 : ! estimated SMs_till & van Genuchten's shape parameter (n)
395 0 : call Genuchten(thetaS_till(i, j, L), Genu_Mual_n, Genu_Mual_alpha, &
396 0 : param(4 : 9), sand(i, j), clay(i, j), Db(i, j, L))
397 : ! estimating field capacity
398 0 : call field_cap(thetaFC_till(i, j, L), Ks_tmp, thetaS_till(i, j, L), Genu_Mual_n)
399 : ! estimating permanent wilting point
400 0 : call PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS_till(i, j, L), thetaPW_till(i, j, L))
401 : end do
402 :
403 : ! *** FOR NON-TILLAGE TYPE OF SOILS ***
404 : ! note j = 1
405 : ! since Ks_tmp has changed earlier ... get the original Ks once again
406 0 : Ks_tmp = Ks_non_till(i, j)
407 : ! estimate SMs & van Genuchten's shape parameter (n)
408 0 : call Genuchten(thetaS(i, j), Genu_Mual_n, Genu_Mual_alpha, param(4 : 9), sand(i, j), clay(i, j), DbM(i, j))
409 : ! estimate field capacity
410 0 : call field_cap(thetaFC(i, j), Ks_tmp, thetaS(i, j), Genu_Mual_n)
411 : ! estimate permanent wilting point
412 0 : call PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS(i, j), thetaPW(i, j))
413 :
414 : end do ! >> HORIZON
415 : end do ! >> SOIL TYPE
416 :
417 : ! calculate other soil properties at each location [L0] for regionalising model parameters
418 0 : do i = 1, size(soilId0, 1) !! over all cells
419 0 : do j = 1, size(soilId0, 2) !! over horizons
420 0 : s = soilId0(i, j)
421 0 : if (j .le. nTillHorizons(1)) then
422 : ! soil properties over the whole soil coloum depth
423 0 : KsVar_H0(i) = KsVar_H0(i) + thetaS_till (s, 1, LCover0(i)) * Ks(s, 1, LCover0(i))
424 0 : KsVar_V0(i) = KsVar_V0(i) + thetaS_till (s, 1, LCover0(i)) / Ks(s, 1, LCover0(i))
425 0 : SMs_FC0(i) = SMs_FC0 (i) + thetaFC_till(s, 1, LCover0(i))
426 0 : SMs_tot0(i) = SMs_tot0(i) + thetaS_till (s, 1, LCover0(i))
427 : else
428 : ! soil_properties over the whole soil column
429 0 : KsVar_H0(i) = KsVar_H0(i) + thetaS (s, 1) * Ks_non_till(s, 1)
430 0 : KsVar_V0(i) = KsVar_V0(i) + thetaS (s, 1) / Ks_non_till(s, 1)
431 0 : SMs_FC0(i) = SMs_FC0 (i) + thetaFC(s, 1)
432 0 : SMs_tot0(i) = SMs_tot0(i) + thetaS (s, 1)
433 : end if
434 : end do
435 : ! ------------------------------------------------------------------
436 : ! DETERMINE RELATIVE VARIABILITIES OF
437 : ! Ks FOR HORIZONTAL FLOW (KsVar_H) & Ks FOR VERTICAL FLOW (KsVar_V)
438 : ! ------------------------------------------------------------------
439 : ! soil moisture saturation deficit relative to the field capacity soil moisture
440 0 : SMs_FC0(i) = (SMs_tot0(i) - SMs_FC0(i)) / SMs_tot0(i)
441 : ! Ks variability over the whole soil coloum depth for
442 : ! both horizontal and vertical flows including relative variabilities
443 0 : KsVar_H0(i) = KsVar_H0(i) / SMs_tot0(i) / param(13)
444 0 : KsVar_V0(i) = SMs_tot0(i) / KsVar_V0(i) / param(13)
445 : end do
446 :
447 : CASE DEFAULT
448 414 : call error_message('***ERROR: iFlag_soilDB option given does not exist. Only 0 and 1 is taken at the moment.')
449 : END SELECT
450 :
451 : ! free space **
452 414 : deallocate(SMs_tot0)
453 414 : if(allocated(Ks_non_till)) deallocate(Ks_non_till)
454 :
455 414 : end subroutine mpr_sm
456 :
457 : ! ------------------------------------------------------------------
458 :
459 : ! NAME
460 : ! PWP
461 :
462 : ! PURPOSE
463 : !> \brief Permanent Wilting point
464 :
465 : !> \details This subroutine calculates the permanent wilting
466 : !> point according to Zacharias et al. (2007, Soil Phy.) and
467 : !> using van Genuchten 1980's equation. For the water retention curve at
468 : !> a matrix potential of -1500 kPa, it is assumed that thetaR = 0.
469 : !> ADDITIONAL INFORMATION
470 : !> Zacharias et al. 2007, Soil Phy.
471 :
472 : ! INTENT(IN)
473 : !> \param[in] "real(dp) :: Genu_Mual_n" - Genuchten shape parameter
474 : !> \param[in] "real(dp) :: Genu_Mual_alpha" - Genuchten shape parameter
475 : !> \param[in] "real(dp) :: thetaS" - saturated water content
476 :
477 : ! INTENT(OUT)
478 : !> \param[out] "real(dp) :: thetaPWP" - Permanent Wilting point
479 :
480 : ! HISTORY
481 : !> \authors Stephan Thober, Rohini Kumar
482 :
483 : !> \date Dec, 2012
484 :
485 : ! Modifications:
486 : ! Matthias Zink Nov 2013 - documentation, moved constants to mhm_constants
487 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
488 :
489 71130 : elemental pure subroutine PWP(Genu_Mual_n, Genu_Mual_alpha, thetaS, thetaPWP)
490 :
491 414 : use mo_mpr_constants, only : PWP_c, PWP_matPot_ThetaR
492 :
493 : implicit none
494 :
495 : ! - Genuchten shape parameter
496 : real(dp), intent(in) :: Genu_Mual_n
497 :
498 : ! - Genuchten shape parameter
499 : real(dp), intent(in) :: Genu_Mual_alpha
500 :
501 : ! - saturated water content
502 : real(dp), intent(in) :: thetaS
503 :
504 : ! - Permanent Wilting point
505 : real(dp), intent(out) :: thetaPWP
506 :
507 71130 : real(dp) :: x
508 :
509 : ! Genuchten shape parameter
510 71130 : real(dp) :: Genu_Mual_m
511 :
512 :
513 71130 : Genu_Mual_m = PWP_c - (PWP_c / Genu_Mual_n)
514 71130 : x = PWP_c + exp(Genu_Mual_n * log(Genu_Mual_alpha * PWP_matPot_ThetaR))
515 71130 : x = exp(Genu_Mual_m * log(x))
516 : ! constrain
517 0 : if (x < 1.0_dp) x = 1.0_dp
518 71130 : thetaPWP = thetaS / x
519 :
520 71130 : end subroutine PWP
521 :
522 : ! ----------------------------------------------------------------------------
523 :
524 : ! NAME
525 : ! field_cap
526 :
527 : ! PURPOSE
528 : !> \brief calculates the field capacity
529 :
530 : !> \details estimate Field capacity; FC -- Flux based
531 : !> approach (Twarakavi, et. al. 2009, WRR)
532 : !> According to the
533 : !> above reference FC is defined as the soil water content at
534 : !> which the drainage from a profile ceases under natural
535 : !> conditions. Since drainage from a soil profile in a simulation
536 : !> never becomes zero, we assume that drainage ceases when the
537 : !> bottom flux from the soil reaches a value that is equivalent to
538 : !> the minimum amount of precipitation that could be recorded
539 : !> (i.e. 0.01 cm/d == 1 mm/d). It is assumed that ThetaR = 0.0_dp
540 : !> ADDITIONAL INFORMATION
541 : !> Twarakavi, et. al. 2009, WRR
542 :
543 : ! INTENT(OUT)
544 : !> \param[out] "real(dp) :: thetaFC" - Field capacity
545 :
546 : ! INTENT(IN)
547 : !> \param[in] "real(dp) :: Ks" - saturated hydraulic conductivity
548 : !> \param[in] "real(dp) :: thetaS" - saturated water content
549 : !> \param[in] "real(dp) :: Genu_Mual_n" - Genuchten shape parameter
550 :
551 : ! HISTORY
552 : !> \authors Stephan Thober, Rohini Kumar
553 :
554 : !> \date Dec 2012
555 :
556 : ! Modifications:
557 : ! Matthias Zink Nov 2013 - documentation, moved constants to mhm_constants
558 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
559 :
560 71130 : elemental pure subroutine field_cap(thetaFC, Ks, thetaS, Genu_Mual_n)
561 :
562 71130 : use mo_mpr_constants, only : field_cap_c1, field_cap_c2
563 :
564 : implicit none
565 :
566 : ! - saturated hydraulic conductivity
567 : real(dp), intent(in) :: Ks
568 :
569 : ! - saturated water content
570 : real(dp), intent(in) :: thetaS
571 :
572 : ! - Genuchten shape parameter
573 : real(dp), intent(in) :: Genu_Mual_n
574 :
575 : ! - Field capacity
576 : real(dp), intent(out) :: thetaFC
577 :
578 71130 : real(dp) :: x
579 :
580 :
581 71130 : x = (field_cap_c1) * (field_cap_c2 + log10(Ks))
582 71130 : thetaFC = thetaS * exp(x * log(Genu_Mual_n))
583 :
584 71130 : end subroutine field_cap
585 :
586 : ! ----------------------------------------------------------------------------
587 :
588 : ! NAME
589 : ! Genuchten
590 :
591 : ! PURPOSE
592 : !> \brief calculates the Genuchten shape parameter
593 :
594 : !> \details estimate SMs_till & van Genuchten's shape parameter (n)
595 : !> (Zacharias et al, 2007, soil Phy.)
596 : !> Global parameters needed (see mhm_parameter.nml):
597 : !> - param( 1) = PTF_lower66_5_constant
598 : !> - param( 2) = PTF_lower66_5_clay
599 : !> - param( 3) = PTF_lower66_5_Db
600 : !> - param( 4) = PTF_higher66_5_constant
601 : !> - param( 5) = PTF_higher66_5_clay
602 : !> - param( 6) = PTF_higher66_5_Db
603 : !> ADDITIONAL INFORMATION
604 : !> Zacharias et al, 2007, soil Phy.
605 :
606 : ! INTENT(OUT)
607 : !> \param[out] "real(dp) :: thetaS" - saturated water content
608 : !> \param[out] "real(dp) :: Genu_Mual_n" - van Genuchten shape parameter
609 : !> \param[out] "real(dp) :: Genu_Mual_alpha" - van Genuchten shape parameter
610 :
611 : ! INTENT(IN)
612 : !> \param[in] "real(dp), dimension(6) :: param" parameters
613 : !> \param[in] "real(dp) :: sand" - [%] sand content
614 : !> \param[in] "real(dp) :: clay" - [%] clay content
615 : !> \param[in] "real(dp) :: Db" - [10^3 kg/m3] bulk density
616 :
617 : ! HISTORY
618 : !> \authors Stephan Thober, Rohini Kumar
619 :
620 : !> \date Dec 2012
621 :
622 : ! Modifications:
623 : ! Rohini Kumar Mar 2014 - ThetaS limit changed from 0 to 0.001
624 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
625 :
626 71130 : subroutine Genuchten(thetaS, Genu_Mual_n, Genu_Mual_alpha, param, sand, clay, Db)
627 :
628 71130 : use mo_mpr_constants, only : vGenuchtenN_c1, vGenuchtenN_c10, vGenuchtenN_c11, vGenuchtenN_c12, vGenuchtenN_c13, &
629 : vGenuchtenN_c14, vGenuchtenN_c15, vGenuchtenN_c16, vGenuchtenN_c17, vGenuchtenN_c18, &
630 : vGenuchtenN_c2, vGenuchtenN_c3, vGenuchtenN_c4, vGenuchtenN_c5, vGenuchtenN_c6, &
631 : vGenuchtenN_c7, vGenuchtenN_c8, vGenuchtenN_c9, vGenuchten_sandtresh
632 :
633 : implicit none
634 :
635 : ! parameters
636 : real(dp), dimension(6), intent(in) :: param
637 :
638 : ! - [%] sand content
639 : real(dp), intent(in) :: sand
640 :
641 : ! - [%] clay content
642 : real(dp), intent(in) :: clay
643 :
644 : ! - [10^3 kg/m3] bulk density
645 : real(dp), intent(in) :: Db
646 :
647 : ! - saturated water content
648 : real(dp), intent(out) :: thetaS
649 :
650 : ! - van Genuchten shape parameter
651 : real(dp), intent(out) :: Genu_Mual_n
652 :
653 : ! - van Genuchten shape parameter
654 : real(dp), intent(out) :: Genu_Mual_alpha
655 :
656 : ! temporal variable
657 71130 : real(dp) :: x
658 :
659 :
660 : ! estimate SMs_till & van Genuchten's parameters (alpha and n)
661 71130 : if (sand < vGenuchten_sandtresh) then
662 48360 : thetaS = param(1) + param(2) * clay + param(3) * Db
663 : Genu_Mual_n = vGenuchtenN_c1 - vGenuchtenN_c2 * (sand**(vGenuchtenN_c3)) + &
664 48360 : vGenuchtenN_c4 * (clay**(vGenuchtenN_c5))
665 : x = vGenuchtenN_c6 + vGenuchtenN_c7 * sand + vGenuchtenN_c8 * clay - &
666 48360 : vGenuchtenN_c9 * Db
667 : else
668 22770 : thetaS = param(4) + param(5) * clay + param(6) * Db
669 : Genu_Mual_n = vGenuchtenN_c10 + vGenuchtenN_c11 * (sand**(vGenuchtenN_c12)) + &
670 22770 : vGenuchtenN_c13 * (clay**(vGenuchtenN_c14))
671 : x = vGenuchtenN_c15 + vGenuchtenN_c16 * sand + vGenuchtenN_c17 * clay - &
672 22770 : vGenuchtenN_c18 * Db
673 : end if
674 :
675 : ! Mualem alpha
676 71130 : Genu_Mual_alpha = exp(x)
677 :
678 : ! hard coded limits, according to (Zacharias et al, 2007, soil Phy.)
679 71130 : if (thetaS < 0.01_dp) then
680 0 : call message('thetaS below threshold limit 1e-2, reset.')
681 : ! Put constrains on theta_S
682 0 : thetaS = 0.01_dp
683 : end if
684 71130 : if (thetaS > 1.0_dp) then
685 0 : call message('thetaS above 1, reset.')
686 : ! Put constrains on theta_S
687 0 : thetaS = 1.0_dp
688 : end if
689 71130 : if (Genu_Mual_n < 1.01000_dp) then
690 0 : call message('Genu_Mual_n below threshold limit 1.01, reset.')
691 0 : Genu_Mual_n = 1.01000_dp
692 : end if
693 71130 : if (Genu_Mual_alpha < 0.00001_dp) then
694 0 : call message('Genu_Mual_alpha below threshold limit 1e-5, reset.')
695 0 : Genu_Mual_alpha = 0.00001_dp
696 : end if
697 :
698 71130 : end subroutine Genuchten
699 :
700 : ! ----------------------------------------------------------------------------
701 :
702 : ! NAME
703 : ! hydro_cond
704 :
705 : ! PURPOSE
706 : !> \brief calculates the hydraulic conductivity Ks
707 :
708 : !> \details By default save this value of Ks, particularly for the
709 : !> deeper layers where OM content plays relatively low or no role
710 : !> Global parameters needed (see mhm_parameter.nml):
711 : !> - param(1) = PTF_Ks_constant
712 : !> - param(2) = PTF_Ks_sand
713 : !> - param(3) = PTF_Ks_clay
714 : !> - param(4) = PTF_Ks_curveSlope
715 : !> ADDITIONAL INFORMATION
716 : !> Written, Stephan Thober, Dec 2012
717 :
718 : ! INTENT(OUT)
719 : !> \param[out] "real(dp) :: KS"
720 :
721 : ! INTENT(IN)
722 : !> \param[in] "real(dp), dimension(4) :: param"
723 : !> \param[in] "real(dp) :: sand" - [%] sand content
724 : !> \param[in] "real(dp) :: clay" - [%] clay content
725 :
726 : ! HISTORY
727 : !> \authors Stephan Thober, Rohini Kumar
728 :
729 : !> \date Dec 2012
730 :
731 : ! Modifications:
732 : ! Matthias Zink Nov 2013 - documentation, moved constants to mhm_constants
733 : ! Matthias Cuntz Jun 2014 - suggested to fix param(4)
734 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
735 :
736 42678 : subroutine hydro_cond(KS, param, sand, clay)
737 :
738 71130 : use mo_mpr_constants, only : Ks_c
739 :
740 : implicit none
741 :
742 : real(dp), dimension(4), intent(in) :: param
743 :
744 : ! - [%] sand content
745 : real(dp), intent(in) :: sand
746 :
747 : ! - [%] clay content
748 : real(dp), intent(in) :: clay
749 :
750 : real(dp), intent(out) :: KS
751 :
752 : ! temporal variable
753 42678 : real(dp) :: x
754 :
755 :
756 : ! saturated vertical hydraulic conductivity, Ks (cm/d)
757 : ! from Cosby et. al. (WRR 1984) Table 4
758 : ! param(4) is the unit conversion from inch/h to cm/d and should be a constant.
759 : ! Fix it in the namelist, i.e. in
760 : ! mhm_parameter.nml set the 4th value (=FLAG) to 0 and the third value to 60.96
761 : ! PTF_Ks_curveSlope = 60.96, 60.96, 60.96, 0, 1
762 42678 : x = param(1) + param(2) * sand - param(3) * clay
763 42678 : Ks = param(4) * exp(X * log(Ks_c))
764 :
765 42678 : if (Ks < 1.10_dp) then
766 0 : call message('JMJMJM-Ks-BAD')
767 : end if
768 :
769 : ! minimum value of Ks = 1.1cm/d
770 42678 : if (Ks < 1.10_dp) Ks = 1.10_dp
771 :
772 42678 : end subroutine hydro_cond
773 :
774 : end module mo_mpr_soilmoist
|