-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathdrains.f
222 lines (205 loc) · 9.22 KB
/
drains.f
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
subroutine drains
!! ~ ~ ~ PURPOSE ~ ~ ~
!! this subroutine finds the effective lateral hydraulic conductivity
!! and computes drainage or subirrigation flux
!! ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! conk(:,:) |mm/hr |lateral saturated hydraulic conductivity for each profile
!! |layer in a give HRU. For example (conk(2,1) is conductivity
!! |of layer from sol_z(1,1) to sol_z(2,1) in HRU1
!! curyr |none |current year in simulation (sequence)
!! drain_co(:) |mm/day |drainage coefficient
!! dg(:,:) |mm |depth of soil layer
!! ddrain(:) |mm |depth of drain tube from the soil surface
!! hru_slp(:) |m/m |average slope steepness in HRU
!! id1 |julian date |first day of simulation in current year
!! ihru |none |HRU number
!! latksatf(:) |none |multiplication factor to determine conk(j1,j) from sol_k(j1,j) for HRU
!! pc(:) |mm/hr |pump capacity (default pump capacity = 1.042mm/hr or 25mm/day)
!! sdrain(:) |mm |distance between two drain tubes or tiles
!! sstmaxd(:) |mm |static maximum depressional storage; read from .sdr
!! sol_k(:,:) |mm/hr |saturated hydraulic conductivity of soil
!! |layer
!! sol_nly(:) |none |number of layers in soil profile
!! sol_z(:,:) |mm |depth to bottom of each profile layer in a given HRU
!! stmaxd(:) |mm |maximum surface depressional storage for the day in a given HRU
!! stor |mm |surface storage for the day in a given HRU
!! storro |mm |surface storage that must b
!! |can move to the tile drain tube
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! lyrtile |mm H2O |drainage tile flow in soil layer for day
!! qtile |mm H2O |drainage tile flow in soil profile for the day
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! above |mm |depth of top layer considered
!! depth |mm |actual depth from surface to impermeable layer
!! cone |mm/hr |effective saturated lateral conductivity - based
!! |on water table depth and conk/sol_k of layers
!! ddarnp |mm |a variable used to indicate distance slightly less
!! |than ddrain. Used to prevent calculating subirrigation
!! |when water table is below drain bottom or when it is empty
!! deep |mm |total thickness of saturated zone
!! depth |mm |effective depth to impermeable layer from soil surface
!! |effective depth may be smaller than actual depth to account
!! |for convergence near drain tubes
!! dflux |mm/hr |drainage flux
!! dot |mm |actual depth from impermeable layer to water level
!! |above drain during subsurface irrigation
!! em |mm |distance from water level in the drains to water table
!! |at midpoint: em is negative during subirrigation
!! gee |none |factor -g- in Kirkham equation
!! hdrain |mm |equivalent depth from water surface in drain tube to
!! |impermeable layer
!! i |none |counter
!! j |none |HRU number
!! j1 |none |counter
!! w |mm |thickness of saturated zone in layer considered
!! y1 |mm |dummy variable for dtwt
!! nlayer |none |number of layers to be used to determine cone
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!! Intrinsic:
!! SWAT: depstor
!! ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~
use parm
integer :: j1, j, m
real:: cone, depth, dg, ad, ap
real:: hdrain, gee, e, gee1, gee2, gee3, pi
real:: k2, k3, k4, k5, k6
!! initialize variables
j = 0
j = ihru
wnan = 0
y1 = dep_imp(j) - wt_shall
if (y1 > dep_imp(j)) y1 = dep_imp(j)
above = 0.
pi = 22./7.
gee1 =0.
!! find number of soil layers
do j1 = 1, mlyr
if(sol_z(j1,j) > 0.) nlayer = j1
end do
!! find effective lateral hydraulic conductivity for the profile in hru j
do j1 = 1, nlayer
if(y1 > sol_z(j1,j)) then
wnan(j1) = 0.
else
wnan(j1) = sol_z(j1,j) - y1
x = sol_z(j1,j) - above
if(wnan(j1) > x) wnan(j1) = x
end if
above = sol_z(j1,j)
end do
sum = 0.
deep = 0.
do j1=1,nlayer
conk(j1,j) = sol_k(j1,j) * latksatf(j) !Daniel 2/26/09
sum = sum + wnan(j1) * conk(j1,j)
deep = deep + wnan(j1)
end do
if((deep <= 0.001).or.(sum <= 0.001)) then
sum = 0.
deep = 0.001
do j1=1,nlayer
!! Compute layer depth ! Daniel 10/05/07
dg = 0.
if(j1 == 1) then
dg = sol_z(j1,j)
else
dg = sol_z(j1,j) - sol_z(j1-1,j)
end if
!! Compute layer depth ! Daniel 10/05/07
sum=sum+conk(j1,j)*dg !Daniel 10/09/07
deep=deep+dg !Daniel 10/09/07
end do
cone=sum/deep
else
cone=sum/deep
end if
!! calculate parameters hdrain and gee1
ad = dep_imp(j) - ddrain(j)
ap = 3.55 - ((1.6 * ad) / sdrain(j)) + 2 * ((2 / sdrain(j))**2)
if (ad/sdrain(j) < 0.3) then
hdrain= ad / (1 + ((ad / sdrain(j)) * (((8 / pi) *
& Log(ad / re(j)) - ap))))
else
hdrain = (sdrain(j) * pi) / (8 * ((log(sdrain(j) / re(j))/
& log(e)) - 1.15))
end if
!! calculate Kirkham G-Factor, gee
k2 = tan((pi * ((2. * ad) - re(j))) / (4. * dep_imp(j)))
k3 = tan((pi * re(j)) / (4. * dep_imp(j)))
do m=1,2
k4 = (pi * m * sdrain(j)) / (2. * dep_imp(j))
k5 = (pi * re(j)) / (2. * dep_imp(j))
k6 = (pi * (2. * ad - re(j))) / (2. * dep_imp(j))
gee2 = (cosh(k4) + cos(k5)) / (cosh(k4) - cos(k5))
gee3 = (cosh(k4) - cos(k6)) / (cosh(k4) + cos(k6))
gee1 = gee1 + Log(gee2 * gee3)
end do
xx = k2 / k3
if (xx < 1.) then
gee = 1.
else
gee = 2 * Log(k2 / k3) + 2 * gee1
end if
if (gee < 1.) gee = 1.
if (gee > 12.) gee = 12.
!! calculate drainage and subirrigation flux section
! drainage flux for ponded surface
depth = ddrain(j) + hdrain
hdmin = depth - ddrain(j)
if (ismax == 1) then
call depstor ! dynamic stmaxd(j): compute current HRU stmaxd based
! on cumulative rainfall and cum. intensity
else
stmaxd(j) = sstmaxd(j)
end if
storro = 0.2 * stmaxd(j) !surface storage that must be filled before surface
!water can move to the tile drain tube
!! Determine surface storage for the day in a given HRU (stor)
!initialize stor on the beginning day of simulation, Daniel 9/20/2007
if (curyr == 1 .and. iida == id1) then
stor= 0.
end if
if (potsa(j) <= 0.) then ! determine stor
stor = precipday - inflpcp - etday !Daniel 10/05/07
if(surfq(j) > 0.0) stor=stmaxd(j)
else
stor = pot_vol(j)/(potsa(j)*1000)
endif
if(hdrain < hdmin) hdrain=hdmin
if((stor > storro).and.(y1 < 5.0)) then
dflux= (12.56637*24.0*cone*(depth-hdrain+stor))/
& (gee*sdrain(j)) !eq.10
if(dflux > drain_co(j)) dflux = drain_co(j) !eq.11
else
! subirrigation flux
em=depth-y1-hdrain
if(em < -1.0) then
ddranp=ddrain(j)-1.0
dot=hdrain+dep_imp(j)-depth
dflux=4.0*24.0*cone*em*hdrain*(2.0+em/dot)/sdrain(j)**2
if((depth-hdrain) >= ddranp) dflux=0.
if(abs(dflux) > pc(j)) then
dflux = -pc(j)*24.0
end if
! drainage flux - for WT below the surface and for ponded depths < storro (S1)
else
dflux=4.0*24.0*cone*em*(2.0*hdrain+em)/sdrain(j)**2 !eq.5
if(dflux > drain_co(j)) dflux=drain_co(j) !eq.11
if(dflux < 0.) dflux=0.
if(em < 0.) dflux=0.
end if
end if
qtile=dflux
! write(222,222) curyr, iida, hdrain, gee1, gee !Daniel 3/1/09
!222 format(1x,4x,i4,4x,i3,4x,3f12.3)
return
end