C     Last change:  C     1 Mar 2001   12:27 pm
      program Landslide_Tadmor
c
c     Stiff source term solution using the second order 
c     central scheme of Tadmor and the modifications
c     suggested by Liotta, Romano and Russo, SIAM
c     J. Numer. Anal. 38(4), 1337-1356, 2000.
c     This is the Randau scheme.
c
c     Solving Sampson's parabolic canal problem
c============================================================
      implicit none
      integer md,nxmax,nxd, mn
      parameter (md=3, nxmax=1601, nxd=nxmax+2*md, mn=2)
      real*8 u(nxd,mn),x(nxd,2),z(nxd,2),dz(nxd,2)

      integer nx, nxx, io, i, isteps
      real*8 toll,cfl,theta,g,h0,dx,dtmax,Tout,Tfinal,tc,a,tau,B,t0,
     .       p,s,zeta_0,zeta_1,zeta,u0

      open(1,file='Stepflow5.out')
      open(2,file='Stepflow5.t')
      open(3,file='Stepflow5.h')
      open(4,file='Stepflow5.uh')
      open(5,file='Stepflow5.n')
      open(8,file='Stepflow5.c')
c
c Parameters
c      
      nx = 501
      nxx = nx+2*md
      if (nx.gt.nxmax) write(*,*)'nx too large'
      toll = 1.d-6
      cfl = 0.5d0
      theta = 1.0d0
      g = 9.81d0
      h0=10.0d0
      dx = 10000.d0/dble(nx-1)
      dtmax  = 0.5d0
      Tout   = 20.0d0
      Tfinal = 10000.0d0
      a = 3000.d0
      tau = 0.001d0
      B = 5.d0
      t0 = 0.d0

c
c Setup Terrain
c
      do io = 1,2
         do i = 1,nxx
            x(i,io) = dble(i-md-1)*dx+(io-1)*dx/2.0
         enddo
         call bed(x(1,io),z(1,io),dz(1,io),nxx)
      end do
      
c
c Initial Conditions
c 
      p = dsqrt(8.d0*g*h0/a/a)
      s = dsqrt(p*p-tau*tau)/2.d0
      u0 = B*exp(-tau*t0/2.d0)*sin(s*t0)
      zeta_0 = a*a*B*B*exp(-tau*t0)/(8.d0*g*g*h0)
     .         *(-s*tau*cos(s*t0)+(tau*tau/4.d0-s*s)*cos(2.d0*s*t0))
     .         -B*B*exp(-tau*t0)/4.d0/g
      zeta_1 = -exp(-tau*t0/2.d0)/g*(B*s*cos(s*t0)
     .         +tau*B/2.d0*sin(s*t0))
      do i = 1,nxx
        u(i,1)  = 20.d0 - h0*(1.d0 - (x(i,1) - 5000.d0)**2/a/a)
        zeta = zeta_0 + (x(i,1) - 5000.d0)*zeta_1 + 20.d0
        if(zeta.lt.u(i,1))then
          u(i,1) = 0.d0
        else
          u(i,1) = zeta - u(i,1)
        end if
        u(i,2) = 0.d0
      end do

c
c Setup for Time looping
c          
      isteps = 0
      tc = 0.0d0
      call dump_solution(u,z,nx,tc,isteps)
c
c Evolve
c
      call momentum_mass(nx,dx,cfl,g,theta,Tfinal,
     . u,dtmax,h0,toll,z,dz,Tout)
c
c Close down
c
      close(1)
      close(2)
      close(3)
      close(4)
      close(5)
      close(8)

      stop
      end

c===============================================================
c The evolver
c===============================================================
      subroutine momentum_mass(nx,dx,cfl,g,theta,tf,u,dtmax,h0,toll,
     .                         z,dz,Tout)
c
c INPUT   nx # of cells in x-direction
c         dx: step sizes in x-direction
c         cfl: CFL #   g: Acceleration due to gravity
c         tf: final time
c         theta=1: MM1 limiter; =2: MM2 limiter; >2: UNO limiter.
c         u: initial cell averages of conservative variables.
c            Supply entries of u((md+1):(nx+md),mn)
c OUTPUT  u: cell averages at final time 'tf'
c REMARK  1. Reset 'ndx' to adjust array dimensions.
c         2. Padded to each side of the computational domain are
c            'md' ghost cells, average values on which are
c            assigned by boundary conditions.
c
      implicit none
      integer md,nxmax,nxd, mn
      parameter (md=3, nxmax=1601, nxd=nxmax+2*md, mn=2)
      real*8 u(nxd,mn), ux(nxd,mn)
      real*8 f(nxd,mn), fx(nxd,mn)
      real*8 v(nxd), du(nxd,2), df(nxd,2), dz(nxd,2), u2(nxd,2), 
     .       u3(nxd,2), h(nxd), z(nxd,2), hz(nxd)
     
      integer nx, nxx, io, i, isteps
      real*8 toll,cfl,theta,g,h0,dx,dt,Tout,Tfinal,tc
      
      integer istop,iout,nt,ii,oi,i2,m
      real*8 dtmax, tf,a,b,xmin,t,xmic,Tleft,em_x,dtcdx2,dtcdx3
      real*8 den,xmt,pre, Tnext, em_old
      real*8 aterm1 ,term1, aterm2 ,term2, aterm3 ,term3  
