PROGRAM MAIN

!     Declaring parameters associated with sizes

      INTEGER IX,IY            ! Loop Variables
      INTEGER NXC,NYC          ! Array length read variables
      INTEGER NX,NY            ! Array size variables
      PARAMETER(NX=34,NY=17)   ! Setting array sizes

!     Declaring Arrays required to read data
      
      REAL HT(NX,NY)           ! Array for initial height field
      REAL PSI(NX,NY)          ! Array for initial streamfunction field 
      REAL U(NX,NY)            ! Array for initial U field
      REAL V(NX,NY)            ! Array for initial V field


!     General parameters
      REAL DX                  ! Grid size parameter
      REAL G                   ! Acceleration due to gravity
      REAL F0                  ! Coriolis Parameter
      PARAMETER(G=9.81,F0=1.0E-4)

!     START OF PROGRAM CODE

!     FIRST THE OPEN STATEMENT ASSOCIATES THE DATA FILE WITH YOUR
!     FORTRAN PROGRAM      

      OPEN(10,FILE='msc11.dat')
      READ (10,100) NXC,NYC,DX  !Header of file contains dimension info
 100  FORMAT(2I4,E15.7)         !Format statement for read
                                
!     NXC and NYC provide a check - they should agree with the values
!     of NX and NY defined in your program

!     SIMPLE ERROR LOOP

      IF ((NXC.NE.NX).OR.(NYC.NE.NY)) THEN
         PRINT*,'Mismatch in array sizes - program aborted'
         STOP
      ENDIF   

!     READ STATEMENTS
      DO IY=1,NY
         READ (10,101) (HT(IX,IY),IX=1,NX) !Read data 
      END DO
 101  FORMAT(5E15.7)

      CLOSE(10)

!     END OF READ LOOP - DATA CONTAINED IN HT ARRAY


!     CALCULATE STREAMFUNCTION

      DO IY=1,NY
         DO IX=1,NX
            PSI(IX,IY)=HT(IX,IY)*(G/F0)
         END DO
      END DO


!     CALCULATE VELOCITY VECTORS

      DO IX=2,NX-1
         DO IY=2,NY-1
            U(IX,IY)=-((PSI(IX,IY+1)-PSI(IX,IY-1))/(2.0*DX))
            V(IX,IY)=(PSI(IX+1,IY)-PSI(IX-1,IY))/(2.0*DX)
         ENDDO   
      ENDDO


!     Your Program Code Should Go Here Including 
!     Including PSOLVE etc.
!     Sample outputting code is included below


!     SAMPLE OUTPUT LOOP - THIS OUTPUTS READ HT DATA TO A FILE TO BE
!     READ BY MATLAB


      OPEN(5,FILE='ht.dat')
      DO IX=1,NX
         DO IY=1,NY
            WRITE(5,*) HT(IX,IY)
         ENDDO   
      ENDDO   
      CLOSE(5)

      OPEN(5,FILE='psi.dat')
      DO IX=1,NX
         DO IY=1,NY
            WRITE(5,*) PSI(IX,IY)
         ENDDO   
      ENDDO   
      CLOSE(5)

      OPEN(5,FILE='u.dat')
      DO IX=1,NX
         DO IY=1,NY
            WRITE(5,*) U(IX,IY)
         ENDDO   
      ENDDO   
      CLOSE(5)

      OPEN(5,FILE='v.dat')
      DO IX=1,NX
         DO IY=1,NY
            WRITE(5,*) V(IX,IY)
         ENDDO   
      ENDDO   
      CLOSE(5)
      END