program main
C*************************************************************
      implicit none
#include "SIZE.h"
#include "SWE.h"
#include "CORIOLIS.h"
c=========================================
c declaration
c=========================================
      REAL*8 dumarray(NXM)

C     read and set some constants
      CALL INI_PARAMS('data.par')
      nout=int(tout/dt) 
C     read initial conditions and H array
      CALL INI_FIELDS
      CALL INI_MASKS  !  initialize masks
      CALL APPLY_MASKS ! applies masks to u,v
      CALL INI_CORIOLIS
      CALL READ_FORCING
      CALL INI_FORCING

      write(*,*)'Model parameters are:'
      write(*,*)'tmax=',tmax,' dt=',dt,' nsteps= ',nsteps
      write(*,*)'tout=',tout,' nout= ',nout
      write(*,*)'dx=',dx,' xmax=',xmax,' nx=',nx
      write(*,*)'dy=',dy,' ymax=',ymax,' ny=',ny
      write(*,*)'epsF=',epsF
      write(*,*)'fo=',fo, ' beta=',beta,' y0=',y0
      write(*,*)'PeriodicInX: ',PeriodicInX
      write(*,*)'PeriodicInY: ',PeriodicInY

C     This is to catch even/odd coriolis time stepping 'problem'.
C     If this is a pickup then expect tstart not equal to 0.
      if (mod(int(tstart/dt),2).NE.0) then
        write(*,*)'WARNING: tstart/dt is NOT even'
      endif

      CALL MODEL

      end

C****************************************************************
 
C****************************************************************
C GRID: C grid in space. u and eta are staggered in time

C x=0          dx                    (i-1)dx       idx        (i+1)dx            Ndx 
C   |     +     |     + ... |     +     |     +     |     +     | ... |     +     |
C  u(0)        u(1)                  u(i-1)        u(i)       u(i+1)            u(N)
C       eta(1)      eta(2)     eta(i-1)     eta(i)      eta(i+1)           eta(N)
C
C water depth H(i), dissipation R(i), forcing f(i,n) are on u grid
C t=0            dt                      (n-1)dt         ndt           (n+1)dt
C  |       +      |      + ...  |     +      |      +     |       +      |
C u(0)          u(1)                      u(n-1)         u(n)          u(n+1) 
C       eta(1/2)     eta(1+1/2)                eta(n-1/2)     eta(n+1/2)  

      SUBROUTINE MODEL
      implicit none
#include "SIZE.h"
#include "SWE.h"

      integer n,iRec
      REAL*8 t

C-------initialize io and write out initial condition
      t=0.d0
      n=0
      iRec=1

      write(*,*)'Writing output at time ',t
      call write_r8_field(nxm+2,nym+2,iRec,u(0,0),uout)
      call write_r8_field(nxm+2,nym+2,iRec,v(0,0),vout)
      call write_r8_field(nxm+2,nym+2,iRec,eta(0,0),etaout)
      iRec=iRec+1
      
C-------loop over time steps       
      do n=1,nsteps
	 t = tstart + dt*n  ! going from n-1 to n; time at n is t

C        Toggle momentum time-stepping for coriolis
	 if (mod(n,2).eq.0) then
	    CALL DO_U_MOMENTUM(t)
	    CALL DO_V_MOMENTUM(t)
	 else
	    CALL DO_V_MOMENTUM(t)
	    CALL DO_U_MOMENTUM(t)
	 endif
	 CALL DO_CONTINUITY(t)
                    
         if (mod(n,nout).eq.0) then   ! if time to write out
            write(*,*)'Writing output at time ',t
            call write_r8_field(nxm+2,nym+2,iRec,u(0,0),uout)
            call write_r8_field(nxm+2,nym+2,iRec,v(0,0),vout)
            call write_r8_field(nxm+2,nym+2,iRec,eta(0,0),etaout)
            iRec=iRec+1
         endif		
      enddo  !  end loop 

      END