Index: inundation/ga/storm_surge/Hobart/Landslide_Tadmor.for
===================================================================
--- inundation/ga/storm_surge/Hobart/Landslide_Tadmor.for	(revision 1451)
+++ inundation/ga/storm_surge/Hobart/Landslide_Tadmor.for	(revision 1451)
@@ -0,0 +1,434 @@
+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
