	program ParabolicWave3

	real*8 grid(150,250,800)
	real*8 load(150,250)

	real*8 coords(4)

	integer allow(150,250)
	integer shape(2)
	real*8 dx,dy,dt,Lx,Ly,Lto,t,x,y
	real*8 e
	integer i,j,k,nx,ny,nt,rank
	
	e=2.718d0
	rank=2


	Lx=10.0d0
	Ly=15.0d0
	Lto=14.0d0

	coords(1)=0.0d0
	coords(2)=Lx
	coords(3)=-Ly/2
	coords(4)=Ly/2

	nx=150
	ny=250
	nt=800

	shape(1)=nx
	shape(2)=ny


	dx=Lx/nx
	dy=Ly/ny
	dt=Lto/nt

C loading logical array

	do i=1,nx

	  x=i*dx

	  do j=1,ny

	    y=j*dy-5

	    if(y .lt. -1.0/4*(x-5)**2+10 .or. y .lt. 0) then
	        allow(i,j)= 1
	    endif


	  enddo
	enddo

C done loading logical

C loading initial waveform with zero initial velocity
	
	do i=1,nx

	   do j=1,ny

	    if(allow(i,j) .eq. 1) then
	     grid(i,j,1)=exp(-10*(j*dy-5)**2)
	     grid(i,j,2)=grid(i,j,1)+dt*2*10*(j*dy-5)
     &	*grid(i,j,1)
	    endif
	   enddo
	enddo


C done loading initial waveform

C boundary conditions are all zero

C now for main loop

	do i=2,nt-1
	   do k=2,ny-1
	    if(allow(j,k) .eq. 1) then
	     grid(1,k,i+1)=2*grid(1,k,i)-grid(1,k,i-1)
     &	+((dt/dx)**2)*(grid(2,k,i)-2*grid(1,k,i)+grid(nx,k,i))
     &	+((dt/dy)**2)*(grid(1,k+1,i)-2*grid(1,k,i)+grid(1,k-1,i))
	    endif

	     do j=2,nx-1
	      if(allow(j,k) .eq. 1) then
	        grid(j,k,i+1)=2*grid(j,k,i)-grid(j,k,i-1)
     &	+((dt/dx)**2)*(grid(j+1,k,i)-2*grid(j,k,i)+grid(j-1,k,i))
     &	+((dt/dy)**2)*(grid(j,k+1,i)-2*grid(j,k,i)+grid(j,k-1,i))
	      endif
	     enddo
	 if(allow(j,k) .eq. 1) then
	grid(nx,k,i+1)=2*grid(nx,k,i)-grid(nx,k,i-1)
     &	+((dt/dx)**2)*(grid(1,k,i)-2*grid(nx,k,i)+grid(nx-1,k,i))
     &	+((dt/dy)**2)*(grid(nx,k+1,i)-2*grid(nx,k,i)+grid(nx,k-1,i))
	 endif
	   enddo
	   write(*,*) i
	enddo

C main body done.  Now load to file


	do i=1,nt
	  t=(i-1)*dt
	   do j=1,nx
	      do k=1,ny
	        load(j,k)=grid(j,k,i)
	      enddo
	   enddo
	   if(mod(i,2).eq.0)then
	  call gft_out_bbox('parabola3',t,shape,rank,coords,load)
	   endif
	  write(*,*) i
	enddo

	write(*,*) (((dt/dx)**2)+((dt/dy)**2))**(0.5)

	stop
	end