c
      xmin(a,b) = 0.5d0*(dsign(1.d0,a)
     .          + dsign(1.d0,b))*dmin1(dabs(a),dabs(b))
      xmic(t,a,b) = xmin(t*xmin(a,b), 0.5d0*(a+b) )
c      
c----------------------------------------------------------
c 
      em_old = 1.0d0
      tc = 0.d0
      istop = 0
c      iplot = 100
      isteps = 0
      iout = 0
      nt =0
      Tleft = dmin1(Tout,tf)
      Tnext = dmin1(Tout,tf)
      dt = dmin1(dtmax,Tleft)
      nxx = nx + 2*md
c
      do while (.true.)
        nt = nt+1
        do ii = 1, 2
          io=ii-1
          oi = 1-io
          i2 = 3-ii
c
c Naive Boundary Conditions
c
          do i = 1, md
            u(i,1) = 0.d0
            u(nx+i,1) = 0.d0
            u(i,2) = 0.d0
            u(nx+i,2) = 0.d0
          end do
c
c Compute f and maximum wave speeds 'em_x'
c see (2.1) and (4.3)
c          
          call swflux(f,u,1,nxx,em_x,g,toll)
          
          call print_twoarray(f,nxx,'flux')
          call print_twoarray(u,nxx,'h and uh')          
c
c Compute numerical derivatives 'ux', 'fx'.
c see (3.1) and (4.1)
c           
          call derivative(u(1,1),ux(1,1),nxx,theta)
          call derivative(u(1,2),ux(1,2),nxx,theta)          
               
          call derivative(f(1,1),fx(1,1),nxx,theta)
          call derivative(f(1,2),fx(1,2),nxx,theta)    

          call print_twoarray(df,nxx,'df')
          call print_twoarray(ux,nxx,'ux')    
c
c Compute time step size according to the input CFL #
c
          if(ii.eq.1) then
            dt = dmin1(cfl*dx/em_x,dtmax,Tout)
c            write(*,*)em_x, dt
            if( ( tc + 2.d0*dt ) .ge. Tnext ) then
              dt = 0.5d0*(Tnext - tc )
              iout = 1
              Tnext = Tnext+Tout
            end if
            if (Tnext.gt.tf ) then
              istop = 1
            endif
          end if
c          write(*,*)'dt = ',dt
c          write(*,*)'Tnext = ',Tnext
          dtcdx2 = dt/dx/2.d0
          dtcdx3 = dt/dx/3.d0
c
c Compute the flux values of f at half time step,
c the conservative values at the half and third time steps.
c see (2.15) and (2.16).
          do i = 3, nxx-2
c Half time step
            u2(i,1) = u(i,1) - dtcdx2*fx(i,1)
            u2(i,2) = u(i,2) - dtcdx2*fx(i,2)
c Third Step
            u3(i,1) = u(i,1) - dtcdx3*fx(i,1)
            u3(i,2) = u(i,2) - dtcdx3*fx(i,2)
            if(u3(i,1).lt.toll)then
              u3(i,1) = 0.d0
            end if
          enddo

          call swflux(f,u2,3,nxx-2,em_x,g,toll)

          call print_twoarray(f,nxx,'flux half')
c
c Compute the values of 'u' at the next time level. see (2.16).
c Continuity equation
          m = 1
          aterm1 = 0.0d0
          aterm2 = 0.0d0
          aterm3 = 0.0d0 
                   
          do i = md + 1 - io, nx + md - io
            term1 =   (0.250d0   * (  u(i,m) +  u(i+1,m) )
     .             + 0.0625d0   * ( ux(i,m) - ux(i+1,m) )
     .             + dtcdx2     * (  f(i,m) -  f(i+1,m) ) )*2.d0
            v(i) = term1
            aterm1 = dmax1(aterm1,abs(term1))
          end do
          do i = md + 1, nx + md
            u(i,m) = v(i-io)
          end do

          call print_twoarray(u,nxx,'u after h update')
c
c Momentum equation
          m = mn

          do i = md + 1 - io, nx + md - io
            term2 =   (0.250d0   * (  u(i,m) +  u(i+1,m) )
     .             + 0.0625d0   * ( ux(i,m) - ux(i+1,m) )
     .             + dtcdx2     * (  f(i,m) -  f(i+1,m) ) )*2.d0
            v(i) = term2
            aterm2 = dmax1(aterm2,abs(term2))
            term3 = 0.5d0*dt*g*( u2(i,1)*dz(i,ii) + 
     .              u2(i+1,1)*dz(i+1,ii))
            aterm3 = dmax1(aterm3,abs(term3))
            v(i) = v(i) - term3
          end do

          do i = md + 1, nx + md
            u(i,m) = v(i-io)
