-
Notifications
You must be signed in to change notification settings - Fork 3
/
bmp_ri_pond.f
138 lines (108 loc) · 5.05 KB
/
bmp_ri_pond.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
subroutine ri_pond(kk,riflw,rised)
!! ~ ~ ~ PURPOSE ~ ~ ~
!! this subroutine routes water through a retention irrigation pond in the subbasin
!! ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! flw(:,:) |mm |stormwater runoff coming in/out of pond at a time step
!! kk |none |pond id number in the subbasin
!! ri_sed(:) |tons |total sediment deposited in the pond
!! sed(:,:) |mm |overland flow sediment coming in/out of pond at a time step
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! flw(:,:) |mm |stormwater runoff coming in/out of pond at a time step
!! ri_sed(:) |tons |total sediment deposited in the pond
!! sed(:,:) |mm |overland flow sediment coming in/out of pond at a time step
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! sb |none |subbasin or reach number
!! ii |none |time step
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!! Intrinsic: Min
!! ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~
use parm
implicit none
integer :: sb, kk, ii
real :: tsa,mxvol,pdia,ksat,dp,sub_ha,mxh,hweir,phead,pipeflow
real :: qin,qout,qpnd,sweir,hpnd,qet
real :: qweir, qseep,qpipe,qpndi,decayexp,splw,qpump
real :: sedconc,sedpndi, sedpnde,ksed,td,sedpump
real, dimension(4,0:nstep), intent(in out) :: riflw,rised
real, dimension(0:nstep) :: inflw,insed,outflw,outsed
sb = inum1
sub_ha = da_ha * sub_fr(sb)
qin = 0.; qout = 0.
outflw = 0.; outsed = 0.
!! Initialize parameters, coefficients, etc
tsa = ri_sa(sb,kk) !total surface area of pond (m^2)
mxvol = ri_vol(sb,kk) !max. capacity of the basin (m^3)
mxh = ri_dep(sb,kk) !max. depth of water, m
ksat = ri_k(sb,kk) !saturated hydraulic conductivity (mm/hr)
!! Get initial values from previous day
qpnd = ri_qi(sb,kk) !m^3
sedpnde = ri_sedi(sb,kk)
inflw(:) = riflw(1,:)
insed(:) = rised(1,:)
do ii=1,nstep
qout = 0.
!inflow = runoff + precipitation
qin = inflw(ii) * 10. * (sub_ha - tsa / 10000.) +
& precipdt(ii) * tsa / 1000. !m^3
!update ponded water volume
qpnd = qpnd + qin
!bypass flow when pond is full
if (qpnd>mxvol) then
qout = qpnd - mxvol !m3
outflw(ii) = qout
outsed(ii) = insed(ii) * qout / qin !tons
qpnd = mxvol
end if
!initial sediment
sedpndi = sedpnde + insed(ii)
if (qpnd>0) then
sedconc = sedpndi / qpnd * 1.e6 !mg/l
else
sedconc = 0.
endif
!Transmission loss through infiltration
qseep = ksat * tsa / 1000./ 60. * idt !m^3
bmp_recharge(sb) = bmp_recharge(sb)
& + qseep / (sub_ha*10000.- tsa) *1000.
!Evapotranspiration loss
qet = ri_evrsv(sb,kk) * tsa * pet_day / 1000. / 1440. * idt !m^3
!water pumped for irrigation
qpump = max(0.,ri_pmpvol(kk,ii))
!mass balance
qpnd = qpnd - qseep - qet - qpump
if(qpnd<0) qpnd = 0
hpnd = qpnd / (mxvol / mxh)
!Estimate TSS removal due to sedimentation
if (sedconc>12.) then ! assume 12mg/l as equilibrium concentration, , Huber et al. 2006
ksed = min(134.8,41.1 * hpnd ** -0.999) !decay coefficient, Huber et al. 2006
td = 1. / nstep !detention time, day
sedconc = (sedconc - 12.) * exp(-ksed * td) + 12.
endif
!sediment pumped
sedpump = qpump * sedconc / 1.e6 !tons
!sediment deposition
sedpnde = qpnd *sedconc / 1.e6 !tons
ri_sed_cumul(sb,kk) = ri_sed(sb,kk) + sedpndi - sedpnde
riflw(1,ii) = qin / (sub_ha *10000. - tsa) * 1000. !mm
riflw(2,ii) = outflw(ii) / (sub_ha *10000. - tsa) * 1000.
riflw(3,ii) = qpump / (sub_ha *10000. - tsa) * 1000.
riflw(4,ii) = qseep / (sub_ha *10000. - tsa) * 1000.
rised(3,:) = sedpump
rised(2,:) = outsed(:)
End do
! Store end-of-day values for next day
ri_qi(sb,kk) = qpnd
ri_qloss(kk,1) = qet
ri_qloss(kk,2) = qseep
ri_sedi(sb,kk) = sedpnde
return
end subroutine