      program gork

c     all unit si

c     size of the ocean in the east-west direction
      parameter (a=2000e3)

c     size of the ocean in the north-south direction
      parameter (b=2000e3)

c     wind is ws*cos(pi*x/a) in N-S direction
      parameter (ws=1.0,pi=3.1416)

c     beta, f parameters
      parameter (f=1e-4,beta=1.6e-11)

c     water depth, Ekman layer depth, density
      parameter (H=4000,delta=40,rho=1000.)

c     two parameters in equation
      parameter (Be=2*beta*H/(f*delta),gamma=-2*pi*ws/(rho*delta*f*a))

c     number of grid points in x
      parameter (npts=41)

c     gridsize (set to be same in x and y)
      parameter (d=a/(npts-1.0))

c     number of grid points in y
      parameter (mpts=b/d+1)

c     solve for the streamfunction
      real*8 psi(npts,mpts)

c     initialize the array to zero, on boundaries values are never changed

      do i=1,npts
         do j=1,mpts
            psi(i,j) = 0
         enddo
      enddo

c     ask user for number of iteraction to make (500 is good)

      read (*,*) nrelax

c     iterate, solving the problem by relaxation
      do k=1,nrelax

c     loop over domain
      do i=2,npts-1
         do j=2,mpts-1
            psi(i,j) = -gamma*sin(pi*(i-1.)/(npts-1.))*d*d/4.+
     >           (psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1))/4.
     >           +Be*(psi(i+1,j)-psi(i-1,j))*d/8.
         enddo
      enddo
      enddo

c     write result to file gork.dat in format suitable for gnplot
c     eg type gnuplot on command line and then
c     > splot 'gork.dat'
c     for contour plot
c     > set nosurface
c     > set view 0,0
c     > set contour
c     > set data style lines

      open (11,file='gork.dat')
      do j=1,mpts
         do i=1,npts
            write (11,*) psi(i,j)
         enddo
         write (11,*)
      enddo
      close (11)

      end
