Index: inundation/ga/storm_surge/Hobart/Landslide_1D.for
===================================================================
--- inundation/ga/storm_surge/Hobart/Landslide_1D.for	(revision 1444)
+++ inundation/ga/storm_surge/Hobart/Landslide_1D.for	(revision 1444)
@@ -0,0 +1,830 @@
+      program sfbasic
+c============================================================
+c Combine evolution of w = h+z and h for shallow flows
+c============================================================
+
+      implicit none
+
+      integer md,nxmax,nxd, mn
+
+      parameter (md=3, nxmax=5000, nxd=nxmax+2*md, mn=2)
+      real*8 h(nxd), uh(nxd), w(nxd)
+      real*8 z(nxd),z2(nxd)
+      real*8 x(nxd),x2(nxd)
+      
+      integer terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+      common /setup/ terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+
+      integer nlow,nhigh,nx,nxx,dn,version,isteps      
+      common /ivariables/ nlow,nhigh,nx,nxx,version,isteps
+      
+      real*8 toll,cfl,theta
+      real*8 g,hf,dx,dt,dtmax
+      real*8 Tout,tf,tc,uhf
+      real*8 K,limith,length
+      real*8 manning,hl,hr,zb      
+      common /rvariables/ toll,cfl,theta,
+     .                    g,hf,dx,dt,dtmax,
+     .                    Tout,tf,tc,uhf,
+     .                    K,limith,length,
+     .                    manning,hl,hr
+          
+      integer io, i
+
+      call util_open_files
+
+      open(1,file='MacDonald.bed')
+c
+c Parameters
+c      
+      nx = 101
+      dtmax  = 1.0d0
+      Tout   = 1.0d0
+      tf = 200.0d0
+      length = 1000.0d0
+
+      nlow = md+2
+      nhigh = md+nx
+      nxx = nx+2*md
+      
+      limith = 1
+      dx = 10.d0
+      K = 0.0d0            
+
+      if (nx.gt.nxmax) write(*,*)'nx too large'
+      toll = 1.0d-10
+      cfl = 0.01d0
+      theta = 1.2d0
+      version = 10
+      g = 9.81d0
+      hl = 5.d0
+      hr = 0.d0
+      tc = 0.d0
+      
+c
+c ---- Setup Terrain, Boundary Values and Initial Values 
+c
+
+c ---- Non-uniform terrain from MacDonald
+      do i = 1,nxx
+        read(1,*)x2(i),z2(i)
+        z2(i) = z2(i)*5.d0
+        x(i) = x2(i) - dx/2.d0
+      end do
+      do i = 2,nxx
+        z(i) = (z2(i-1)+z2(i))/2.d0
+      end do
+      z(1) = (3.d0*z2(1) - z2(2))/2.d0
+      z(nxx) = (3.d0*z2(nxx-1) - z2(nxx-2))/2.d0
+
+      call boundary_apply(x,h,uh,z)
+
+c ---- Rectangular pulse initial condition
+      do i = md+1, nx+md
+        if (x(i).lt.100.d0.or.x(i).gt.200.d0) then
+          h(i) = hr
+        else
+          h(i) = hl
+        endif
+        uh(i) = 0.0d0
+      end do  
+c
+c ---- Setup for Time looping
+c          
+      isteps = 0
+      tc = 0.0d0
+      call util_dump_solution(h,uh,z)
+c
+c ---- Evolve
+c
+      call evolver_euler(x,h,uh,z,z2)
+c
+c ---- Close down
+c
+      write(15,*)isteps,nxx,version
+      do i=1,nxx
+         write(16,*) x(i)
+      enddo
+      
+      do i=5,nxx-3
+         write(25,*) x(i)
+      enddo
+
+      call util_close_files
+
+      stop
+      end
+c===========================================================
+c Shallow Water Flux
+c===========================================================
+      subroutine flux_huh(f1,f2,h0,uh0,em_x,l1,l2,g,toll)
+      implicit none
+
+      real*8 f1,f2,h0,uh0,z,l1,l2
+      real*8 em_x,g,toll
+      
+      real*8 u,cvel,h,uh,w
+      integer i
+
+c-----------------------------------
+      
+      if(em_x.lt.1.0d-15) then
+         em_x = 1.d-15
+      endif
+  
+      uh = uh0
+      h = h0
+      
+      if(h.le.toll)then
+         u = 0.d0
+         h = 0.0d0
+         uh = 0.0d0
+      else
+         u = uh/h
+      endif
+    
+      cvel = dsqrt( g*h )
+      l1 = u - cvel
+      l2 = u + cvel
+      em_x = dmax1( em_x, dabs(u) + cvel )
+      f1 = uh
+      f2 = uh*u + 0.5d0*g*h**2
+      
+      return
+      end
+c===========================================================
+c Shallow Water Flux
+c===========================================================
+      subroutine flux_huh_old(f1,f2,h0,uh0,z,em_x,l1,l2,g,toll)
+      implicit none
+
+      real*8 f1,f2,h0,uh0,z,l1,l2
+      real*8 em_x,g,toll
+      
+      real*8 u,cvel,h,uh,w
+      integer i
+
+c-----------------------------------
+      
+      if(em_x.lt.1.0d-15) then
+         em_x = 1.d-15
+      endif
+  
+      w = h0+z
+      uh = uh0
+      h = h0
+      
+      if(h.le.toll)then
+         u = 0.d0
+         h = 0.0d0
+         uh = 0.0d0
+      else
+         u = uh / h
+      endif 
+      if (abs(u).gt.1.0d1) then
+         u = dsign(1.0d0,u)*1.0d3
+         uh = u*h
+      end if     
+      cvel = dsqrt( g*h )
+      l1 = u - cvel
+      l2 = u + cvel
+      em_x = dmax1( em_x, dabs(u) + cvel )
+      f1 = uh
+      f2 = uh*u + 0.5d0*g*h**2
+      
+      return
+      end
+c===========================================================
+c Shallow Water Flux using h+z and uh variables
+c===========================================================
+      subroutine flux_wuh(f1,f2,w0,uh0,z,em_x,l1,l2,g,toll)
+      implicit none
+
+      real*8 f1,f2,w0,uh0,z
+      real*8 g,toll
+      
+      real*8 h,u,cvel,uh,w,l1,l2,em_x
+      integer i
+
+c-----------------------------------
+
+      if(em_x.lt.1.0d-15) then
+         em_x = 1.d-15
+      endif
+  
+      w = w0
+      uh = uh0
+      
+      h = w - z
+      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 )
+      l1 = u - cvel
+      l2 = u + cvel
+      em_x = dmax1( em_x, dabs(u) + cvel )
+      f1 = uh
+      f2 = uh*u + 0.5d0*g*h**2
+  
+      return
+      end
+c=============================================================== 
+c array limiter
+c===============================================================
+      subroutine array_limit_h(h,uh,z2,z,uh_l,uh_r,h_l,h_r,h_2)
+
+      implicit none
+      integer md,nxmax,nxd, mn
+      parameter (md=3, nxmax=5000, nxd=nxmax+2*md, mn=2)
+      real*8 h(nxd), uh(nxd), w(nxd)
+      real*8 z(nxd),z2(nxd)
+      real*8 x(nxd),x2(nxd)
+  
+      integer terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+      common /setup/ terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+
+      integer nlow,nhigh,nx,nxx,dn,version,isteps      
+      common /ivariables/ nlow,nhigh,nx,nxx,version,isteps
+      
+      real*8 toll,cfl,theta
+      real*8 g,hf,dx,dt,dtmax
+      real*8 Tout,tf,tc,uhf
+      real*8 K,limith,length
+      real*8 manning,hl,hr     
+      common /rvariables/ toll,cfl,theta,
+     .                    g,hf,dx,dt,dtmax,
+     .                    Tout,tf,tc,uhf,
+     .                    K,limith,length,
+     .                    manning,hl,hr
+          
+      real*8 h_2(nxd)
+      real*8 h_l(nxd),h_r(nxd)
+      real*8 h1_l(nxd), h1_r(nxd)
+      real*8 uh_l(nxd), uh_r(nxd)
+      
+      integer sslope(nxd)
+      real*8 xmin,xmic,xminmod,a,b,c,t
+
+      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) )
+      xminmod(a,b,c) = xmin(a,xmin(b,c))
+
+      
+      integer i    
+      real*8 dh, hmax, hmin , tt, cc
+      real*8 w0,w1,w2,d0h,d0w,d1h,d1w,d2h,d2w,ddh,ddw
+      real*8 h0,h1,h2,hz
+      real*8 ss0,ss1,ss
+      real*8 dz0,dz1,dz2
+      real*8 hx1,hx2,h_x
+      real*8 froud,d1,d2
+      real*8 gh3,uh2,duh, w_x,dz
+      
+c------------------------------------ 
+      g = 9.81d0     
+      cc = 0.1d0
+      ss0 = 0.01d0
+      ss1 = 0.01d0
+
+      do i=md,nx+md+2
+
+         h_2(i) = h(i)
+         h0 = h(i-1)
+         h1 = h(i)
+         h2 = h(i+1)
+         
+         dz0 = abs(z2(i-1)-z2(i-2))
+         dz1 = abs(z2(i)  -z2(i-1))
+         dz2 = abs(z2(i+1)-z2(i)  )
+
+         w0 = h0 + z(i-1)
+         w1 = h1 + z(i)
+         w2 = h2 + z(i+1)
+         
+         d1w = w1 - w0
+         d2w = w2 - w1
+         dz  = (z2(i) - z2(i-1))
+                
+c ------- W = H+Z           
+         h_x = xmin( d1w, d2w) - dz
+         
+c ------- Return the left and right values of h in an interval
+         h_l(i)  = h(i)   - 0.5d0*h_x
+         h_r(i)  = h(i)   + 0.5d0*h_x 
+         
+c ------- Test for h_l or h_r negative
+         hmin = dmin1(h0,h1,h2)
+         if ( h_l(i) .le. hmin ) then
+            h_l(i) = hmin
+            h_r(i) = h(i) + (h(i) - h_l(i))
+         elseif ( h_r(i) .le. hmin ) then
+            h_r(i) = hmin
+            h_l(i) = h(i) + (h(i) - h_r(i))
+         endif
+      enddo
+c
+c ---- Enforce a minimum height
+c
+      do i=md,nx+md+2
+         hmax = dmax1(h_l(i),h_r(i))
+         if (hmax .lt. 1.0d-3 ) then
+           h_r(i)  = 0.0d0
+           uh_r(i) = 0.0d0
+           h_l(i)  = 0.0d0
+           uh_l(i) = 0.0d0
+           h(i)    = 0.0d0
+           uh(i)   = 0.0d0
+         endif
+      enddo
+c
+c ---- Try for a flat subinterval 
+c
+      do i=md,nx+md+2
+         hmax = dmax1(h_l(i),h_r(i))
+         dz  = (z2(i) - z2(i-1))
+
+         if (h_l(i) .lt. 0.01d0*h_r(i) ) then
+           h_r(i)  = sqrt(abs(2.0d0*dz*h(i)))
+         endif
+
+         if (h_r(i) .lt. 0.01d0*h_l(i) ) then         
+           h_l(i)  = sqrt(abs(2.0d0*dz*h(i)))       
+         endif         
+      enddo
+            
+      return
+      end
+c=============================================================== 
+c array limiter
+c===============================================================
+      subroutine array_limit(u,u_l,u_r)
+
+      implicit none
+      integer md,nxmax,nxd, mn
+      parameter (md=3, nxmax=5000, nxd=nxmax+2*md, mn=2)
+      real*8 h(nxd), uh(nxd), w(nxd)
+      real*8 z(nxd),z2(nxd)
+      real*8 x(nxd),x2(nxd)
+   
+      integer terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+      common /setup/ terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+
+      integer nlow,nhigh,nx,nxx,dn,version,isteps      
+      common /ivariables/ nlow,nhigh,nx,nxx,version,isteps
+      
+      real*8 toll,cfl,theta
+      real*8 g,hf,dx,dt,dtmax
+      real*8 Tout,tf,tc,uhf
+      real*8 K,limith,length
+      real*8 manning,hl,hr     
+      common /rvariables/ toll,cfl,theta,
+     .                    g,hf,dx,dt,dtmax,
+     .                    Tout,tf,tc,uhf,
+     .                    K,limith,length,
+     .                    manning,hl,hr
+        
+      real*8 u(nxd), u_l(nxd),u_r(nxd)
+      real*8 u_x,d1,d2,d0,ux1
+      real*8 u1_l(nxd),u1_r(nxd)
+      integer i
+      
+      real*8 a,b,t,xmin,xmic
+         
+      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------------------------------------
+      
+      do i=md,nx+md+2
+          d1 = u(i)   - u(i-1)
+          d2 = u(i+1) - u(i)
+          u_x = xmic( theta, d1, d2 ) 
+          u_l(i) = u(i) - 0.5d0*u_x
+          u_r(i) = u(i) + 0.5d0*u_x          
+      enddo
+             
+      return
+      end 
+c===========================================================
+c Output Routine
+c===========================================================
+      subroutine util_dump_solution(h,uh,z)
+
+      implicit none
+      integer md,nxmax,nxd, mn
+      parameter (md=3, nxmax=5000, nxd=nxmax+2*md, mn=2)
+      real*8 h(nxd), uh(nxd), w(nxd)
+      real*8 z(nxd),z2(nxd)
+      real*8 x(nxd),x2(nxd)
+      
+      integer terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+      common /setup/ terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+
+      integer nlow,nhigh,nx,nxx,dn,version,isteps      
+      common /ivariables/ nlow,nhigh,nx,nxx,version,isteps
+      
+      real*8 toll,cfl,theta
+      real*8 g,hf,dx,dt,dtmax
+      real*8 Tout,tf,tc,uhf
+      real*8 K,limith,length
+      real*8 manning,hl,hr     
+      common /rvariables/ toll,cfl,theta,
+     .                    g,hf,dx,dt,dtmax,
+     .                    Tout,tf,tc,uhf,
+     .                    K,limith,length,
+     .                    manning,hl,hr
+     
+      integer ii,i
+
+      write(12,*)tc
+      
+      isteps = isteps + 1
+      write(*,*)isteps, tc
+
+      do i = 1, nx+2*md
+c------ Dump h
+           write(13,200)sngl(h(i))
+c------ Dump uh
+           write(14,200)sngl(uh(i))
+c------ Dump w
+           write(11,200)sngl(h(i)+z(i))
+  200 format(e20.8)
+      end do
+
+      return
+      end
+c=============================================================== 
+c The evolver using euler's method
+c===============================================================
+      subroutine evolver_euler(x,h,uh,z,z2)
+
+      implicit none
+      integer md,nxmax,nxd, mn
+      parameter (md=3, nxmax=5000, nxd=nxmax+2*md, mn=2)
+      real*8 h(nxd), uh(nxd), w(nxd)
+      real*8 z(nxd),z2(nxd)
+      real*8 x(nxd),x2(nxd)
+      
+      integer terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+c      common /setup/ terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+
+      integer nlow,nhigh,nx,nxx,dn,version,isteps      
+      common /ivariables/ nlow,nhigh,nx,nxx,version,isteps
+      
+      real*8 toll,cfl,theta
+      real*8 g,hf,dx,dt,dtmax
+      real*8 Tout,tf,tc,uhf
+      real*8 K,limith,length
+      real*8 manning,hl,hr     
+      common /rvariables/ toll,cfl,theta,
+     .                    g,hf,dx,dt,dtmax,
+     .                    Tout,tf,tc,uhf,
+     .                    K,limith,length,
+     .                    manning,hl,hr
+     
+      real*8 h_new(nxd),uh_new(nxd)
+      real*8 h_h(nxd),uh_h(nxd)
+      real*8 h_h_x(nxd),uh_h_x(nxd)
+      real*8 h_x(nxd), uh_x(nxd)
+      real*8 f_h(nxd), f_uh(nxd)
+      real*8 f_h_l(nxd), f_uh_l(nxd)
+      real*8 f_h_r(nxd), f_uh_r(nxd)
+      real*8 z_l(nxd),z_r(nxd)
+      real*8 h_2(nxd), uh_2(nxd)
+
+      real*8 h_l(nxd), h_r(nxd), uh_l(nxd), uh_r(nxd)     
+      real*8 z2l, z2r, zl, zr
+       
+      integer io, i, jj, ipos,ip
+      real*8 factor, Ki
+
+      integer istop,iout,nt,ii,oi,i2,m
+      real*8 a,b,t,Tleft,em_x,dtcdx,dtcdx2,dtcdx3
+      real*8 den, xmt,pre, Tnext, em_old, hh
+      real*8 aterm1 ,term1, aterm2 ,term2, aterm3 ,term3
+
+      integer kk
+      real*8 umax, umax2, vmin, q1,q2,q3,h1,h2
+      real*8 Kcal,k1,k2,za,zb,pi
+
+      real*8 fh_l, fuh_l, fh_r, fuh_r
+      real*8 l1_l, l2_l, l1_r, l2_r
+      real*8 a_max, a_min, aa, axa
+ 
+c---------------------------------------------------------- 
+c
+      pi = 2.d0*dtan(1.d0)
+      factor = 1.0d0
+      em_old = 1.0d0
+      tc = 0.d0
+      istop = 0
+      iout = 0
+      nt =0
+      Tleft = dmin1(Tout,tf)
+      Tnext = dmin1(Tout,tf)
+      dt = dmin1(dtmax,Tleft)
+      
+      do i=1,nxx
+         f_h(i)  = 0.0d0
+         f_uh(i) = 0.0d0
+         h_l(i)  = h(i)
+         h_r(i)  = h(i)
+         uh_l(i) = uh(i)
+         uh_r(i) = uh(i)
+      enddo
+
+      call boundary_apply(x,h,uh,z)
+      call array_limit  (uh,uh_l,uh_r)         
+      call array_limit_h(h,uh,z2,z,uh_l,uh_r,h_l,h_r,h_2) 
+            
+      do i=1,nxx
+         if ( h(i) .lt. 0.0d0 ) then
+           write(*,*)'i=',i,' h(i) = ',h(i)
+         endif
+      enddo
+
+c------------------------------------------------------
+c
+      do while (.true.)
+c          write(*,*)nt
+          nt = nt+1
+          ii = 1
+
+          call boundary_apply(x,h,uh,z)
+ 
+          call array_limit  (uh,uh_l,uh_r)         
+          call array_limit_h(h,uh,z2,z,uh_l,uh_r,h_l,h_r,h_2) 
+     
+          do i=1,nxx
+            if ( h(i) .lt. 0.0d0 ) then
+              write(*,*)'i=',i,' h(i) = ',h(i),' nt = ',nt
+            endif
+          enddo   
+
+c-------------------------
+c Flux Calculation
+c-------------------------
+          call numerical_flux_central(h,uh,z2,z,uh_l,uh_r,h_l,h_r,h_2,
+     .                       f_h,f_uh, em_x)
+
+c-----------------------------------------------------
+c Compute time step size according to the input CFL # 
+c-----------------------------------------------------
+ 1002     dt = dmin1(factor*cfl*dx/em_x,dtmax,Tout)
+c          write(*,*)'dt = ',dt,' factor = ',factor
+          if( ( tc + dt ) .ge. Tnext ) then
+              dt = (Tnext - tc )
+              iout = 1
+              Tnext = Tnext+Tout
+          end if
+          if (Tnext.gt.tf ) then
+              istop = 1
+          endif
+          
+c-------------------------------------------------
+c Compute the values of 'u' at the next time level
+c-------------------------------------------------
+          ipos = 0
+
+          do i = md+2, md+nx
+             
+             h_new(i)  =  h(i)  - dt/dx*( f_h(i)  - f_h(i-1) )
+             uh_new(i) =  uh(i) - dt/dx*( f_uh(i) - f_uh(i-1)) 
+     .              - dt/dx*g*(z2(i)-z2(i-1))*(h(i))
+
+             if (h_new(i).lt.-1.0d-5) then
+                write(*,*)'h_new(',i,') = ',h_new(i)
+                ipos=1
+             endif
+             
+             if (h_new(i).le.0.0d0) then
+                h_new(i) = 0.0d0
+c - my statement
+                uh_new(i) = 0.d0
+             endif
+             
+          end do
+
+          if(ipos.eq.1) then
+             factor = factor/2.0d0
+             write(*,*)'factor = ',factor
+             goto 1002
+          endif
+
+          do i = md+2,nx+md
+             h(i) = h_new(i)
+             uh(i) = uh_new(i)
+          enddo
+          
+          factor = 1.0d0
+          tc = tc + dt
+
+          if(iout.eq.1) then
+             call util_dump_solution(h,uh,z)
+             
+             call boundary_apply(x,h,uh,z)
+             call array_limit  (uh,uh_l,uh_r)         
+             call array_limit_h(h,uh,z2,z,uh_l,uh_r,h_l,h_r,h_2) 
+             iout = 0
+          endif          
+          if(istop.eq.1) go to 1001
+c 
+c End timestep 
+c
+      end do
+c
+ 1001 write(*,*)isteps
+
+      return
+
+      end
+c===========================================================
+c Boundary Conditions
+c
+c Naive Boundary Conditions
+c Inflow, or reflective
+c
+c===========================================================
+      subroutine boundary_apply(x,h,uh,z)
+
+      implicit none
+      integer md,nxmax,nxd, mn
+      parameter (md=3, nxmax=5000, nxd=nxmax+2*md, mn=2)
+      real*8 h(nxd), uh(nxd), w(nxd)
+      real*8 z(nxd),z2(nxd)
+      real*8 x(nxd),x2(nxd)
+      
+      integer terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+      common /setup/ terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+
+      integer nlow,nhigh,nx,nxx,dn,version,isteps      
+      common /ivariables/ nlow,nhigh,nx,nxx,version,isteps
+      
+      real*8 toll,cfl,theta
+      real*8 g,hf,dx,dt,dtmax
+      real*8 Tout,tf,tc,uhf
+      real*8 K,limith,length
+      real*8 manning,hl,hr     
+      common /rvariables/ toll,cfl,theta,
+     .                    g,hf,dx,dt,dtmax,
+     .                    Tout,tf,tc,uhf,
+     .                    K,limith,length,
+     .                    manning,hl,hr
+           
+      real*8 du,dh,uf,h0,u0,uh0
+      integer ii,i
+      
+c--------------------------------------------------
+
+      do i = 1, md      
+c==================
+c Left Boundary
+c==================
+
+c --- Zero
+        h(i)  =  0.d0
+        uh(i) =  0.d0
+      end do
+              
+      do i = 1, md         
+c==================
+c Right Boundary
+c==================
+
+c --- Zero
+        h(nx+md+i)  =  0.d0
+        uh(nx+md+i) =  0.d0
+c        h(nx+md+1) = h(nx-1)
+c        uh(nx+md+1) = uh(nx-1)
+      end do
+
+      return
+      end
+c===========================================================
+c Numerical Central Flux
+c===========================================================
+      subroutine numerical_flux_central(h,uh,z2,z,uh_l,uh_r,h_l,h_r,
+     .          h_2,f_h,f_uh, em_x)
+
+      implicit none
+      integer md,nxmax,nxd, mn
+      parameter (md=3, nxmax=5000, nxd=nxmax+2*md, mn=2)
+      real*8 h(nxd), uh(nxd), w(nxd)
+      real*8 z(nxd),z2(nxd)
+      real*8 x(nxd),x2(nxd)
+
+c      integer terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+c      common /setup/ terrain,bc_l,bc_r,initial,s_h,s_uh,s_semi_uh
+
+      integer nlow,nhigh,nx,nxx,dn,version,isteps      
+      common /ivariables/ nlow,nhigh,nx,nxx,version,isteps
+      
+      real*8 toll,cfl,theta
+      real*8 g,hf,dx,dt,dtmax
+      real*8 Tout,tf,tc,uhf
+      real*8 K,limith,length
+      real*8 manning,hl,hr     
+      common /rvariables/ toll,cfl,theta,
+     .                    g,hf,dx,dt,dtmax,
+     .                    Tout,tf,tc,uhf,
+     .                    K,limith,length,
+     .                    manning,hl,hr
+     
+      real*8 f_h(nxd), f_uh(nxd)
+      real*8 h_2(nxd), uh_2(nxd)
+      real*8 h_l(nxd), h_r(nxd), uh_l(nxd), uh_r(nxd)
+           
+
+      real*8 fh_l, fuh_l, fh_r, fuh_r
+      real*8 l1_l, l2_l, l1_r, l2_r
+      real*8 a_max, a_min, aa, axa
+      
+      integer istop,iout,nt,ii,oi,i2,m,i
+      real*8 a,b,t,Tleft,em_x,dtcdx,dtcdx2,dtcdx3
+c      real*8 den, xmt,pre, Tnext, em_old, hh
+c      real*8 aterm1 ,term1, aterm2 ,term2, aterm3 ,term3
+            
+c-----------------------------------
+   
+      em_x = 1.0d-15
+
+      do i = md,nxx-md+1
+
+c
+c ---- Calculate flux
+c
+
+         call flux_huh(fh_l,fuh_l,h_r(i),uh_r(i),
+     .                 em_x,l1_l,l2_l,g,toll)
+         call flux_huh(fh_r,fuh_r,h_l(i+1),uh_l(i+1),
+     .                 em_x,l1_r,l2_r,g,toll)
+c
+c --- Calculate max wave speed
+c         
+         a_max = dmax1(l2_l,l2_r,0.0d0)
+         a_min = dmin1(l1_l,l1_r,0.0d0)         
+c 
+c --- Compute the numerical flux 
+c 
+         aa = a_max - a_min
+         axa = a_max*a_min
+         if ( aa .gt.  dx*dt/1000.0d0) then
+            h_2(i)  =(a_max*h_l(i+1) - a_min*h_r(i) -
+     .                fh_r + fh_l )/aa
+            f_h(i)  =(a_max*fh_l  - a_min*fh_r  + 
+     .                axa*(h_l(i+1)-h_r(i)))/aa
+            f_uh(i) =(a_max*fuh_l - a_min*fuh_r + 
+     .                axa*(uh_l(i+1)-uh_r(i)))/aa
+         else
+            h_2(i)  = 0.0d0
+            f_h(i)  = 0.0d0
+            f_uh(i) = 0.0d0
+         endif 
+            
+      enddo
+      
+      return
+      end
+c===========================================================
+c Open Files
+c===========================================================
+      subroutine util_open_files
+      
+      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(25,file='data.xx') 
+      
+      return
+      end
+c===========================================================
+c Close Files
+c===========================================================
+      subroutine util_close_files
+      
+      close(11)
+      close(12)
+      close(13)
+      close(14)
+      close(15)
+      close(16)
+
+      close(25)
+                  
+      return
+      end
Index: inundation/ga/storm_surge/Hobart/MacDonald.bed
===================================================================
--- inundation/ga/storm_surge/Hobart/MacDonald.bed	(revision 1444)
+++ inundation/ga/storm_surge/Hobart/MacDonald.bed	(revision 1444)
@@ -0,0 +1,107 @@
+   -30.0         4.393
+   -20.0         4.393
+   -10.0         4.393
+     0.0         4.393      
+     10.0        4.316   
+     20.0        4.241   
+     30.0        4.168   
+     40.0        4.095   
+     50.0        4.024   
+     60.0        3.953   
+     70.0        3.884   
+     80.0        3.815   
+     90.0        3.747   
+     100.0       3.679   
+     110.0       3.611   
+     120.0       3.544   
+     130.0       3.477   
+     140.0       3.410   
+     150.0       3.343   
+     160.0       3.276   
+     170.0       3.209   
+     180.0       3.142   
+     190.0       3.074   
+     200.0       3.007   
+     210.0       2.938   
+     220.0       2.870   
+     230.0       2.801   
+     240.0       2.731   
+     250.0       2.661   
+     260.0       2.590   
+     270.0       2.519   
+     280.0       2.447   
+     290.0       2.375   
+     300.0       2.302   
+     310.0       2.228   
+     320.0       2.154   
+     330.0       2.079   
+     340.0       2.004   
+     350.0       1.928   
+     360.0       1.852   
+     370.0       1.776   
+     380.0       1.699   
+     390.0       1.623   
+     400.0       1.546   
+     410.0       1.469   
+     420.0       1.393   
+     430.0       1.317   
+     440.0       1.241   
+     450.0       1.165   
+     460.0       1.091   
+     470.0       1.017   
+     480.0       0.944   
+     490.0       0.872   
+     500.0       0.801   
+     510.0       0.737   
+     520.0       0.686   
+     530.0       0.646   
+     540.0       0.608   
+     550.0       0.574   
+     560.0       0.543   
+     570.0       0.514   
+     580.0       0.487   
+     590.0       0.462   
+     600.0       0.439   
+     610.0       0.418   
+     620.0       0.398   
+     630.0       0.379   
+     640.0       0.361   
+     650.0       0.345   
+     660.0       0.329   
+     670.0       0.314   
+     680.0       0.300   
+     690.0       0.287   
+     700.0       0.274   
+     710.0       0.262   
+     720.0       0.250   
+     730.0       0.239   
+     740.0       0.228   
+     750.0       0.218   
+     760.0       0.208   
+     770.0       0.198   
+     780.0       0.189   
+     790.0       0.179   
+     800.0       0.170   
+     810.0       0.161   
+     820.0       0.152   
+     830.0       0.144   
+     840.0       0.135   
+     850.0       0.127   
+     860.0       0.119   
+     870.0       0.110   
+     880.0       0.102   
+     890.0       0.094   
+     900.0       0.085   
+     910.0       0.077   
+     920.0       0.069   
+     930.0       0.061   
+     940.0       0.052   
+     950.0       0.044   
+     960.0       0.035   
+     970.0       0.027   
+     980.0       0.018   
+     990.0       0.009   
+     1000.0      0.000
+     1010.0      0.000
+     1020.0      0.000
+     1030.0      0.000   
