c23456789012345678901234567890123456789012345678901234567890123456789012
c-----------------------------------------------------------------------
c This program is the solution of homework 3, ce632, spring 2000
c-----------------------------------------------------------------------
program ce632hw3sp00
parameter (mr=100,mt=1000)
parameter (lu5=5,lu6=6,lu7=7,lu8=8)
dimension q3(0:mr,0:mt),q3m(0:mr,0:mt)
dimension q4(0:mr,0:mt),q4m(0:mr,0:mt)
dimension q3no(0:mr,0:mt),q4no(0:mr,0:mt)
dimension q3nom(0:mr,0:mt),q4nom(0:mr,0:mt)
dimension time(-1:mt)
open(lu7,file='ce632hw3sp00.dat',status='unknown')
open(lu8,file='ce632hw3sp00.res',status='unknown')
data pi /3.1415926/
data qbase,qpeak,period,drop /50.,200.,96.,1./
data time(0) /0./
data xn,beta /0.02972,1.666667/
data tol/ 0.001/
data time(-1) /-6./
c-----------------------------------------------------------------------
c input data
c-----------------------------------------------------------------------
read(lu7,*) nr,nt,rl,tt
write(lu8,'(1x,a,i10)')'number of space steps= ',nr
write(lu8,'(1x,a,i10)')'number of time steps= ',nt
write(lu8,'(1x,a,f10.2)')'reach length (miles)= ',rl
write(lu8,'(1x,a,f10.2)')'total simulation time (hr)= ',tt
c-----------------------------------------------------------------------
c initial conditions
c-----------------------------------------------------------------------
dx= rl*5280./nr
dt= tt*3600./nt
pi2ovperiod= pi*2/period
dth= dt/3600.
slope= drop/5280.
alpha= (1.486/xn)*sqrt(slope)
rebeta= 1./beta
dtov3= dt/3.
do j= 0,nr
q3(j,0)= qbase
q4(j,0)= qbase
q3m(j,0)= qbase
q4m(j,0)= qbase
enddo
c-----------------------------------------------------------------------
c boundary conditions
c-----------------------------------------------------------------------
qa= 0.5*(qpeak + qbase)
qd= 0.5*(qpeak - qbase)
do n= 0,nt
q3(0,n)= qbase
q4(0,n)= qbase
q3m(0,n)= qbase
q4m(0,n)= qbase
enddo
t= -dth
n= -1
dowhile (t.lt.period)
t= t + dth
n= n+1
q3(0,n)= qa - qd*cos(pi2ovperiod*t)
q4(0,n)= q3(0,n)
q3m(0,n)= q3(0,n)
q4m(0,n)= q3(0,n)
q3no(0,n)= q3(0,n) - qbase
enddo
c-----------------------------------------------------------------------
c three-point variable-parameter muskingum-cunge routing
c-----------------------------------------------------------------------
do n= 1,nt
do j= 1,nr
cel0= beta*q3(j-1,n)*(q3(j-1,n)/alpha)**(-rebeta)
cel1= beta*q3(j-1,n-1)*(q3(j-1,n-1)/alpha)**(-rebeta)
cel2= beta*q3(j,n-1)*(q3(j,n-1)/alpha)**(-rebeta)
qref= 0.3333333*(q3(j-1,n) + q3(j-1,n-1) + q3(j,n-1))
cel= 0.3333333*(cel0 + cel1 + cel2)
cou= cel*dt/dx
rey= qref/(slope*cel*dx)
c3= 1./(1. + cou + rey)
c0= c3*(-1. + cou + rey)
c1= c3*( 1. + cou - rey)
c2= c3*( 1. - cou + rey)
q3(j,n)= c0*q3(j-1,n) + c1*q3(j-1,n-1) + c2*q3(j,n-1)
c-----------------------------------------------------------------------
c four-point variable-parameter muskingum-cunge routing
c-----------------------------------------------------------------------
qtrial= q3(j,n)
iflag= 0
dowhile (iflag.eq.0)
cel0= beta*q4(j-1,n)*(q4(j-1,n)/alpha)**(-rebeta)
cel1= beta*q4(j-1,n-1)*(q4(j-1,n-1)/alpha)**(-rebeta)
cel2= beta*q4(j,n-1)*(q4(j,n-1)/alpha)**(-rebeta)
cel3= beta*qtrial*(qtrial/alpha)**(-rebeta)
qref= 0.25*(q4(j-1,n) + q4(j-1,n-1) + q4(j,n-1) + qtrial)
cel= 0.25*(cel0 + cel1 + cel2 + cel3)
cou= cel*dt/dx
rey= qref/(slope*cel*dx)
c3= 1./(1. + cou + rey)
c0= c3*(-1. + cou + rey)
c1= c3*( 1. + cou - rey)
c2= c3*( 1. - cou + rey)
q4(j,n)= c0*q4(j-1,n) + c1*q4(j-1,n-1) + c2*q4(j,n-1)
if(abs(q4(j,n) - qtrial).lt.tol) then
iflag= 1
else
qtrial= q4(j,n)
endif
enddo
enddo
enddo
do n= 0,nt
do j= 0,nr
q3no(j,n)= q3(j,n) - qbase
q4no(j,n)= q4(j,n) - qbase
enddo
enddo
c-----------------------------------------------------------------------
c three-point modified variable-parameter muskingum-cunge routing
c-----------------------------------------------------------------------
do n= 1,nt
do j= 1,nr
qref= 0.3333333*(q3m(j-1,n) + q3m(j-1,n-1) + q3m(j,n-1))
cel= beta*qref*(qref/alpha)**(-rebeta)
cou= cel*dt/dx
rey= qref/(slope*cel*dx)
c3= 1./(1. + cou + rey)
c0= c3*(-1. + cou + rey)
c1= c3*( 1. + cou - rey)
c2= c3*( 1. - cou + rey)
q3m(j,n)= c0*q3m(j-1,n) + c1*q3m(j-1,n-1) + c2*q3m(j,n-1)
c-----------------------------------------------------------------------
c four-point modified variable-parameter muskingum-cunge routing
c-----------------------------------------------------------------------
qtrial= q3m(j,n)
iflag= 0
dowhile (iflag.eq.0)
qref= 0.25*(q4m(j-1,n) + q4m(j-1,n-1) + q4m(j,n-1) + qtrial)
cel= beta*qref*(qref/alpha)**(-rebeta)
cou= cel*dt/dx
rey= qref/(slope*cel*dx)
c3= 1./(1. + cou + rey)
c0= c3*(-1. + cou + rey)
c1= c3*( 1. + cou - rey)
c2= c3*( 1. - cou + rey)
q4m(j,n)= c0*q4m(j-1,n) + c1*q4m(j-1,n-1) + c2*q4m(j,n-1)
if(abs(q4m(j,n) - qtrial).lt.tol) then
iflag= 1
else
qtrial= q4m(j,n)
endif
enddo
enddo
enddo
do n= 0,nt
do j= 0,nr
q3nom(j,n)= q3m(j,n) - qbase
q4nom(j,n)= q4m(j,n) - qbase
enddo
enddo
c-----------------------------------------------------------------------
c writing output
c-----------------------------------------------------------------------
write(lu6,*)
1' step time inflow inflow outflow(3p) outflow(4p) outflow
2(3pm) outflow(4pm)'
write(lu8,*)
1' step time inflow inflow outflow(3p) outflow(4p) outflow
2(3pm) outflow(4pm)'
do n= 0,nt
time(n)= time(n-1) + dth
write(lu6,'(x,i5,f10.2,6f10.3)')
1n,time(n),q3(0,n),q3no(0,n),q3(nr,n),q4(nr,n),q3m(nr,n),q4m(nr,n)
write(lu8,'(x,i5,f10.2,6f10.3)')
1n,time(n),q3(0,n),q3no(0,n),q3(nr,n),q4(nr,n),q3m(nr,n),q4m(nr,n)
enddo
c-----------------------------------------------------------------------
c calculating mass conservation/with baseflow/conventional methods
c-----------------------------------------------------------------------
sumin23= 0.
sumout23= 0.
sumin43= 0.
sumout43= 0.
sumin24= 0.
sumout24= 0.
sumin44= 0.
sumout44= 0.
do n= 1,nt-1,2
sumin43= sumin43 + q3(0,n)
sumout43= sumout43 + q3(nr,n)
sumin44= sumin44 + q4(0,n)
sumout44= sumout44 + q4(nr,n)
enddo
do n= 2,nt-2,2
sumin23= sumin23 + q3(0,n)
sumout23= sumout23 + q3(nr,n)
sumin24= sumin24 + q4(0,n)
sumout24= sumout24 + q4(nr,n)
enddo
volin3= dtov3*(q3(0,0) + 4*sumin43 + 2*sumin23 + q3(0,nt))
volin4= dtov3*(q4(0,0) + 4*sumin44 + 2*sumin24 + q4(0,nt))
volout3= dtov3*(q3(nr,0) + 4*sumout43 + 2*sumout23 + q3(nr,nt))
volout4= dtov3*(q4(nr,0) + 4*sumout44 + 2*sumout24 + q4(nr,nt))
perc3= 100.*(volout3/volin3)
perc4= 100.*(volout4/volin4)
write(lu6,'(x,a,f10.3)')
1'percentage mass conservation/with baseflow/3p= ',perc3
write(lu8,'(x,a,f10.3)')
1'percentage mass conservation/with baseflow/3p= ',perc3
write(lu6,'(x,a,f10.3)')
1'percentage mass conservation/with baseflow/4p= ',perc4
write(lu8,'(x,a,f10.3)')
1'percentage mass conservation/with baseflow/4p= ',perc4
c-----------------------------------------------------------------------
c calculating mass conservation/with baseflow/modified methods
c-----------------------------------------------------------------------
sumin23= 0.
sumout23= 0.
sumin43= 0.
sumout43= 0.
sumin24= 0.
sumout24= 0.
sumin44= 0.
sumout44= 0.
do n= 1,nt-1,2
sumin43= sumin43 + q3m(0,n)
sumout43= sumout43 + q3m(nr,n)
sumin44= sumin44 + q4m(0,n)
sumout44= sumout44 + q4m(nr,n)
enddo
do n= 2,nt-2,2
sumin23= sumin23 + q3m(0,n)
sumout23= sumout23 + q3m(nr,n)
sumin24= sumin24 + q4m(0,n)
sumout24= sumout24 + q4m(nr,n)
enddo
volin3= dtov3*(q3m(0,0) + 4*sumin43 + 2*sumin23 + q3m(0,nt))
volin4= dtov3*(q4m(0,0) + 4*sumin44 + 2*sumin24 + q4m(0,nt))
volout3= dtov3*(q3m(nr,0) + 4*sumout43 + 2*sumout23 + q3m(nr,nt))
volout4= dtov3*(q4m(nr,0) + 4*sumout44 + 2*sumout24 + q4m(nr,nt))
perc3= 100.*(volout3/volin3)
perc4= 100.*(volout4/volin4)
write(lu6,'(x,a,f10.3)')
1'percentage mass conservation/with baseflow/3pm= ',perc3
write(lu8,'(x,a,f10.3)')
1'percentage mass conservation/with baseflow/3pm= ',perc3
write(lu6,'(x,a,f10.3)')
1'percentage mass conservation/with baseflow/4pm= ',perc4
write(lu8,'(x,a,f10.3)')
1'percentage mass conservation/with baseflow/4pm= ',perc4
c-----------------------------------------------------------------------
c calculating mass conservation/without baseflow/conventional method
c-----------------------------------------------------------------------
sumin23= 0.
sumout23= 0.
sumin43= 0.
sumout43= 0.
sumin24= 0.
sumout24= 0.
sumin44= 0.
sumout44= 0.
do n= 0,nt
q3no(0,n)= q3(0,n) - qbase
q4no(0,n)= q4(0,n) - qbase
q3no(nr,n)= q3(nr,n) - qbase
q4no(nr,n)= q4(nr,n) - qbase
enddo
do n= 1,nt-1,2
sumin43= sumin43 + q3no(0,n)
sumout43= sumout43 + q3no(nr,n)
sumin44= sumin44 + q4no(0,n)
sumout44= sumout44 + q4no(nr,n)
enddo
do n= 2,nt-2,2
sumin23= sumin23 + q3no(0,n)
sumout23= sumout23 + q3no(nr,n)
sumin24= sumin24 + q4no(0,n)
sumout24= sumout24 + q4no(nr,n)
enddo
volin3= dtov3*(q3no(0,0) + 4*sumin43 + 2*sumin23 + q3no(0,nt))
volin4= dtov3*(q4no(0,0) + 4*sumin44 + 2*sumin24 + q4no(0,nt))
volout3= dtov3*(q3no(nr,0) +4*sumout43 + 2*sumout23 + q3no(nr,nt))
volout4= dtov3*(q4no(nr,0) +4*sumout44 + 2*sumout24 + q4no(nr,nt))
perc3= 100.*(volout3/volin3)
perc4= 100.*(volout4/volin4)
write(lu6,'(x,a,f10.3)')
1'percentage mass conservation/without baseflow/3p= ',perc3
write(lu8,'(x,a,f10.3)')
1'percentage mass conservation/without baseflow/3p= ',perc3
write(lu6,'(x,a,f10.3)')
1'percentage mass conservation/without baseflow/4p= ',perc4
write(lu8,'(x,a,f10.3)')
1'percentage mass conservation/without baseflow/4p= ',perc4
c-----------------------------------------------------------------------
c calculating mass conservation/without baseflow/modified method
c-----------------------------------------------------------------------
sumin23= 0.
sumout23= 0.
sumin43= 0.
sumout43= 0.
sumin24= 0.
sumout24= 0.
sumin44= 0.
sumout44= 0.
do n= 0,nt
q3nom(0,n)= q3m(0,n) - qbase
q4nom(0,n)= q4m(0,n) - qbase
q3nom(nr,n)= q3m(nr,n) - qbase
q4nom(nr,n)= q4m(nr,n) - qbase
enddo
do n= 1,nt-1,2
sumin43= sumin43 + q3nom(0,n)
sumout43= sumout43 + q3nom(nr,n)
sumin44= sumin44 + q4nom(0,n)
sumout44= sumout44 + q4nom(nr,n)
enddo
do n= 2,nt-2,2
sumin23= sumin23 + q3nom(0,n)
sumout23= sumout23 + q3nom(nr,n)
sumin24= sumin24 + q4nom(0,n)
sumout24= sumout24 + q4nom(nr,n)
enddo
volin3= dtov3*(q3nom(0,0) + 4*sumin43 + 2*sumin23 + q3nom(0,nt))
volin4= dtov3*(q4nom(0,0) + 4*sumin44 + 2*sumin24 + q4nom(0,nt))
volout3= dtov3*(q3nom(nr,0) + 4*sumout43 +2*sumout23+q3nom(nr,nt))
volout4= dtov3*(q4nom(nr,0) + 4*sumout44 +2*sumout24+q4nom(nr,nt))
perc3= 100.*(volout3/volin3)
perc4= 100.*(volout4/volin4)
write(lu6,'(x,a,f10.3)')
1'percentage mass conservation/without baseflow/3pm= ',perc3
write(lu8,'(x,a,f10.3)')
1'percentage mass conservation/without baseflow/3pm= ',perc3
write(lu6,'(x,a,f10.3)')
1'percentage mass conservation/without baseflow/4pm= ',perc4
write(lu8,'(x,a,f10.3)')
1'percentage mass conservation/without baseflow/4pm= ',perc4
end
c-----------------------------------------------------------------------
c end
c-----------------------------------------------------------------------
|