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 Landslide 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============================================================
      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,dx,dtmax,Tout,Tfinal,tc,a,tau,B,t0,
     .       p,s,zeta_0,zeta_1,zeta,u0,version,hl,hr

      open(11,file='data.w')
      open(12,file='data.t')
      open(13,file='data.h')
      open(14,file='data.uh')
      open(15,file='data.n')
      open(16,file='data.x')

      open(1,file='MacDonald.bed')
c
c Parameters
c   
      version = 1.0   
      nx = 101
      nxx = nx+2*md
      if (nx.gt.nxmax) write(*,*)'nx too large'
      toll = 1.d-6
      cfl = 0.2d0
      theta = 1.0d0
      g = 9.81d0
      hl= 5.0d0
      hr = 0.d0
      dx = 10000.d0/dble(nx-1)
      dtmax  = 0.5d0
      Tout   = 1.0d0
      Tfinal = 1500.0d0

      t0 = 0.d0

c
c Setup Terrain
c
c ---- Non-uniform terrain from MacDonald
      io = 1
      do i = 1,nxx
        read(1,*)x(i,io),z(i,io)
        z(i,io) = z(i,io)*5.d0
      end do
      do i = 2,nxx-1
        dz(i,io) = (z(i+1,io) - z(i-1,io))/2.d0/dx
      end do
      dz(1,io) = dz(2,io)
      dz(nxx,io) = dz(nxx-1,io)

      io = 2
      do i = 1,nxx
         x(i,io) = x(i,io-1)
         z(i,io) = z(i,io-1)
         dz(i,io) = dz(i,io-1)
      end do
      
c
c Initial Conditions
c 
c ---- Rectangular pulse initial condition
      do i = md+1, nx+md
        if (x(i,io).lt.100.d0.or.x(i,io).gt.200.d0) then
          u(i,1) = hr
        else
          u(i,1) = hl
        endif
        u(i,2) = 0.0d0
      end do  

c
c Setup for Time looping
c          

      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,toll,z,dz,Tout,isteps)

      write(15,*)isteps,nxx,version
      do i=1,nxx
         write(16,*) x(i,1)
      enddo


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,toll,
     .                         z,dz,Tout,isteps)
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,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) 
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)    

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

          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)

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

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
c            u(i,m) = u(i,m)*exp(-0.001d0*dt)
          end do
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

      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(12,*)tc
      io = 1
      i = 2
      isteps = isteps + 1
      write(*,*)isteps, tc

      do i = 1, nx + 2*md
         write(13,200)u(i,1)
         write(14,200)u(i,2)
         write(11,200)z(i,io) + u(i,1)
  200 format(e20.8)
      end do

      return
      end