c-----------------------------------------------------------------------
c program1.f
c PHY387N Spring 96
c
c Example program to illustrate
c  o Compilation and execution of programs 
c  o How to read and write from and to files from a program
c  o How to plot data using gnuplot (After the program is run)
c
c Program sets up a 2D grid function and dumps it to a file.
c-----------------------------------------------------------------------
		program first

		implicit none

		integer nmax, nx, ny
      parameter(nmax=100)

c		Declare arrays
      real*8	u(nmax,nmax)
		real*8   xx(nmax), yy(nmax)

c		Parameters
		real*8	hx, hy, max(2), min(2)

c		Open input file for parameters

		open(unit=11,file='input.dat')
		  read(11,*)nx, ny 
		  read(11,*)min(1), max(1)
		  read(11,*)min(2), max(2)
      close(unit=11)

		write(6,*)'Grid size ',nx,' by ', ny
		write(6,*)'min/max in x direction',min(1),max(1)
		write(6,*)'min/max in y direction',min(2),max(2)

c		Calculate hx, hy

		hx = (max(1) - min(1))/float(nx)
		hy = (max(2) - min(2))/float(ny)

c		Generate coordinates

		write(6,*)'Generating coordinates '
		call gencrds(nx, ny, xx, yy, min, max, hx, hy)

c		Calculate simple function over x and y

		write(6,*)'Calculating simple function'
		call getfunc(nx, ny, xx, yy, u)

c		Dump the data

		write(6,*)'Write data to file'
		call dumpdat(nx, ny, xx, yy, u, 'output.dat')

		stop
		end
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c	subroutine to generate coordinates
		subroutine gencrds(nx, ny, xx, yy, min, max, hx, hy)
		  implicit none
       
		  integer	nx, ny
   
        real*8    hx, hy,
     *            min(2), max(2),
     *            xx(nx), yy(ny)

c     Local variables
		  integer i

		  do i = 1, nx
		    xx(i) = min(1) + float(i-1) * hx
		  end do

		  do i = 1, ny
		    yy(i) = min(2) + float(i-1) * hy
		  end do

        return
		end 
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c	subroutine to set up simple function 
		subroutine getfunc(nx, ny, xx, yy, u)
		  implicit none

		  integer 	nx, ny

		  real*8		xx(nx), yy(ny), u(nx,ny)

c		Local variables

		  integer i, j

		  do i = 1, nx
          do j = 1, ny
		      u( i, j) =  sin(0.5 * (xx(i)**2 + yy(j)**2) )
          end do
		  end do

		  return
		end
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c	subroutine to dump data to output file
		subroutine dumpdat(nx, ny, xx, yy, u, fname)
		  implicit none

		  integer	nx, ny

		  real*8		xx(nx), yy(ny), u(nx,ny)
		  character*32	fname

c		Local Variables
		  integer 	i, j
		

c		Open output file
		  open(unit=12,file=fname)

c     Write out data
		  do i = 1, nx
		    do j = 1, ny

		      write(12,'(3F20.15)')xx(i),yy(j),u(i,j)

		    end do
c			Put in empty line for parametric surface plot by gnuplot
			 write(12,*)
		  end do

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