forked from ARVE-Research/makesoil
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpedotransfermod.f90
297 lines (205 loc) · 8.55 KB
/
pedotransfermod.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
module pedotransfermod
! module containing the pedotransfer functions of Balland et al. (2008) as adjusted and modified by Sandoval et al. (2024)
use parametersmod, only : sp
implicit none
public :: fDp
public :: fDb
public :: fTsat
public :: fT33
public :: fT1500
public :: fKsat
! SoilGrids map legend for USDA soil orders
! Histosols: 5,10-13
! Spodosols: 15-19
! Inceptisols: 85,86,98-94
! Entisols: 95-99
! Aridisols: 50-56
! Mollisols: 69,70,72,73,74,77
! Aquoll, argiaquoll, argiudoll, cryaquoll, duraquoll, natrustoll and paleustoll 71,75,76
! Alfisols: 80-84
! Oxisols: 30-34
! Ultisols: 60-64
! Vertisols: 40-45
! Andisols: 20-27
! Typical bulk densities of coarse fragments
! Sedimentary rocks Siliceous Mafic Peridotite
! 1.8 2.0 2.2 2.4 2.6 2.7 2.8 3.0 3.2
real(sp), parameter :: Dcf = 2.7 ! bulk density of coarse fragments (g cm-3)
contains
! ----------------------------
real(sp) function fDp(orgm,cfvo)
! function to estimate soil particle density (g cm-3) (Sandoval et al., 2024, eqn A26)
use parametersmod, only : sp
implicit none
! argument
real(sp), intent(in) :: orgm ! soil organic matter (mass fraction)
real(sp), intent(in) :: cfvo ! coarse fragment content (volume fraction)
! local variable
real(sp) :: Dp
! ---
Dp = 1. / (orgm / 1.3 + (1. - orgm) / 2.65)
fDp = Dp * (1. - cfvo) + cfvo * Dcf
end function fDp
! ----------------------------
real(sp) function fDb(usda,clay,cfvo,zpos,orgm,Dp)
! function to estimate soil bulk density (Balland et al., 2008, eqn 18) in g cm-3
! This function requires information about the USDA soil type to determine the parameter set to use,
! e.g., the field SoilGrids250m 2017-03 - Predicted USDA 2014 suborder classes
! NB not all of the soil suborders listed in Table 6 of Balland et al. (2008)
! are present in the SoilGrids250m dataset.
use parametersmod, only : sp
implicit none
! arguments
integer, intent(in) :: usda ! USDA 2014 suborder class (code, complete legend in SoilGrids TAXOUSDA file)
real(sp), intent(in) :: clay ! clay content (mass fraction)
real(sp), intent(in) :: cfvo ! coarse fragment content (volume fraction)
real(sp), intent(in) :: zpos ! below surface (cm)
real(sp), intent(in) :: orgm ! organic matter content (mass fraction)
real(sp), intent(in) :: Dp ! particle density (g cm-3)
! Parameters specific to USDA soil order and suborder
real(sp), dimension(4,19), parameter :: parstable = reshape( &
[ 1.18, 0.89, 0.022, 6.27, & ! 1 Histosols
1.18, 0.70, 0.022, 6.27, & ! 2 Spodosols
1.18, 0.89, 0.022, 6.27, & ! 3 Inceptisols General
1.50, 1.10, 0.010, 6.27, & ! 4 Inceptisols Exceptions: dystropept, eutropept, halaquept, tropept, ustropept and xerumbrept
1.50, 1.10, 0.022, 6.27, & ! 5 Entisols General
1.50, 1.10, 0.002, 6.27, & ! 6 Entisols Exceptions: cryorthent
1.50, 1.10, 0.022, 6.27, & ! 7 Aridisols
1.50, 1.10, 0.022, 6.27, & ! 8 Mollisols General
1.50, 1.40, 0.022, 6.27, & ! 9 Mollisols Exceptions: aquoll, argiaquoll, argiudoll, cryaquoll, duraquoll, natrustoll and paleustoll
1.65, 1.10, 0.022, 6.27, & ! 10 Alfisols General
1.18, 0.89, 0.022, 6.27, & ! 11 Alfisols Exceptions: cryoboralf, fragiboralf, kandiudalf, kanhaplustalf and paleudalf
1.40, 1.20, 0.002, 6.27, & ! 12 Oxisols General
1.40, 1.20, 0.001, 6.27, & ! 13 Oxisols Exceptions: kandiudox
1.55, 1.00, 0.005, 6.27, & ! 14 Ultisols General
1.55, 1.00, 0.001, 6.27, & ! 15 Ultisols Exceptions: haploxerult, kandihumult, kanhaplustult, palehumult, rhodudult and udult
1.90, 1.20, 0.010, 6.27, & ! 16 Vertisols
1.40, 1.10, 0.050, 6.27, & ! 17 Andisols General
1.00, 0.00, 0.001, 2., & ! 18 Andisols Exceptions: haplaquands, hapludands, haplustands, melanudands, udivitrands, vitraquands, vitricryands, vitritorrands, vitrixerands
1.00, 0.00, 0.001, 2. ], & ! 19 Andisols Exceptions: melanudands, udivitrands, vitraquands, vitricryands, vitritorrands, vitrixerands
shape(parstable))
! local variables
real(sp), dimension(4) :: pars
real(sp) :: a
real(sp) :: b
real(sp) :: c
real(sp) :: d
real(sp) :: Dbs
! ----
select case(usda)
case(5,10:13) ! Histosols, including Gelisols: Histels
pars = parstable(:,1)
case(15:19) ! Spodosols
pars = parstable(:,2)
case(85:86,89:94) ! Inceptisols general
pars = parstable(:,3)
case(95:99) ! Entisols general
pars = parstable(:,5)
case(6:7) ! Gelisols (Orthels and Turbels), using table definition for Cryorthent
pars = parstable(:,6)
case(50:56) ! Aridisols
pars = parstable(:,7)
case(69:70,72:74,77) ! Mollisols general
pars = parstable(:,8)
case(71,75:76) ! Mollisols exceptions
pars = parstable(:,9)
case(80:84) ! Alfisols general
pars = parstable(:,10)
case(30:34) ! Oxisols general
pars = parstable(:,12)
case(60:64) ! Ultisols general
pars = parstable(:,14)
case(40:45) ! Vertisols
pars = parstable(:,16)
case(20:27) ! Andisols general
pars = parstable(:,17)
case default
write(0,*)'USDA type',usda,' was not classified, using average parameter set'
pars = [ 1.40, 0.95, 0.015, 5.82 ]
end select
a = pars(1)
b = pars(2)
c = pars(3)
d = pars(4)
Dbs = (a + (Dp - a - b * (1. - clay)) * (1. - exp(-c * zpos))) / (1. + d * orgm)
fDb = Dbs * (1. - cfvo) + cfvo * Dcf ! Balland et al. (2008) eqn A.13
end function fDb
! ----------------------------
real(sp) function fTsat(Dp,Db)
! function to estimate soil porosity Theta-sat (Sandoval et al., 2024, eqn A25c) (units)
use parametersmod, only : sp
implicit none
real(sp), intent(in) :: Dp ! particle density (g cm-3)
real(sp), intent(in) :: Db ! bulk density (g cm-3)
fTsat = 1. - Db / Dp
end function fTsat
! ----------------------------
real(sp) function fT33(Tsat,clay,sand,orgm)
! function to estimate soil water content at field capacity Theta-33 (Sandoval et al., 2024, eqn A25b) (units)
use parametersmod, only : sp
implicit none
! arguments
real(sp), intent(in) :: Tsat ! soil porosity (units)
real(sp), intent(in) :: clay ! clay content (mass fraction)
real(sp), intent(in) :: sand ! sand content (mass fraction)
real(sp), intent(in) :: orgm ! soil organic matter content (mass fraction)
! parameters
real(sp), parameter :: a = -0.0547
real(sp), parameter :: b = -0.0010
real(sp), parameter :: c = 0.4760
real(sp), parameter :: d = 0.9402
! ----
fT33 = Tsat * (c + (d - c) * clay**0.5) * exp((a * sand - b * orgm) / Tsat)
end function fT33
! ----------------------------
real(sp) function fT1500(T33,clay)
! function to estimate soil water content at wilting point Theta-1500 (Sandoval et al., 2024, eqn A25a) (units)
use parametersmod, only : sp
implicit none
! arguments
real(sp), intent(in) :: T33 ! field capacity water content (volume fraction)
real(sp), intent(in) :: clay ! clay content (mass fraction)
! parameters
real(sp), parameter :: c = 0.2018
real(sp), parameter :: d = 0.7809
! ----
fT1500 = T33 * (c + (d - c) * clay**0.5)
end function fT1500
! ----------------------------
real(sp) function fKsat(sand,clay,orgm,Db,Tsat,T33,T1500)
! function to estimate saturated hydraulic conductivity (Sandoval et al., 2024) (units)
! NB this equation comes from the code on github in the file splash.point.R, lines 351-363
! because eqn A25d in the GMD paper does not appear to produce valid results.
use parametersmod, only : sp
implicit none
! arguments
real(sp), intent(in) :: sand ! sand content (mass fraction)
real(sp), intent(in) :: clay ! clay content (mass fraction)
real(sp), intent(in) :: orgm ! sand content (mass fraction)
real(sp), intent(in) :: Db ! bulk density (g cm-3)
real(sp), intent(in) :: Tsat ! sand content (mass fraction)
real(sp), intent(in) :: T33 ! sand content (mass fraction)
real(sp), intent(in) :: T1500 ! sand content (mass fraction)
! parameters
real(sp), parameter :: Ksmax = 857.48454
real(sp), parameter :: k2 = -2.70927
real(sp), parameter :: k3 = 3.62264
real(sp), parameter :: k4 = 7.33398
real(sp), parameter :: k5 = -8.11795
real(sp), parameter :: k6 = 18.75552
real(sp), parameter :: k7 = 1.03319
real(sp), parameter :: l1500 = log(1500.)
real(sp), parameter :: l33 = log(33.)
real(sp), parameter :: num = l1500 - l33
! local variables
real(sp) :: Tdrain
real(sp) :: B
real(sp) :: lambda
! ----
Tdrain = Tsat - T33
B = num / (log(T33) - log(T1500))
lambda = 1. / B
fKsat = Ksmax / (1. + exp(k2 * sand + k3 * Db + k4 * clay + k5 * Tdrain + k6 * orgm + k7 * lambda))
end function fKsat
! ----------------------------
end module pedotransfermod