c
c Friction term
            u(i,m) = u(i,m)*exp(-0.001d0*dt)
          end do
          call print_twoarray(u,nxx,'u after uh update')          
c
          tc = tc + dt
        end do

        if(iout.eq.1) then
           call dump_solution(u,z,nx,tc,isteps)
           iout = 0
        endif
        if(istop.eq.1) go to 1001
      end do
c
 1001 write(*,*)isteps
      write(5,*)isteps,nx
      return
 
      end

c===========================================================
c Setup Terrain with parabolic canal
c===========================================================
      subroutine bed(x,z,dz,n)

      implicit none
      integer md, nxmax,nxd,mn
      parameter (md=3, nxmax=1601, nxd=nxmax+2*md, mn=2)
      integer n
      real*8 z(n),dz(n),x(n)
      
      integer i
      real*8 h0, a

      h0 = 10.d0
      a = 3000.d0

      do i=1,n
         z(i)  = 20.d0 - h0*(1.d0 - (x(i) - 5000.d0)**2/a/a)
         dz(i) = 2.d0*h0*(x(i) - 5000.d0)/a/a
      enddo
      return
      end
c===========================================================
c Shallow Water Flux
c===========================================================
      subroutine swflux(f,q,n0,n1,em_x,g,toll)
      implicit real*8 (a-h)
      implicit real*8 (o-z)
      parameter (md=3, nxmax=1601, nxd=nxmax+2*md, mn=2)
      real*8 f(nxd,mn),q(nxd,mn)
      integer n0,n1
      
      
      em_x = 1.d-15
      do i = n0,n1
         h  = q(i,1)
         uh = q(i,2)
         if(h.le.toll)then
            u = 0.d0
            h = 0.0d0
            uh = 0.0d0
         else
            u = uh / h
         endif 
         if (abs(u).gt.1.0d3) then
            u = dsign(1.0d0,u)*1.0d3
            uh = u*h
         end if     
         cvel = dsqrt( g*h )
         em_x = dmax1( em_x, dabs(u) + cvel )
         f(i,1) = uh
         f(i,2) = uh*u + 0.5d0*g*h**2
c         q(i,1) = h
c         q(i,2) = uh
      enddo
      
      return
      end
c===========================================================
c Calculate limited derivatives
c===========================================================
      subroutine derivative(q,qx,n,theta)
      implicit real*8 (a-h)
      implicit real*8 (o-z)

      parameter (md=3, nxmax=1601, nxd=nxmax+2*md, mn=2)
      real*8 dq(nxd,2)         
      real*8 q(n),qx(n) 

      xmin(a,b) = 0.5d0*(dsign(1.d0,a)
     .          + dsign(1.d0,b))*dmin1(dabs(a),dabs(b))
      xmic(z,a,b) = xmin(z*xmin(a,b), 0.5d0*(a+b) )
      
      
      do i = 1,n-1
         dq(i,1) = q(i+1) - q(i)
      end do
      do i = 1,n-2
         dq(i,2) = dq(i+1,1) - dq(i,1)
      end do
      if( theta .lt. 2.5d0 ) then
         do i = 3,n-2
            qx(i) = xmic( theta, dq(i-1,1), dq(i,1) )
         end do
      else
         do i = 3,n-2
           qx(i) = xmin(dq(i-1,1) 
     .                  + 0.5d0*xmin(dq(i-2,2),dq(i-1,2)),
     .                    dq(i,1) - 0.5d0*xmin(dq(i-1,2),dq(i,2)))
        end do
      end if
      return
      end      
      
      
         
c===========================================================
c Output Routine
c===========================================================
      subroutine dump_solution(u,z,nx,tc,isteps)
      implicit real*8 (a-h)
      implicit real*8 (o-z)
      parameter (md=3, nxmax=1601, nxd=nxmax+2*md, mn=2)
      real*8 u(nxd,mn),z(nxd,2)

      write(2,*)tc
      io = 1
      isteps = isteps + 1
      write(*,*)isteps, tc
      do ii = 1,2
        do i = md + 1, nx + md
           write(3,200)u(i,1)
           write(4,200)u(i,2)
           write(1,200)z(i,io) + u(i,1)
  200 format(e20.8)
        end do
      end do

      iflag = 0
      do i = md + 1,nx + md
        if(u(i,1).gt.1.d-04.and.iflag.eq.0)then
c          write(8,200)u(i,1)
          write(8,200)i*20.d0
          iflag = 1
        end if
      end do
      return
      end
c===========================================================
c Output Routine
c===========================================================
      subroutine print_twoarray(u,nxx,msg)
      implicit real*8 (a-h)
      implicit real*8 (o-z)
      parameter (md=3, nxmax=1601, nxd=nxmax+2*md, mn=2)
      real*8 u(nxd,mn)
      character*(*) msg
      return
      write(*,*)msg

      do i = 1,nxx
         write(*,*)u(i,1),u(i,2)
      end do     
      return
      end