c-----------------------------------------------------------------------
c FORTRAN 77 routines for mike's cv code
c-----------------------------------------------------------------------
c Routine to read in grid properties and initial data parameters
c
                subroutine readinput(nx, ny, xmin, xmax, ymin, ymax,cfl,
     .                               niters, ml, thresh, re, mineff,
     .                               bufw, blkw, oe)
                        implicit none

                        integer nx, ny
                        real*8  xmin, xmax, ymin, ymax,
     .                          cfl, thresh,
     .                          mineff
                        integer niters, ml, re, bufw, blkw, oe, d, li
                        integer chk

c               Open file and read in parameters

                        open(unit=1,file='input.par')

                                read(1,*)nx, ny
                                read(1,*)xmin, xmax
                                read(1,*)ymin, ymax
                                read(1,*)cfl
                                read(1,*)niters
                                read(1,*)ml
                                read(1,*)thresh
                                read(1,*)re
                                read(1,*)mineff
                                read(1,*)bufw
                                read(1,*)blkw
                                read(1,*)oe

                        close(unit=1)

                        return
                end
c-----------------------------------------------------------------------

c-----------------------------------------------------------------------
c Routine to initialize data
c-----------------------------------------------------------------------
		subroutine initial(lb, ub, shape, u)
			implicit none

			integer	lb(2), ub(2), shape(2)

			real*8	u(shape(1), shape(2))

                        integer i, j
          do 20 i = 1, shape(1) 
          do 20 j = 1, shape(2) 
c          do 20 k = 1, shape(3) 
             u(i,j) = 0.
 20       continue

          if(lb(1).eq.0.and.lb(2).eq.0) then
             u(1,1) = 1.
          endif

		return
		end
c-----------------------------------------------------------------------

c-----------------------------------------------------------------------
c     Routine to update 
c-----------------------------------------------------------------------
      subroutine update(lb, ub, shape, u, fu, gu, hu, dt, h)
         implicit none

         integer  lb(2), ub(2), shape(2)

         real*8   u( shape(1), shape(2)),
     .            fu( shape(1), shape(2)),
     .            gu(  shape(1), shape(2)),
     .            hu(  shape(1), shape(2))

         real*8 unp1(shape(1),shape(2))
         integer  gshape(2)

        real*8   dt, h

        real ratio(3) 
        real fdir(3) 
        integer me
c
c
        integer i,j,k,ix,jy,xmax,ymax,xmin,ymin,zmin,zmax
        real bl
        real*8 flux
C
        ratio(1) = dt/h
        ratio(2) = dt/h

        xmin = 1
        ymin = 1
        zmin = 1
C
         xmax = shape(1)
         ymax = shape(2)
         zmax = 1

C        define flux assuming simple upwind
C 

C        write (*,*) gshape(1), gshape(2)
C        write (*,*) shape(1), shape(2)
C        write (*,*) me,':',xmin,ymin,xmax,ymax
       
          do 20 i = xmin, xmax 
          do 20 j = ymin, ymax 
             unp1(i,j) = u(i,j)
 20       continue
c
          do 30 i = xmin, xmax-1 
          do 30 j = ymin, ymax
c
c assume +ve wave
             fu(i,j) = ratio(1) * flux ( u(i,j) )
             if (lb(1).eq.0 .and. i .eq. 1) then
                bl = 0.
             else
                bl = 1.
             endif
             unp1(i,j) = - bl * fu(i,j) + unp1(i,j) 
             unp1(i+1,j) = fu(i,j) + unp1(i+1,j) 
 30       continue
c
          i = xmax
          do 35 j = ymin, ymax
c
c assume +ve wave
             fu(i,j) = ratio(1) * flux ( u(i,j) )
             unp1(i,j) = - fu(i,j) + unp1(i,j) 
 35       continue
c
          do 40 i = xmin, xmax
          do 40 j = ymin, ymax-1 
c
c assume +ve wave
             gu(i,j) = ratio(2) * flux ( u(i,j) )
             if (lb(2).eq.0 .and. j .eq. 1) then
                bl = 0.
             else
                bl = 1.
             endif
             unp1(i,j) = - bl * gu(i,j) + unp1(i,j) 
             unp1(i,j+1) = gu(i,j) + unp1(i,j+1) 
 40       continue
c
          j = ymax 
          do 45 i = xmin, xmax 
c
c assume +ve wave
             gu(i,j) = ratio(2) * flux ( u(i,j) )
             unp1(i,j) = - gu(i,j) + unp1(i,j) 
 45       continue
c
        do 70 i = xmin, xmax
        do 70 j = ymin, ymax
        do 70 k = zmin, zmax 
c
c inflow b.c.
          if ( lb(1) .eq. 0 .and. lb(2) .eq. 0 .and.
     .         i .eq. 1 .and. j .eq. 1 ) go to 70
c
           u(i,j) = unp1(i,j)
 70     continue
        
c
        return 
	end
c
        real*8 function flux ( u )

         real*8 u, u2
c        flux = u

        u2 = u * u
        flux = u2 /( u2 + ( 1. - u ) ** 2 )

        return
        end

c-----------------------------------------------------------------------

c-----------------------------------------------------------------------
c     Routine to update boundary
c-----------------------------------------------------------------------
      subroutine update_boundary(lb, ub, shape, u, fu, gu, hu, dt, h)
         implicit none

         integer  lb(2), ub(2), shape(2)

         real*8   u( shape(1), shape(2)),
     .            fu( shape(1), shape(2)),
     .            gu(  shape(1), shape(2)),
     .            hu(  shape(1), shape(2))

        real*8   dt, h
c
c
        integer i,j,k,ix,jy,xmax,ymax,xmin,ymin,zmin,zmax
        real bl
        real*8 flux
C
        if (lb(1).eq.0 .and. lb(2).eq.0) then

           u(1,1) = 1.

        endif
c
        return
        end

c-----------------------------------------------------------------------

c-----------------------------------------------------------------------
      subroutine pythnorm(lb, ub, shape, truncU, truncerr,
     .                    wgtu)
c-----------------------------------------------------------------------
        implicit none

        integer lb(2), ub(2), shape(2)

        real*8  truncU(shape(1), shape(2)),
     .          truncerr(shape(1), shape(2))

        real*8  wgtu

c       Local variables

        integer i, j, k

        do i = 1, shape(1)
          do j = 1, shape(2)
            truncerr(i,j) = sqrt( wgtu * truncU(i,j)**2 )
          end do
        end do

        return
       end
c-----------------------------------------------------------------------
