SUBROUTINE DO_U_MOMENTUM(t)
      implicit none
      REAL*8 t  !  going from n-1 to n; time at n is t

#include "SIZE.h"
#include "SWE.h"
#include "MASKS.h"
#include "CORIOLIS.h"

      integer i,j
      REAL*8 xu,yu,fu,fcor,Ravg
      REAL*8 calc_u_forcing

      do j=1,ny
         yu=ycor(j)  !  y for u grid
         do i=1,nx
            xu=(i-1)*dx  !  x for u grid
C           calculate forcing in momentum eqn. at x=xu,y=yu,time=(n-1/2)dt
            fu = calc_u_forcing(xu,yu,t-0.5*dt)
            Ravg = epsF
            fcor = avgU(i,j)*
     &            (f1v(j)*(maskV(i,j)*v(i,j)+maskV(i-1,j)*v(i-1,j)) 
     &            + (f2v(j)*(maskV(i,j+1)*v(i,j+1) 
     &              +maskV(i-1,j+1)*v(i-1,j+1))))
            u(i,j) = maskU(i,j)*
     &                (u(i,j)*(1-(0.5*dt*Ravg)) 
     &              - grav*alx*(eta(i,j)-eta(i-1,j)) 
     &              + dt*fu 
     &              + dt*fcor)/(1+(0.5*dt*Ravg))
         enddo
      enddo

      if (PeriodicInX) then
         do j=0,ny+1
C            eta(0,j)=eta(nx,j)
            u(nx+1,j)=u(1,j)  !  actual periodic BC
C            v(0,j)=v(nx,j)  !  for coriolis
         enddo
      endif
      if (PeriodicInY) then
         do i=0,nx+1           
C            eta(i,0)=eta(i,ny)
            u(i,0)=u(i,ny)  !  for coriolis
C            v(i,ny+1)=v(i,1)  !  actual periodic BC
         enddo
      endif

      END


      SUBROUTINE DO_V_MOMENTUM(t)
      implicit none
      REAL*8 t  !  going from n-1 to n; time at n is t

#include "SIZE.h"
#include "SWE.h"
#include "MASKS.h"
#include "CORIOLIS.h"

      integer i,j
      REAL*8 xv,yv,fv,fcor,Ravg
      REAL*8 calc_v_forcing

      do j=1,ny
         yv=ycor(j) - 0.5*dy  !  y for v grid
         do i=1,nx
            xv=(i-0.5)*dx  !  x for v grid
C           calculate forcing in momentum eqn. at x=xv,y=yv,time=(n-1/2)dt
            fv = calc_v_forcing(xv,yv,t-0.5*dt)
            Ravg = epsF
            fcor = -avgV(i,j)*
     &            (f1u(j)*(maskU(i+1,j-1)*u(i+1,j-1) 
     &              +maskU(i,j-1)*u(i,j-1)) 
     &            + (f2u(j)*(maskU(i+1,j)*u(i+1,j) 
     &              +maskU(i,j)*u(i,j))))
            v(i,j) = maskV(i,j)*
     &                (v(i,j)*(1-(0.5*dt*Ravg)) 
     &              - grav*aly*(eta(i,j)-eta(i,j-1)) 
     &              + dt*fv 
     &              + dt*fcor)/(1+(0.5*dt*Ravg))
         enddo
      enddo

      if (PeriodicInX) then
         do j=0,ny+1
C            eta(0,j)=eta(nx,j)
C            u(nx+1,j)=u(1,j)  !  actual periodic BC
            v(0,j)=v(nx,j)  !  for coriolis
         enddo
      endif
      if (PeriodicInY) then
         do i=0,nx+1           
C            eta(i,0)=eta(i,ny)
C            u(i,0)=u(i,ny)  !  for coriolis
            v(i,ny+1)=v(i,1)  !  actual periodic BC
         enddo
      endif

      END

      SUBROUTINE DO_CONTINUITY(t)
      implicit none
      REAL*8 t  !  going from n-1 to n; time at n is t

#include "SIZE.h"
#include "SWE.h"
#include "MASKS.h"

      integer i,j
C      REAL*8 Hup1,Hum1,Hvp1,Hvm1

      do j=1,ny
         do i=1,nx
C            Hup1 = min(H(i,j),H(i+1,j))
C            Hum1 = min(H(i,j),H(i-1,j))
C            Hvp1 = min(H(i,j),H(i,j+1))
C            Hvm1 = min(H(i,j),H(i,j-1))
C           update eta from time t=(n-1/2)dt to (n+1/2)dt
            eta(i,j) = eta(i,j) 
     &         - alx*(Hup1(i,j)*u(i+1,j)-Hum1(i,j)*u(i,j)) 
     &         - aly*(Hvp1(i,j)*v(i,j+1)-Hvm1(i,j)*v(i,j))
         enddo
      enddo

      if (PeriodicInX) then
         do j=0,ny+1
            eta(0,j)=eta(nx,j)
C            u(nx+1,j)=u(1,j)  !  actual periodic BC
C            v(0,j)=v(nx,j)  !  for coriolis
         enddo
      endif
      if (PeriodicInY) then
         do i=0,nx+1           
            eta(i,0)=eta(i,ny)
C            u(i,0)=u(i,ny)  !  for coriolis
C            v(i,ny+1)=v(i,1)  !  actual periodic BC
         enddo
      endif

      END