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