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  = 0.5d0
      Tout   = 1.0d0
      tf = 1500.0d0
      length = 1000.0d0

      nlow = md+2
      nhigh = md+nx
      nxx = nx+2*md
      
      limith = 1
      
      dx = 10000.d0/dble(nx-1)
      K = 0.0d0            

      if (nx.gt.nxmax) write(*,*)'nx too large'
      toll = 1.0d-10
      cfl = 0.2d0
      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)
 
         d1h = h1 - h0
         d2h = h2 - h1
         d1w = w1 - w0
         d2w = w2 - w1
         dz  = (z2(i) - z2(i-1))
                
c ------- W = H+Z           
c         h_x = xmin( d1w, d2w) - dz
         h_x = xmic(theta, d1h, d2h)
         
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
c$$$         hmin = dmin1(h0,h1,h2)
c$$$         if ( h_l(i) .le. hmin ) then
c$$$            h_l(i) = hmin
c$$$            h_r(i) = h(i) + (h(i) - h_l(i))
c$$$         elseif ( h_r(i) .le. hmin ) then
c$$$            h_r(i) = hmin
c$$$            h_l(i) = h(i) + (h(i) - h_r(i))
c$$$         endif
      enddo
c
c ---- Enforce a minimum height
c
c$$$      do i=md,nx+md+2
c$$$         hmax = dmax1(h_l(i),h_r(i))
c$$$         if (hmax .lt. 1.0d-3 ) then
c$$$c           h_r(i)  = 0.0d0
c$$$           uh_r(i) = 0.0d0
c$$$c           h_l(i)  = 0.0d0
c$$$           uh_l(i) = 0.0d0
c$$$c           h(i)    = 0.0d0
c$$$           uh(i)   = 0.0d0
c$$$         endif
c$$$      enddo
c
c ---- Try for a flat subinterval 
cc
c      do i=md,nx+md+2
c         hmax = dmax1(h_l(i),h_r(i))
c         dz  = (z2(i) - z2(i-1))
c
c         if (h_l(i) .lt. 0.01d0*h_r(i) ) then
c           h_r(i)  = sqrt(abs(2.0d0*dz*h(i)))
c         endif
c
c         if (h_r(i) .lt. 0.01d0*h_l(i) ) then         
c           h_l(i)  = sqrt(abs(2.0d0*dz*h(i)))       
c         endif         
c      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 ) 
c          u_x = xmin( 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