-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbmp_det_pond.f
212 lines (170 loc) · 8.7 KB
/
bmp_det_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
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
subroutine det_pond
!! ~ ~ ~ PURPOSE ~ ~ ~
!! the purpose of this program is to read in data from the detention pond
!! input file (.dtp) and perform computations
!! ~ ~ ~ INCOMING VARIABLES ~ ~ ~ ~ ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! dtp_evrsv |none |detention pond evaporation coefficient
!! dtp_weirtype(:)|none |Type of weir: 1 rectangular and 2 circular
!! dtp_numweir(:) |none |Total number of weirs in the BMP
!! dtp_numstage(:)|none |Total number of stages in the weir
!! dtp_wdratio(:,:)|none |Width depth ratio (rectangular wier) at different stages
!! dtp_depweir(:,:)|m |Depth of rectangular wier at different stages
!! dtp_diaweir(:,:)|m |Diameter of circular wier at different stages
!! dtp_retperd(:,:)|years |Return period at different stages
!! dtp_pcpret(:,:)|mm |precipitation for different return periods (not used)
!! dtp_cdis(:,:) |none |coeffieicne of discharge at different stages
!! hhvaroute(2,:,:) |m^3 H2O |water
!! hhvaroute(3,:,:) |metric tons |sediment or suspended solid load
!! i_mo |none |current month of simulation
!! sub_subp_dt(:,:) |mm H2O |precipitation for time step in subbasin
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ OUTGOING VARIABLES ~ ~ ~ ~ ~ ~ ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! hhvaroute(2,:,:) |m^3 H2O |water
!! hhvaroute(3,:,:) |metric tons|sediment or suspended solid load
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!! name |units |definition
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ii |none |time step counter
!! k |none |weir stage counter
!! titldum |NA |dummy string
!! ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!! ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!! SWAT: surf_seep,surf_evap,water_depth
!! Intrinsic: Log,exp
!! ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~
use parm
implicit none
character (len=80) :: titldum
integer :: ii, k, sb
real :: qin,qout,qpnd,sedin,sedout,sedpnd,spndconc,qdepth,
& watdepact,qstage,backup_length,seep_sa,evap_sa,pcp_vol,
& evap_vol,seep_vol,warea,pi,qovmax,qaddon,depaddon
pi = 3.14159
sb = inum1
qout = 0.; sedout = 0.; depaddon = 0.
if (iyr<dtp_iyr(sb) .or.
&(iyr==dtp_iyr(sb) .and. i_mo<dtp_imo(sb))) then
return
endif
!! Get initial values from previous day
qpnd = dtp_ivol(sb) !m^3
sedpnd = dtp_ised(sb) !tons
!! Storage capacity under addon
qaddon = dtp_addon(sb,1)**3./(3.*dtp_lwratio(sb)*ch_s(2,sb)**2.) !m3 Note: V = d^3 / (3*R*S^2) |Modify by J. Osorio (3/19/2013)
!! iterate for subdaily flow/sediment routing
do ii=1,nstep
qout = 0.; qovmax = 0; depaddon = 0.
qin = hhvaroute(2,ihout,ii) + qpnd !m^3
sedin = hhvaroute(3,ihout,ii) + sedpnd !tons
if (qin>1e-6) then
spndconc = sedin / qin !initial sed conc, tons/m3
else
cycle
end if
!! Estimate water depth
qdepth = (3.*qin*dtp_lwratio(sb)*ch_s(2,sb)**2)**0.33333 !! Note: h = (3*V*R*S^2)^3 |Modify by J. Osorio (3/19/2013)
!! skip to next time step if no ponding occurs
if (qdepth<=0.0001) cycle
if (dtp_stagdis(sb)==0) then
!! Calculate weir outflow
do k = 1, dtp_numstage(sb)
qstage = 0.
!! calculate weir discharge
if (dtp_weirtype(sb,k)==2) then
!! Circular weir
dtp_depweir(sb,k) = dtp_diaweir(sb,k) + dtp_addon(sb,k)
if (qdepth>dtp_depweir(sb,k)) then
!! Fully submerged
qdepth = qdepth - dtp_depweir(sb,k)
watdepact = qdepth + dtp_diaweir(sb,k) / 2
warea = 3.14159 * dtp_diaweir(sb,k) ** 2 / 4.
!! orifice equation
qstage = dtp_cdis(sb,k) * 0.6 * warea *
& sqrt(19.6 * watdepact) !m3/s/unit
qstage = qstage * dtp_numweir(sb) * 60. * idt !m^3
else
!! Partially submerged
watdepact = max(qdepth - dtp_addon(sb,k),0.)
dtp_wrwid(sb,k) = dtp_diaweir(sb,k) * 0.667
!! weir/orifice discharge
qstage = dtp_cdis(sb,k) * 1.84 * dtp_wrwid(sb,k) *
& watdepact ** 1.5 !m3/s
qstage = qstage * dtp_numweir(sb) * 60. * idt !m^3
end if
else
!! Rectangular weir
watdepact = max(qdepth - dtp_addon(sb,k),0.)
!! Estimate weir/orifice discharge
qstage = dtp_cdis(sb,k) * 1.84 * dtp_wrwid(sb,k) *
& watdepact ** 1.5 !m3/s The Bureau of Reclamation, in their Water Measurement Manual, SI units
qstage = qstage * dtp_numweir(sb) * 60. * idt !m^3
end if
qout = qout + qstage
end do
!Limit total outflow amount less than available water above addon
if(qout>qin-qaddon) qout = max(qin - qaddon,0.)
!! Flow over the emergency weir
watdepact = qdepth - (dtp_depweir(sb,1) + dtp_addon(sb,1))
if (dtp_weirtype(sb,k)==1 .and. watdepact>0.) then
qstage = dtp_cdis(sb,k) * 1.84 *
& dtp_totwrwid(sb) * (watdepact ** 1.5) * 60. * idt !m3/s
qout = qout + qstage
end if
else
!! Use stage-discharge relationship if available
if (dtp_stagdis(sb)==1) then
select case(dtp_reltype(sb))
case(1) !! 1 is exponential function
qout = dtp_coef1(sb) * exp(dtp_expont(sb) * qdepth) +
& dtp_intcept(sb)
case(2) !! 2 is Linear function
qout = dtp_coef1(sb) * qdepth + dtp_intcept(sb)
case(3) !! 3 is logarthmic function
qout = dtp_coef1(sb) * log(qdepth) + dtp_intcept(sb)
case(4) !! 4 is power function
qout = dtp_coef1(sb) * (qdepth**3) + dtp_coef2(sb) *
& (qdepth**2) + dtp_coef3(sb) * qdepth + dtp_intcept(sb)
case(5)
qout = dtp_coef1(sb)*(qdepth**dtp_expont(sb))+
& dtp_intcept(sb)
end select
qout = qout * 60. * idt
end if !! end of stage-discharge calculation
end if
!! Check mass balance for flow
if (qout>qin) then !no detention occurs
qout = qin
qpnd = 0.
else !detention occurs
!! Estimating surface area of water
backup_length = qdepth / ch_s(2,sb)
seep_sa = backup_length/dtp_lwratio(sb)
& + (4. * dtp_lwratio(sb) * qdepth**2) / 3. !! Note: SSA = w + (4*l*d^2)/(3*w) |Modify by J. Osorio (3/20/2013)
evap_sa = 2. * backup_length**2 / 3. * dtp_lwratio(sb) !! Note: ESA = 2 * w * l / 3 |Modify by J. Osorio (3/20/2013)
!! converting surface area to hectares (ha)
seep_sa = seep_sa / 10000.0
evap_sa = evap_sa / 10000.0
!! Estimate rainfall, evapotranspiration, and seepage
pcp_vol = 10.0 * sub_subp_dt(sb,ii) * evap_sa !m^3
evap_vol = 10.0 * dtp_evrsv(sb) * pet_day * evap_sa !m^3
seep_vol = 10.0 * ch_k(2,sb) * seep_sa * idt / 60. !m^3
!! Check mass balance for water in the pond
qpnd = qin + pcp_vol - qout - evap_vol - seep_vol
if (qpnd<0) qpnd = 0.
end if
!! Mass balance for sediment
sedout = spndconc * qout !tons
sedpnd = max(0.,sedin - sedout) !tons
!! Store flow/sediment out of the pond at the subbasin outlet
hhvaroute(2,ihout,ii) = max(0.,qout)
hhvaroute(3,ihout,ii) = max(0.,sedout)
end do !! Outermost do loop ends here
! Store end-of-day values for next day
dtp_ivol(sb) = qpnd !m^3
dtp_ised(sb) = sedpnd !tons
end subroutine