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----------------------------------------------------------------------- |