c-----------------------------------------------------------------------
c FORTRAN 77 routines for 2D wave equation solution
c-----------------------------------------------------------------------
c Routine to read in grid properties and initial data parameters
c
		subroutine readinput(nx, ny, xmin, xmax, ymin, ymax,cfl,
     .		  amplitude, width, vx, vy, niters, ml, thresh)
			implicit none

			integer	nx, ny
         integer niters, ml

			real*8	xmin, xmax, ymin, ymax,
     .		     cfl, amplitude, width, vx, vy, thresh

c		Open file and read in parameters

			open(unit=1,file='wave2d.in')

				read(1,*)nx, ny
				read(1,*)xmin, xmax
				read(1,*)ymin, ymax
				read(1,*)width
				read(1,*)amplitude
				read(1,*)vx, vy
				read(1,*)cfl
				read(1,*)niters
				read(1,*)ml
				read(1,*)thresh
			
			close(unit=1)

			return
		end 
c-----------------------------------------------------------------------
c Routine to initialize data
		subroutine initial(lb, ub, shape, phi0, A0, B0, C0,
     .				 h, dt, ampl, width, vx, vy, xmin, ymin)
			implicit none

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

			real*8	phi0(shape(1), shape(2)), 
     .				A0(shape(1), shape(2)), 
     .				B0(shape(1), shape(2)), 
     .				C0(shape(1), shape(2))

			real*8	h, dt, ampl, width, vx, vy, xmin, ymin

c		Local variables

			integer	i, j

			real*8	x, y, t

c     Initialize a0, b0, c0 based on what phi0 and dtphi0 are
c     Carry out simple minded initial data for dtphi0 for now and
c     Gaussian wave-packet for phi0

         t = 0.0

         do i = 1, shape(1)
           do j = 1, shape(2)
             x = h * (lb(1)/2 +i-1) + xmin 
             y = h * (lb(2)/2 +j-1) + ymin

             phi0(i,j) = ampl * exp(-((x - vx*t)**2 +
     .                  (y - vy*t)**2)/width**2)

           end do
         end do
c     Initialize a0, b0, c0
         do i = 1, shape(1)
           do j = 1, shape(2)
             x = h * (lb(1)/2  + i-1) + xmin
             y = h * (lb(2)/2  + j-1) + ymin

             a0(i,j) = 2.0 * ((x - vx*t)*vx + (y-vy*t)*vy)*
     .                  phi0(i,j)/width**2

             b0(i,j) = -2.0 * (x-vx*t) * phi0(i,j) / width**2


             c0(i,j) = -2.0 * (y-vy*t) * phi0(i,j) / width**2
           end do
         end do

			return
		end
c-----------------------------------------------------------------------
c Interior Update using the Leap frog scheme
c
c     anp => A at n+1
c     anm => A at n-1
c     an  => A at n
c
c     bnp => B at n+1
c     bnm => B at n-1
c     bn  => B at n
c
c     cnp => C at n+1
c     cnm => C at n-1
c     cn  => C at n
      subroutine leapfrog2d(lb, ub, shape, anp, anm, an,
     .            bnp, bnm, bn, cnp, cnm, cn, phip, phim, phi,
     .            dt, h)
         implicit none

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


         real*8   anp( shape(1), shape(2)),
     .            anm( shape(1), shape(2)),
     .            an(  shape(1), shape(2)),
     .            bnp( shape(1), shape(2)),
     .            bnm( shape(1), shape(2)),
     .            bn(  shape(1), shape(2)),
     .            cnp( shape(1), shape(2)),
     .            cnm( shape(1), shape(2)),
     .            cn(  shape(1), shape(2)),
     .            phip(shape(1), shape(2)),
     .            phim(shape(1), shape(2)),
     .            phi( shape(1), shape(2))

         real*8   dt, h, epsilon

c        Local variables

         integer  i, j, s_sten
c        Loop over interior grid points

	 s_sten = 1
         epsilon = 0.25

c           update A
         do i = s_sten+1, shape(1)-s_sten
            do j = s_sten+1, shape(2)-s_sten
               anp(i,j) = anm(i,j) +
     .            (dt /h) * (bn(i+1,j) - bn(i-1,j)) +
     .            (dt /h) * (cn(i,j+1) - cn(i,j-1)) 

            end do
         end do

c           update B
         do i = s_sten+1, shape(1)-s_sten
            do j = s_sten+1, shape(2)-s_sten
               bnp(i,j) = bnm(i,j) +
     .            (dt / h) * (an(i+1,j) - an(i-1,j)) 
            end do
         end do

c           update C
         do i = s_sten+1, shape(1)-s_sten
            do j = s_sten+1, shape(2)-s_sten
               cnp(i,j) = cnm(i,j) +
     .            (dt / h) * (an(i,j+1) - an(i,j-1)) 

            end do
         end do

         do i = s_sten+1, shape(1) - s_sten
            do j = s_sten+1, shape(2) - s_sten
               phip(i,j) = phim(i,j) + 2.0 * dt * an(i,j)
            end do
         end do

		

         return
      end
