Line data Source code
1 : !> \file mo_mpr_smhorizons.f90
2 : !> \brief \copybrief mo_mpr_smhorizons
3 : !> \details \copydetails mo_mpr_smhorizons
4 :
5 : !> \brief setting up the soil moisture horizons
6 : !> \details This module sets up the soil moisture horizons
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_SMhorizons
13 :
14 : use mo_kind, only : i4, dp
15 : use mo_common_constants, only : nodata_dp
16 :
17 : implicit none
18 :
19 : public :: mpr_SMhorizons
20 :
21 : private
22 :
23 : contains
24 :
25 : ! ----------------------------------------------------------------------------
26 :
27 : ! NAME
28 : ! mpr_SMhorizons
29 :
30 : ! PURPOSE
31 : !> \brief upscale soil moisture horizons
32 :
33 : !> \details calculate soil properties at the level 1.
34 : !> Global parameters needed (see mhm_parameter.nml):
35 : !> - param(1) = rootFractionCoefficient_forest
36 : !> - param(2) = rootFractionCoefficient_impervious
37 : !> - param(3) = rootFractionCoefficient_pervious
38 : !> - param(4) = infiltrationShapeFactor
39 :
40 : ! INTENT(IN)
41 : !> \param[in] "real(dp), dimension(:) :: param" parameters
42 : !> \param[in] "integer(i4), dimension(:, :) :: processMatrix" - matrix specifying user defined processes
43 : !> \param[in] "integer(i4) :: iFlag_soil" - flags for handling multiple soil databases
44 : !> \param[in] "integer(i4) :: nHorizons_mHM" - number of horizons to model
45 : !> \param[in] "real(dp), dimension(:) :: HorizonDepth" [10^-3 m] horizon depth from
46 : !> surface, postive downwards
47 : !> \param[in] "integer(i4), dimension(:) :: LCOVER0" Land cover at level 0
48 : !> \param[in] "integer(i4), dimension(:, :) :: soilID0" soil ID at level 0
49 : !> \param[in] "integer(i4), dimension(:) :: nHorizons" horizons per soil type
50 : !> \param[in] "integer(i4), dimension(:) :: nTillHorizons" Number of Tillage horizons
51 : !> \param[in] "real(dp), dimension(:, :, :) :: thetaS_till" saturated water content of soil
52 : !> horizons upto tillage depth,
53 : !> f(OM, management)
54 : !> \param[in] "real(dp), dimension(:, :, :) :: thetaFC_till" Field capacity of tillage
55 : !> layers; LUC dependent,
56 : !> f(OM, management)
57 : !> \param[in] "real(dp), dimension(:, :, :) :: thetaPW_till" Permament wilting point of
58 : !> tillage layers; LUC dependent,
59 : !> f(OM, management)
60 : !> \param[in] "real(dp), dimension(:, :) :: thetaS" saturated water content of soil
61 : !> horizons after tillage depth
62 : !> \param[in] "real(dp), dimension(:, :) :: thetaFC" Field capacity of deeper layers
63 : !> \param[in] "real(dp), dimension(:, :) :: thetaPW" Permanent wilting point of
64 : !> deeper layers
65 : !> \param[in] "real(dp), dimension(:, :, :) :: Wd" weights of mHM Horizons
66 : !> according to horizons provided
67 : !> in soil database
68 : !> \param[in] "real(dp), dimension(:, :, :) :: Db" Bulk density
69 : !> \param[in] "real(dp), dimension(:, :) :: DbM" mineral Bulk density
70 : !> \param[in] "real(dp), dimension(:) :: RZdepth" [mm] Total soil depth
71 : !> \param[in] "logical, dimension(:, :) :: mask0" mask at L0
72 : !> \param[in] "integer(i4), dimension(:) :: cell_id0" Cell ids of hi res field
73 : !> \param[in] "integer(i4), dimension(:) :: upp_row_L1" Upper row of hi res block
74 : !> \param[in] "integer(i4), dimension(:) :: low_row_L1" Lower row of hi res block
75 : !> \param[in] "integer(i4), dimension(:) :: lef_col_L1" Left column of hi res block
76 : !> \param[in] "integer(i4), dimension(:) :: rig_col_L1" Right column of hi res block
77 : !> \param[in] "integer(i4), dimension(:) :: nL0_in_L1" Number of L0 cells within a L1 cel
78 :
79 : ! INTENT(INOUT)
80 : !> \param[inout] "real(dp), dimension(:, :) :: L1_beta" Parameter that determines the
81 : !> relative contribution to SM, upscaled
82 : !> Bulk density
83 : !> \param[inout] "real(dp), dimension(:, :) :: L1_SMs" [10^-3 m] depth of saturated SM cont
84 : !> \param[inout] "real(dp), dimension(:, :) :: L1_FC" [10^-3 m] field capacity
85 : !> \param[inout] "real(dp), dimension(:, :) :: L1_PW" [10^-3 m] permanent wilting point
86 : !> \param[inout] "real(dp), dimension(:, :) :: L1_fRoots" fraction of roots in soil horizons
87 :
88 : ! HISTORY
89 : !> \authors Luis Samaniego, Rohini Kumar, Stephan Thober
90 :
91 : !> \date Dec 2012
92 :
93 : ! Modifications:
94 : ! Stephan Thober Jan 2013 - updated calling sequence for upscaling operators
95 : ! Juliane Mai Oct 2013 - OLD parametrization
96 : ! --> param(1) = rootFractionCoefficient_forest
97 : ! --> param(2) = rootFractionCoefficient_impervious
98 : ! --> param(3) = rootFractionCoefficient_pervious
99 : ! --> param(4) = infiltrationShapeFactor
100 : ! -------------------------------
101 : ! rootFractionCoeff_perv = rootFractionCoeff_forest - delta_1
102 : ! -------------------------------
103 : ! NEW parametrization
104 : ! --> param(1) = rootFractionCoefficient_forest
105 : ! --> param(2) = rootFractionCoefficient_impervious
106 : ! --> param(3) = delta_1
107 : ! --> param(4) = infiltrationShapeFactor
108 : ! ! if processMatrix(3,1) = 3 and 4 additionally
109 : ! --> param(5) = rootFractionCoefficient_sand
110 : ! --> param(6) = rootFractionCoefficient_clay
111 : ! --> param(7) = FCmin_glob
112 : ! --> param(8) = FCdelta_glob
113 : ! Stephan Thober Mar 2014 - added omp parallelization
114 : ! Rohini Kumar Mar 2016 - changes for handling multiple soil database options
115 : ! M. Cuneyd Demirel, Simon Stisen Apr 2017 - added FC dependency on root fraction coefficient
116 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
117 : ! M. Cuneyd Demirel, Simon Stisen Jun 2020 - added Feddes and FC dependency on root fraction coefficient processCase(3) = 4
118 : ! M. Cuneyd Demirel, Simon Stisen Feb 2021 - Bug fix normalization of FCnorm
119 :
120 828 : subroutine mpr_SMhorizons(param, processMatrix, iFlag_soil, nHorizons_mHM, HorizonDepth, LCOVER0, soilID0, nHorizons, &
121 1242 : nTillHorizons, thetaS_till, thetaFC_till, thetaPW_till, thetaS, thetaFC, thetaPW, Wd, Db, &
122 828 : DbM, RZdepth, mask0, cell_id0, upp_row_L1, low_row_L1, lef_col_L1, rig_col_L1, nL0_in_L1, &
123 828 : L1_beta, L1_SMs, L1_FC, L1_PW, L1_fRoots, &
124 : ! neutron count
125 414 : latWat_till , & ! lattice water upto tillage depth
126 828 : COSMIC_L3_till, & ! COSMIC parameter L3 upto tillage depth
127 414 : latWat , & ! lattice water
128 414 : COSMIC_L3 , & ! COSMIC paramter L3
129 414 : L1_bulkDens , & ! L bulk density
130 414 : L1_latticeWater,& ! L1 lattice water content
131 414 : L1_COSMICL3 ) ! L1 COSMIC L3 parameter from neutron module
132 :
133 : use mo_message, only : message, error_message
134 : use mo_string_utils, only : num2str
135 : use mo_upscaling_operators, only : upscale_harmonic_mean
136 : !$ use omp_lib
137 :
138 : implicit none
139 :
140 : ! parameters
141 : real(dp), dimension(:), intent(in) :: param
142 :
143 : ! - matrix specifying user defined processes
144 : integer(i4), dimension(:, :), intent(in) :: processMatrix
145 :
146 : ! - flags for handling multiple soil databases
147 : integer(i4), intent(in) :: iFlag_soil
148 :
149 : ! - number of horizons to model
150 : integer(i4), intent(in) :: nHorizons_mHM
151 :
152 : ! [10^-3 m] horizon depth from
153 : ! surface, postive downwards
154 : real(dp), dimension(:), intent(in) :: HorizonDepth
155 :
156 : ! Land cover at level 0
157 : integer(i4), dimension(:), intent(in) :: LCOVER0
158 :
159 : ! soil ID at level 0
160 : integer(i4), dimension(:, :), intent(in) :: soilID0
161 :
162 : ! horizons per soil type
163 : integer(i4), dimension(:), intent(in) :: nHorizons
164 :
165 : ! Number of Tillage horizons
166 : integer(i4), dimension(:), intent(in) :: nTillHorizons
167 :
168 : ! saturated water content of soil
169 : ! horizons upto tillage depth,
170 : ! f(OM, management)
171 : real(dp), dimension(:, :, :), intent(in) :: thetaS_till
172 :
173 : ! Field capacity of tillage
174 : ! layers; LUC dependent,
175 : ! f(OM, management)
176 : real(dp), dimension(:, :, :), intent(in) :: thetaFC_till
177 :
178 : ! Permament wilting point of
179 : ! tillage layers; LUC dependent,
180 : ! f(OM, management)
181 : real(dp), dimension(:, :, :), intent(in) :: thetaPW_till
182 :
183 : ! saturated water content of soil
184 : ! horizons after tillage depth
185 : real(dp), dimension(:, :), intent(in) :: thetaS
186 :
187 : ! Field capacity of deeper layers
188 : real(dp), dimension(:, :), intent(in) :: thetaFC
189 :
190 : ! Permanent wilting point of
191 : ! deeper layers
192 : real(dp), dimension(:, :), intent(in) :: thetaPW
193 :
194 : ! weights of mHM Horizons
195 : ! according to horizons provided
196 : ! in soil database
197 : real(dp), dimension(:, :, :), intent(in) :: Wd
198 :
199 : ! Bulk density
200 : real(dp), dimension(:, :, :), intent(in) :: Db
201 :
202 : ! mineral Bulk density
203 : real(dp), dimension(:, :), intent(in) :: DbM
204 :
205 : ! [mm] Total soil depth
206 : real(dp), dimension(:), intent(in) :: RZdepth
207 :
208 : ! mask at L0
209 : logical, dimension(:, :), intent(in) :: mask0
210 :
211 : ! Cell ids of hi res field
212 : integer(i4), dimension(:), intent(in) :: cell_id0
213 :
214 : ! Upper row of hi res block
215 : integer(i4), dimension(:), intent(in) :: upp_row_L1
216 :
217 : ! Lower row of hi res block
218 : integer(i4), dimension(:), intent(in) :: low_row_L1
219 :
220 : ! Left column of hi res block
221 : integer(i4), dimension(:), intent(in) :: lef_col_L1
222 :
223 : ! Right column of hi res block
224 : integer(i4), dimension(:), intent(in) :: rig_col_L1
225 :
226 : ! Number of L0 cells within a L1 cel
227 : integer(i4), dimension(:), intent(in) :: nL0_in_L1
228 :
229 : ! Parameter that determines the
230 : ! relative contribution to SM, upscaled
231 : ! Bulk density
232 : real(dp), dimension(:, :), intent(inout) :: L1_beta
233 :
234 : ! [10^-3 m] depth of saturated SM cont
235 : real(dp), dimension(:, :), intent(inout) :: L1_SMs
236 :
237 : ! [10^-3 m] field capacity
238 : real(dp), dimension(:, :), intent(inout) :: L1_FC
239 :
240 : ! [10^-3 m] permanent wilting point
241 : real(dp), dimension(:, :), intent(inout) :: L1_PW
242 :
243 : ! fraction of roots in soil horizons
244 : real(dp), dimension(:, :), intent(inout) :: L1_fRoots
245 :
246 : ! neutron count
247 : real(dp), dimension(:,:,:), intent(in) :: latWat_till ! lattice water
248 : real(dp), dimension(:,:,:), intent(in) :: COSMIC_L3_till! COSMIC parameter L3
249 : real(dp), dimension(:,:), intent(in) :: latWat ! lattice water
250 : real(dp), dimension(:,:), intent(in) :: COSMIC_L3 ! COSMIC paramter L3
251 : ! out
252 : real(dp), dimension(:,:), intent(inout) :: L1_bulkDens
253 : real(dp), dimension(:,:), intent(inout) :: L1_latticeWater
254 : real(dp), dimension(:,:), intent(inout) :: L1_COSMICL3
255 :
256 :
257 : ! loop index
258 : integer(i4) :: h
259 :
260 : ! loop index
261 : integer(i4) :: k
262 :
263 : ! loop index
264 : integer(i4) :: l
265 :
266 : ! loop index
267 : integer(i4) :: s
268 :
269 19167964 : real(dp) :: dpth_f
270 :
271 19167964 : real(dp) :: dpth_t
272 :
273 19167964 : real(dp) :: fTotRoots
274 :
275 : ! beta 0
276 19167964 : real(dp), dimension(size(LCOVER0, 1)) :: beta0
277 :
278 : ! [10^3 kg/m3] Bulk density
279 19168378 : real(dp), dimension(size(LCOVER0, 1)) :: Bd0
280 :
281 : ! [10^-3 m] depth of saturated SM cont
282 19167550 : real(dp), dimension(size(LCOVER0, 1)) :: SMs0
283 :
284 : ! [10^-3 m] field capacity
285 19167964 : real(dp), dimension(size(LCOVER0, 1)) :: FC0
286 :
287 : ! [10^-3 m] permanent wilting point
288 19168378 : real(dp), dimension(size(LCOVER0, 1)) :: PW0
289 :
290 : ! neutron count
291 19168378 : real(dp), dimension(size(LCOVER0,1)) :: LW0 ! lattice water
292 19168378 : real(dp), dimension(size(LCOVER0,1)) :: L30 ! COSMIC parameter L3
293 :
294 :
295 : ! fraction of roots in soil horizons
296 19168378 : real(dp), dimension(size(LCOVER0, 1)) :: fRoots0
297 :
298 19167964 : real(dp) :: tmp_rootFractionCoefficient_forest
299 :
300 19167964 : real(dp) :: tmp_rootFractionCoefficient_impervious
301 :
302 19167964 : real(dp) :: tmp_rootFractionCoefficient_pervious
303 :
304 : ! Field capacity dependent
305 : ! root frac coeffiecient
306 19167964 : real(dp) :: tmp_rootFractionCoefficient_perviousFC
307 :
308 : ! Model parameter describing the threshold for
309 : ! actual ET reduction for sand
310 19167964 : real(dp) :: tmp_rootFractionCoefficient_sand
311 :
312 : ! Model parameter describing the threshold for actual
313 : ! ET reduction for clay
314 19167964 : real(dp) :: tmp_rootFractionCoefficient_clay
315 :
316 19167964 : real(dp) :: FCmin_glob
317 :
318 : ! real(dp) :: FCdelta_glob
319 :
320 19167964 : real(dp) :: FCmax_glob
321 :
322 19167964 : real(dp) :: FCnorm
323 : ! the minimum number of till horizons
324 : integer(i4) :: min_nTH
325 :
326 :
327 611064 : min_nTH = minval(nTillHorizons(:))
328 414 : tmp_rootFractionCoefficient_forest = param(1) ! min(1.0_dp, param(2) + param(3) + param(1))
329 414 : tmp_rootFractionCoefficient_impervious = param(2)
330 :
331 : ! decide which parameterization should be used for route fraction:
332 414 : select case (processMatrix(3, 1))
333 : case(1,2)
334 : ! 1 and 2 - dependent on land cover
335 410 : tmp_rootFractionCoefficient_pervious = param(1) - param(3) ! min(1.0_dp, param(2) + param(3))
336 : !write(*,*) 'tmp_rootFractionCoefficient_forest = ', tmp_rootFractionCoefficient_forest
337 : !write(*,*) 'tmp_rootFractionCoefficient_impervious = ', tmp_rootFractionCoefficient_impervious
338 : !write(*,*) 'tmp_rootFractionCoefficient_pervious = ', tmp_rootFractionCoefficient_pervious
339 : case(3,4)
340 : ! 3 and 4 - dependent on land cover and additionally soil texture
341 4 : tmp_rootFractionCoefficient_pervious = param(3) ! min(1.0_dp, param(2) + param(3))
342 : !delta approach is used as in tmp_rootFractionCoefficient_pervious
343 4 : tmp_rootFractionCoefficient_sand = param(6) - param(5)
344 : !the value in parameter namelist is before substraction i.e. param(5)
345 4 : tmp_rootFractionCoefficient_clay = param(6)
346 4 : FCmin_glob=param(7)
347 418 : FCmax_glob=param(7)+param(8)
348 : !write(*,*) 'tmp_rootFractionCoefficient_forest = ', tmp_rootFractionCoefficient_forest
349 : !write(*,*) 'tmp_rootFractionCoefficient_impervious = ', tmp_rootFractionCoefficient_impervious
350 : !write(*,*) 'tmp_rootFractionCoefficient_pervious = ', tmp_rootFractionCoefficient_pervious
351 : !write(*,*) 'tmp_rootFractionCoefficient_sand = ', tmp_rootFractionCoefficient_sand
352 : !write(*,*) 'tmp_rootFractionCoefficient_clay = ', tmp_rootFractionCoefficient_clay
353 : !write(*,*) 'FCmin_glob = ', FCmin_glob
354 : !write(*,*) 'FCmax_glob = ', FCmax_glob
355 : end select
356 :
357 :
358 : ! select case according to a given soil database flag
359 828 : SELECT CASE(iFlag_soil)
360 : ! classical mHM soil database format
361 : CASE(0)
362 1242 : do h = 1, nHorizons_mHM
363 :
364 38335928 : Bd0 = nodata_dp
365 38335928 : SMs0 = nodata_dp
366 38335928 : FC0 = nodata_dp
367 38335928 : PW0 = nodata_dp
368 38335928 : fRoots0 = nodata_dp
369 38335928 : tmp_rootFractionCoefficient_perviousFC = nodata_dp
370 :
371 : ! neutron count
372 38335928 : LW0 = nodata_dp
373 38335928 : L30 = nodata_dp
374 :
375 : ! Initalise mHM horizon depth
376 : ! Last layer depth is soil type dependent, and hence it assigned within the inner loop
377 : ! by default for the first soil layer
378 828 : dpth_f = 0.0_dp
379 828 : dpth_t = HorizonDepth(H)
380 : ! check for the layer (2, ... n-1 layers) update depth
381 828 : if( (H .gt. 1) .and. (H .lt. nHorizons_mHM) ) then
382 0 : dpth_f = HorizonDepth(H-1)
383 : dpth_t = HorizonDepth(H)
384 : end if
385 :
386 : !$OMP PARALLEL
387 : !$OMP DO PRIVATE( l, s ) SCHEDULE( STATIC )
388 38335928 : cellloop0 : do k = 1, size(LCOVER0, 1)
389 38335100 : l = LCOVER0(k)
390 38335100 : s = soilID0(k, 1) ! >> in this case the second dimension of soilId0 = 1
391 : ! depth weightage bulk density
392 191675500 : Bd0(k) = sum(Db(s, : nTillHorizons(s), L) * Wd(S, H, 1 : nTillHorizons(S)), &
393 : Wd(S, H, 1 : nTillHorizons(S)) > 0.0_dp) &
394 38335100 : + sum(dbM(S, nTillHorizons(S) + 1 : nHorizons(S)) &
395 0 : * Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)), &
396 345015900 : Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)) >= 0.0_dp)
397 : ! depth weightage thetaS
398 76670200 : SMs0(k) = sum(thetaS_till(S, : nTillHorizons(s), L) &
399 38335100 : * Wd(S, H, 1 : nTillHorizons(S)), &
400 : Wd(S, H, 1 : nTillHorizons(S)) > 0.0_dp) &
401 0 : + sum(thetaS(S, nTillHorizons(S) + 1 - min_nTH : nHorizons(s) - min_nTH) &
402 0 : * Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)), &
403 230010600 : Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)) > 0.0_dp)
404 : ! depth weightage FC
405 76670200 : FC0(k) = sum(thetaFC_till(S, : nTillHorizons(s), L) &
406 38335100 : * Wd(S, H, 1 : nTillHorizons(S)), &
407 : Wd(S, H, 1 : nTillHorizons(S)) > 0.0_dp) &
408 0 : + sum(thetaFC(S, nTillHorizons(S) + 1 - min_nTH : nHorizons(s) - min_nTH) &
409 0 : * Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)), &
410 230010600 : Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)) > 0.0_dp)
411 : ! depth weightage PWP
412 76670200 : PW0(k) = sum(thetaPW_till(S, : nTillHorizons(s), L) &
413 38335100 : * Wd(S, H, 1 : nTillHorizons(S)), &
414 : Wd(S, H, 1 : nTillHorizons(S)) > 0.0_dp) &
415 0 : + sum(thetaPW(S, nTillHorizons(S) + 1 - min_nTH : nHorizons(s) - min_nTH) &
416 0 : * Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)), &
417 230010600 : Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)) > 0.0_dp)
418 :
419 : ! neutron count --> depth weightage LW and L30
420 76670200 : LW0(k) = sum( latWat_till(S, : nTillHorizons(s), L) &
421 38335100 : * Wd(S, H, 1 : nTillHorizons(S) ), &
422 : Wd(S, H, 1 : nTillHorizons(S)) > 0.0_dp ) &
423 0 : + sum( latWat(S,nTillHorizons(S) + 1 - min_nTH : nHorizons(s) - min_nTH) &
424 0 : * Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)), &
425 230010600 : Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)) > 0.0_dp )
426 :
427 76670200 : L30(k) = sum( COSMIC_L3_till(S, : nTillHorizons(s), L) &
428 38335100 : * Wd(S, H, 1 : nTillHorizons(S) ), &
429 : Wd(S, H, 1 : nTillHorizons(S)) > 0.0_dp ) &
430 0 : + sum( COSMIC_L3(S,nTillHorizons(S) + 1 - min_nTH : nHorizons(s) - min_nTH) &
431 0 : * Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)), &
432 230010600 : Wd(S, H, nTillHorizons(S) + 1 : nHorizons(S)) > 0.0_dp )
433 :
434 :
435 : ! Horizon depths: last soil horizon is varying, and thus the depth
436 : ! of the horizon too...
437 38335100 : if(H .eq. nHorizons_mHM) then
438 19167550 : dpth_f = HorizonDepth(nHorizons_mHM - 1)
439 19167550 : dpth_t = RZdepth(S)
440 : end if
441 : ! other soil properties [SMs, FC, PWP in mm]
442 38335100 : SMs0(k) = SMs0(k) * (dpth_t - dpth_f)
443 38335100 : FC0(k) = FC0(k) * (dpth_t - dpth_f)
444 38335100 : PW0(k) = PW0(k) * (dpth_t - dpth_f)
445 38335928 : LW0(k) = LW0(k) * (dpth_t - dpth_f)
446 : end do cellloop0
447 : !$OMP END DO
448 : !$OMP END PARALLEL
449 :
450 :
451 : !$OMP PARALLEL
452 : !$OMP DO PRIVATE( l, tmp_rootFractionCoefficient_perviousFC, FCnorm ) SCHEDULE( STATIC )
453 38335928 : celllloop0 : do k = 1, size(LCOVER0, 1)
454 38335100 : l = LCOVER0(k)
455 : !---------------------------------------------------------------------
456 : ! Effective root fractions in soil horizon...
457 : ! as weightage sum (according to LC fraction)
458 : !---------------------------------------------------------------------
459 : ! vertical root distribution = f(LC), following asymptotic equation
460 : ! [for refrence see, Jackson et. al., Oecologia, 1996. 108(389-411)]
461 :
462 : ! Roots(H) = 1 - beta^d
463 : ! where,
464 : ! Roots(H) = cumulative root fraction [-], range: 0-1
465 : ! beta = fitted extinction cofficient parameter [-], as a f(LC)
466 : ! d = soil surface to depth [cm]
467 :
468 : ! NOTES **
469 : ! sum(fRoots) for soil horions = 1.0
470 :
471 : ! if [sum(fRoots) for soil horions < 1.0], then
472 : ! normalise fRoot distribution such that all roots end up
473 : ! in soil horizon and thus satisfying the constrain that
474 : ! sum(fRoots) = 1
475 :
476 : ! The above constrains means that there are not roots below the soil horizon.
477 : ! This may or may not be realistic but it has been coded here to satisfy the
478 : ! conditions of the EVT vales, otherwise which the EVT values would be lesser
479 : ! than the acutal EVT from whole soil layers.
480 :
481 : ! Code could be modified in a way that a portion of EVT comes from the soil layer
482 : ! which is in between unsaturated and saturated zone or if necessary the saturated
483 : ! layer (i.e. Groundwater layer) can also contribute to EVT. Note that the above
484 : ! modification should be done only if and only if [sum(fRoots) for soil horions < 1.0].
485 : ! In such cases, you have to judiciously decide which layers (either soil layer between
486 : ! unsaturated and saturated zone or saturated zone) will contribute to EVT and in which
487 : ! proportions. Also note that there are no obervations on the depth avialable ata a
488 : ! moment on these layers.
489 : !------------------------------------------------------------------------
490 828 : select case(L)
491 : case(1)
492 : ! forest
493 13734852 : fRoots0(k) = (1.0_dp - tmp_rootFractionCoefficient_forest**(dpth_t * 0.1_dp)) &
494 13734852 : - (1.0_dp - tmp_rootFractionCoefficient_forest**(dpth_f * 0.1_dp))
495 : case(2)
496 : ! impervious
497 2412708 : fRoots0(k) = (1.0_dp - tmp_rootFractionCoefficient_impervious**(dpth_t * 0.1_dp)) &
498 2412708 : - (1.0_dp - tmp_rootFractionCoefficient_impervious**(dpth_f * 0.1_dp))
499 : case(3)
500 : ! permeable
501 22187540 : select case (processMatrix(3, 1))
502 : case(1,2)
503 : ! permeable
504 21971588 : fRoots0(k) = (1.0_dp - tmp_rootFractionCoefficient_pervious**(dpth_t * 0.1_dp)) &
505 21971588 : - (1.0_dp - tmp_rootFractionCoefficient_pervious**(dpth_f * 0.1_dp))
506 : case(3,4)
507 : ! permeable
508 : ! introducing global FC dependency on root frac. coef. by Simon Stisen and M. Cuneyd Demirel from GEUS.dk
509 : ! The normalization is based on Demirel et al 2018 (doi: 10.5194/hess-22-1299-2018)
510 : ! Case 3 is based on Jarvis (doi: 10.1016/0022-1694(89)90050-4)
511 : ! Case 4 is based on Feddes (doi: 10.1016/0022-1694(76)90017-2)
512 :
513 215952 : FCnorm = (((FC0(k) / (dpth_t - dpth_f)) - FCmin_glob) / (FCmax_glob - FCmin_glob))
514 :
515 215952 : if(FCnorm .lt. 0.0_dp) then
516 : ! print*, "FCnorm is below 0, will become 0", FCnorm
517 : FCnorm=0.0_dp
518 215216 : else if(FCnorm .gt. 1.0_dp) then
519 : ! print*, "FCnorm is above 1, will become 1", FCnorm
520 8804 : FCnorm=1.0_dp
521 : end if
522 :
523 : tmp_rootFractionCoefficient_perviousFC = (FCnorm * tmp_rootFractionCoefficient_clay) &
524 215952 : + ((1 - FCnorm) * tmp_rootFractionCoefficient_sand)
525 :
526 : fRoots0(k) = (1.0_dp - tmp_rootFractionCoefficient_perviousFC**(dpth_t * 0.1_dp)) &
527 215952 : - (1.0_dp - tmp_rootFractionCoefficient_perviousFC**(dpth_f * 0.1_dp))
528 :
529 : end select
530 :
531 22187540 : if((fRoots0(k) .lt. 0.0_dp) .OR. (fRoots0(k) .gt. 1.0_dp)) then
532 : ! why is this not stopping here?
533 : call message('***ERROR: Fraction of roots out of range [0,1]. Cell', &
534 0 : num2str(k), ' has value ', num2str(fRoots0(k)))
535 : ! stop
536 : end if
537 : end select
538 :
539 : end do celllloop0
540 : !$OMP END DO
541 : !$OMP END PARALLEL
542 :
543 38335928 : beta0 = Bd0 * param(4)
544 :
545 : !---------------------------------------------
546 : ! Upscale the soil related parameters
547 : !---------------------------------------------
548 828 : L1_SMs(:, h) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
549 32088 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, SMs0)
550 828 : L1_beta(:, h) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
551 32088 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, beta0)
552 828 : L1_PW(:, h) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
553 32088 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, PW0)
554 828 : L1_FC(:, h) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
555 32088 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, FC0)
556 828 : L1_fRoots(:, h) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
557 32088 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, fRoots0)
558 :
559 : ! !> neutron count
560 828 : L1_bulkDens(:,h) = upscale_harmonic_mean( nL0_in_L1, Upp_row_L1, Low_row_L1, &
561 32088 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, Bd0 )
562 828 : L1_latticeWater(:,h) = upscale_harmonic_mean( nL0_in_L1, Upp_row_L1, Low_row_L1, &
563 32088 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, LW0 )
564 828 : L1_COSMICL3(:,h) = upscale_harmonic_mean( nL0_in_L1, Upp_row_L1, Low_row_L1, &
565 32502 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, L30 )
566 :
567 : end do
568 :
569 : ! to handle multiple soil horizons with unique soil class
570 : CASE(1)
571 : ! horizon wise calculation
572 0 : do h = 1, nHorizons_mHM
573 0 : Bd0 = nodata_dp
574 0 : SMs0 = nodata_dp
575 0 : FC0 = nodata_dp
576 0 : PW0 = nodata_dp
577 0 : fRoots0 = nodata_dp
578 0 : tmp_rootFractionCoefficient_perviousFC = nodata_dp
579 :
580 : ! neutron count
581 0 : LW0 = nodata_dp
582 0 : L30 = nodata_dp
583 :
584 : ! initalise mHM horizon depth
585 0 : if (h .eq. 1) then
586 0 : dpth_f = 0.0_dp
587 0 : dpth_t = HorizonDepth(h)
588 : ! check for the layer (2, ... n-1 layers) update depth
589 : else
590 0 : dpth_f = HorizonDepth(h - 1)
591 0 : dpth_t = HorizonDepth(h)
592 : end if
593 : ! need to be done for every layer to get fRoots
594 : !$OMP PARALLEL
595 : !$OMP DO PRIVATE( l, s ) SCHEDULE( STATIC )
596 0 : cellloop1 : do k = 1, size(LCOVER0, 1)
597 0 : L = LCOVER0(k)
598 0 : s = soilID0(k, h)
599 0 : if (h .le. nTillHorizons(1)) then
600 0 : Bd0(k) = Db(s, 1, L)
601 0 : SMs0(k)= thetaS_till (s, 1, L) * (dpth_t - dpth_f) ! in mm
602 0 : FC0(k) = thetaFC_till(s, 1, L) * (dpth_t - dpth_f) ! in mm
603 0 : PW0(k) = thetaPW_till(s, 1, L) * (dpth_t - dpth_f) ! in mm
604 0 : LW0(k) = latWat_till(s, 1, L) * (dpth_t - dpth_f) ! in mm ! >> neutron count
605 : else
606 0 : Bd0(k) = DbM(s, 1)
607 0 : SMs0(k)= thetaS (s, 1) * (dpth_t - dpth_f) ! in mm
608 0 : FC0(k) = thetaFC(s, 1) * (dpth_t - dpth_f) ! in mm
609 0 : PW0(k) = thetaPW(s, 1) * (dpth_t - dpth_f) ! in mm
610 0 : LW0(k) = latWat(s, 1) * (dpth_t - dpth_f) ! in mm ! >> neutron count
611 : end if
612 : end do cellloop1
613 : !$OMP END DO
614 : !$OMP END PARALLEL
615 :
616 :
617 : !$OMP PARALLEL
618 : !$OMP DO PRIVATE( l, tmp_rootFractionCoefficient_perviousFC, FCnorm ) SCHEDULE( STATIC )
619 :
620 0 : celllloop1 : do k = 1, size(LCOVER0, 1)
621 0 : l = LCOVER0(k)
622 : !================================================================================
623 : ! fRoots = f[LC] --> (fRoots(H) = 1 - beta^d)
624 : ! see below for comments and references for the use of this simple equation
625 : ! NOTE that in this equation the unit of soil depth is in cm
626 : !================================================================================
627 :
628 0 : select case(L)
629 : case(1)
630 : ! forest
631 0 : fRoots0(k) = (1.0_dp - tmp_rootFractionCoefficient_forest**(dpth_t * 0.1_dp)) &
632 0 : - (1.0_dp - tmp_rootFractionCoefficient_forest**(dpth_f * 0.1_dp))
633 : case(2)
634 : ! impervious
635 0 : fRoots0(k) = (1.0_dp - tmp_rootFractionCoefficient_impervious**(dpth_t * 0.1_dp)) &
636 0 : - (1.0_dp - tmp_rootFractionCoefficient_impervious**(dpth_f * 0.1_dp))
637 : case(3)
638 :
639 0 : select case (processMatrix(3, 1))
640 : case(1,2)
641 : ! permeable
642 0 : fRoots0(k) = (1.0_dp - tmp_rootFractionCoefficient_pervious**(dpth_t * 0.1_dp)) &
643 0 : - (1.0_dp - tmp_rootFractionCoefficient_pervious**(dpth_f * 0.1_dp))
644 : case(3,4)
645 : ! permeable
646 : ! introducing global FC dependency on root frac. coef. by Simon Stisen and M. Cuneyd Demirel from GEUS.dk
647 : ! The normalization is based on Demirel et al 2018 (doi: 10.5194/hess-22-1299-2018)
648 : ! Case 3 is based on Jarvis (doi: 10.1016/0022-1694(89)90050-4)
649 : ! Case 4 is based on Feddes (doi: 10.1016/0022-1694(76)90017-2)
650 0 : FCnorm = (((FC0(k) / (dpth_t - dpth_f)) - FCmin_glob) / (FCmax_glob - FCmin_glob))
651 :
652 0 : if(FCnorm .lt. 0.0_dp) then
653 : ! print*, "FCnorm is below 0, will become 0", FCnorm
654 : FCnorm=0.0_dp
655 0 : else if(FCnorm .gt. 1.0_dp) then
656 : ! print*, "FCnorm is above 1, will become 1", FCnorm
657 0 : FCnorm=1.0_dp
658 : end if
659 :
660 : tmp_rootFractionCoefficient_perviousFC = (FCnorm * tmp_rootFractionCoefficient_clay) &
661 0 : + ((1 - FCnorm) * tmp_rootFractionCoefficient_sand)
662 :
663 :
664 : fRoots0(k) = (1.0_dp - tmp_rootFractionCoefficient_perviousFC**(dpth_t * 0.1_dp)) &
665 0 : - (1.0_dp - tmp_rootFractionCoefficient_perviousFC**(dpth_f * 0.1_dp))
666 :
667 : end select
668 :
669 0 : if((fRoots0(k) .lt. 0.0_dp) .OR. (fRoots0(k) .gt. 1.0_dp)) then
670 : ! why is this not stopping here?
671 : call message('***ERROR: Fraction of roots out of range [0,1]. Cell', &
672 0 : num2str(k), ' has value ', num2str(fRoots0(k)))
673 : ! stop
674 : end if
675 : end select
676 :
677 : end do celllloop1
678 : !$OMP END DO
679 : !$OMP END PARALLEL
680 :
681 : ! beta parameter
682 0 : beta0 = Bd0 * param(4)
683 :
684 : ! Upscale the soil related parameters
685 0 : L1_SMs(:, h) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
686 0 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, SMs0)
687 0 : L1_beta(:, h) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
688 0 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, beta0)
689 0 : L1_PW(:, h) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
690 0 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, PW0)
691 0 : L1_FC(:, h) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
692 0 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, FC0)
693 0 : L1_fRoots(:, h) = upscale_harmonic_mean(nL0_in_L1, Upp_row_L1, Low_row_L1, &
694 0 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, fRoots0)
695 :
696 : ! neutron count
697 0 : L1_bulkDens(:,h) = upscale_harmonic_mean( nL0_in_L1, Upp_row_L1, Low_row_L1, &
698 0 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, Bd0 )
699 0 : L1_latticeWater(:,h) = upscale_harmonic_mean( nL0_in_L1, Upp_row_L1, Low_row_L1, &
700 0 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, LW0 )
701 0 : L1_COSMICL3(:,h) = upscale_harmonic_mean( nL0_in_L1, Upp_row_L1, Low_row_L1, &
702 0 : Lef_col_L1, Rig_col_L1, cell_id0, mask0, nodata_dp, L30 )
703 :
704 : end do
705 : ! anything else
706 : CASE DEFAULT
707 414 : call error_message('***ERROR: iFlag_soilDB option given does not exist. Only 0 and 1 is taken at the moment.')
708 : END SELECT
709 :
710 :
711 : ! below operations are common to all soil databases flags
712 : !$OMP PARALLEL
713 : !------------------------------------------------------------------------
714 : ! CHECK LIMITS OF PARAMETERS
715 : ! [PW <= FC <= ThetaS]
716 : ! units of all variables are now in [mm]
717 : ! If by any means voilation of this rule appears (e.g. numerical errors)
718 : ! than correct it --> threshold limit = 1% of the upper ones
719 : !------------------------------------------------------------------------
720 31674 : L1_FC = merge(L1_SMs - 0.01_dp * L1_SMs, L1_FC, L1_FC .gt. L1_SMs)
721 31674 : L1_PW = merge(L1_FC - 0.01_dp * L1_FC, L1_PW, L1_PW .gt. L1_FC )
722 : ! check the physical limit
723 31674 : L1_SMs = merge(0.0001_dp, L1_SMs, L1_SMs .lt. 0.0_dp)
724 31674 : L1_FC = merge(0.0001_dp, L1_FC, L1_FC .lt. 0.0_dp)
725 31674 : L1_PW = merge(0.0001_dp, L1_PW, L1_PW .lt. 0.0_dp)
726 : ! Normalise the vertical root distribution profile such that [sum(fRoots) = 1.0]
727 : !$OMP DO PRIVATE( fTotRoots ) SCHEDULE( STATIC )
728 15630 : do k = 1, size(L1_fRoots, 1)
729 45648 : fTotRoots = sum(L1_fRoots(k, :), L1_fRoots(k, :) .gt. 0.0_dp)
730 : ! This if clause is necessary for test program but may be redundant in actual program
731 15630 : if (fTotRoots .gt. 0.0_dp) then
732 45648 : L1_fRoots(k, :) = L1_fRoots(k, :) / fTotRoots
733 : else
734 0 : L1_fRoots(k, :) = 0.0_dp
735 : end If
736 : end do
737 :
738 : !$OMP END DO
739 : !$OMP END PARALLEL
740 : !close(1)
741 414 : end subroutine mpr_SMhorizons
742 :
743 : end module mo_mpr_SMhorizons
|