c-----------------------------------------------------------------------
c     Routine to update phi directly using 2nd order PDE form
      subroutine update2nd(lb, ub, shape, phip,phim,phi, dt, h)
         implicit none

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

         real*8   phip( shape(1), shape(2)),
     .            phim( shape(1), shape(2)),
     .            phi(  shape(1), shape(2))

         real*8   dt, h

c     Local variables

         integer  i, j

         do i = 2, shape(1) - 1
            do j = 2, shape(2) - 1
               phip(i,j) = phim(i,j)
     .   +     (dt/h)**2 * (phi(i+1,j) + phi(i-1,j) - 2.0*phi(i,j))
     .   +     (dt/h)**2 * (phi(i,j+1) + phi(i,j+1) - 2.0*phi(i,j))
            end do
         end do

         return
      end
c-----------------------------------------------------------------------
      subroutine updatephi(lb, ub, shape, phip,phim,an, dt, h)
         implicit none

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

         real*8   phip( shape(1), shape(2)),
     .            phim( shape(1), shape(2)),
     .            an(  shape(1), shape(2))

         real*8   dt, h

c     Local variables

         integer  i, j

         do i = 2, shape(1) - 1
            do j = 2, shape(2) - 1
               phip(i,j) = phim(i,j) + 2.0 * dt * an(i,j)
            end do
         end do

         return
      end
c-----------------------------------------------------------------------
c  Boundary update routine
c  For now set simple boundary condition

      subroutine bdyupdt2d(lb, ub, shape, a, b, c, phi,
     .				dt, h, glb, gub)
			implicit none

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


         real*8   a( shape(1), shape(2)), 
     .            b( shape(1), shape(2)), 
     .            c( shape(1), shape(2)), 
     .				phi( shape(1), shape(2))

         real*8   dt, h

c        Local variables

         integer  i, j, k

c		Possible combinations 

c       lb(1) = glb(1) then we are at the x=xmin boundary

			if( lb(1) .eq. glb(1) )then
            do j = 1, shape(2)
			     a(lb(1), j) = 0.0		     
			     b(lb(1), j) = 0.0		     
			     c(lb(1), j) = 0.0		     

			     phi(lb(1), j) = 0.0		     
			   end do
			end if

c       ub(1) = gub(1) then we are at the x=xmax boundary

			if( ub(1) .eq. gub(1) )then
            do j = 1, shape(2)
			     a(ub(1), j) = 0.0		     
			     b(ub(1), j) = 0.0		     
			     c(ub(1), j) = 0.0		     

			     phi(ub(1), j) = 0.0		     
			   end do
			end if

c       lb(2) = glb(2) then we are at the y=ymin boundary

			if( lb(2) .eq. glb(2) )then
            do i = 1, shape(1)
			     a(i,lb(2)) = 0.0		     
			     b(i,lb(2)) = 0.0		     
			     c(i,lb(2)) = 0.0		     

			     phi(i, lb(2)) = 0.0		     
			   end do
			end if

c       ub(2) = gub(2) then we are at the y=ymax boundary

			if( ub(2) .eq. gub(2) )then
            do i = 1, shape(1)
			     a(i, ub(2)) = 0.0		     
			     b(i, ub(2)) = 0.0		     
			     c(i, ub(2)) = 0.0		     

			     phi(i, ub(2)) = 0.0		     
			   end do
			end if

         return
      end
c-----------------------------------------------------------------------
      subroutine pythnorm(lb, ub, shape, truncA,truncB,truncC,truncerr,
     .           wgta, wgtb, wgtc)
	implicit none

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

	real*8	trunca(shape(1), shape(2)), 
     .		truncb(shape(1), shape(2)), 
     .		truncc(shape(1), shape(2)), 
     .		truncerr(shape(1), shape(2))

	real*8	wgta, wgtb, wgtc

c	Local variables

        integer i, j, k

c	write(0,*)'weights ',wgtA, wgtB, wgtC

        do i = 1, shape(1)
          do j = 1, shape(2)
            truncerr(i,j) = sqrt( wgtA * trunca(i,j)**2 +
     .            wgtB * truncB(i,j)**2 + wgtC * truncC(i,j)**2 )
          end do
        end do

c	do i = 1, shape(1)
c          do j = 1, shape(2)
c             write(12,*)i,j,truncerr(i,j)
c          end do
c          write(12,*)
c        end do

        return
       end
c-----------------------------------------------------------------------
c  Routine to do smoothing in 2d
c   Take four neighboring points (i,j) 
c
c   (i,j)  weight of w. The rest are a weight of (1-w)/4
c          o    
c
c      o   o   o
c
c          o
c
       subroutine smooth2(lb, ub, shape, uf,w)
          implicit none

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

        real*8  uf(shape(1), shape(2)), w

c     Local variables

        integer i, j


c     Smoothing using four points in the interior only


        do i = 2, shape(1)-1
          do j = 2, shape(2)-1
            uf(i,j) =  w*uf(i,j) + (1.0d0 - w) * 0.25d0 * 
     .            (uf(i+1, j) + uf(i-1, j) + uf(i, j+1) + uf(i, j-1) )
          end do
        end do

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