! source file: /Users/csomes/Research/Models/UVic_ESCM/2.9/source/mom/tropic.F
      subroutine tropic (c2dtsf, acor, f, itt, dtts)

!=======================================================================

!        S O L V E   T H E   B A R O T R O P I C   E Q U A T I O N

!   There are several significant changes made in MOM_2 in the
!   calculation of the vertically averaged velocities.

!   1. Based on a 1994 rederivation of the finite difference equations
!     for the stream function by Charles Goldberg (GFDL/Princeton
!     University/Trenton State College), the coefficients in
!     the Poisson equations for dpsi differ slighty from those
!     used in MOM_1.  Designed to reduce the residuals in the
!     finite difference momentum equations, these new equations
!     for dpsi seem to be less stiff.  Tests on a variety of
!     geometries and topographies indicate that the new equations
!     converge more rapidly to the same tolerances in all solvers
!     at all resolutions.

!   2. In all three Poisson solvers for stream function, the MOM_2
!     default mode is that the values of psi and dpsi are not held
!     fixed at zero on the boundary of the "main continent".  They
!     are allowed to float in the same way that other "island
!     boundary" values float.  Although this requires an island
!     integral to be computed on the boundary of every land mass, it
!     also makes the iterative system less stiff, and again results
!     in fewer iterations to converge to the same tolerances
!     in all solvers at all resolutions.

!     Tests indicate that except on architectures where island
!     integrals are prohibitively expensive and the stream function
!     formulation therefore computationally infeasible, the cost
!     of the extra island integrals is small and the savings due
!     to reduced numbers of iterations significant.

!     The user may specify either a land mass on whose boundary the
!     stream function is later normalized to zero or that no post
!     solver normalization is to take place.  Surface pressures are
!     always normalized to have mean zero.  Options are also provided
!     for turning off the island integrals on the boundaries of
!     selected land masses; however, this practice is not recommended.

!   3. The convergence criterion in all three solvers has been
!     changed to the following:  "stop when the predicted maximum error
!     in the solved variable (dpsi or surface pressure) is less than
!     the user specified tolerance."  Convergence tolerances are now in
!     the same units as the variable solved for, and tell the user how
!     many digits of the answer are correct.Thus, if one expects a
!     maximum dpsi of 1.0e12 and desires convergence to 5 significants
!     digits, one chooses a tolerance of 1.0e7. Note that the tolerance
!     used in MOM 2 is NOT the same as "crit" in MOM 1.

!     The maximum error in the solution is predicted as follows:
!     First, convergence of the solver is assumed to be "geometric" in
!     the sense that the maximum absolute correction added to the
!     solved variable in iteration k is modeled as

!           step(k) = step(1)*(convergence_rate)**(k-1)

!     The estimated maximum error in the solution after k iterations
!     is then bounded by the sum of the missing terms in the geometric
!     series truncated after k terms:

!     sum {step(i)} = step(k)*convergence_rate/(1.0 - convergence_rate)
!       i=k+1,infinity

!     Experimental evidence indicates that the convergence rate
!     of an iterative solver remains essentially stable over many
!     iterations, and that the maximum errors when the solvers are
!     stopped by this criterion as compared to solutions obtained
!     by allowing the solvers to run to machine precision are indeed
!     less than the stated tolerances.
!=======================================================================

      implicit none

      integer i, jrow, itt, nconv, luptdb, luptd, npt

      real c2dtsf, dtts, acor, fxa, dpsi1, absmax

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "emode.h"
      include "grdvar.h"
      include "iounit.h"
      include "mw.h"
      include "switch.h"

      character(8) :: bc_symm

      real f(imt,jmt)

      data nconv /0/
      save nconv

      if (eoyear) nconv = 0

      call setbcx (zu(1,1,1), imt, jmt)
      call setbcx (zu(1,1,2), imt, jmt)

!-----------------------------------------------------------------------
!     construct the forcing for the stream function equation
!-----------------------------------------------------------------------

      call sfforc (zu, dxu, dyu, csu, ztd)

!     filter forcing at high latitudes

      call filz (ztd, cf)

!-----------------------------------------------------------------------
!     read in solution from 2 previous timesteps for the
!     purpose of computing an initial guess for the present solution.
!-----------------------------------------------------------------------

      luptdb = nkflds - mod(itt,2)
      luptd  = nkflds - 1 + mod(itt,2)
      call oget (kflds, nwds, luptdb, res)
      call oget (kflds, nwds, luptd, ptd)

      fxa=c1
      if (.not. leapfrog) fxa=p5
      do jrow=1,jmt
        do i=1,imt
          ptd(i,jrow) = fxa*(c2*ptd(i,jrow)-res(i,jrow))
        enddo
      enddo
      do jrow=2,jmtm1
        ptd(1,jrow)   = ptd(imtm1,jrow)
        ptd(imt,jrow) = ptd(2,jrow)
      enddo

!-----------------------------------------------------------------------
!     choose 5 or 9 point numerics
!-----------------------------------------------------------------------

!     initialize coefficients using 5 point numerics

      call sfc5pt (acor, f, c2dtsf, dxu, dyu, csu, hr, cf)
      npt = 5

!-----------------------------------------------------------------------
!     choose a method for solving for the "tau+1" stream function change
!-----------------------------------------------------------------------

      variable   = 'dpsi'
      bc_symm    = 't odd'

      call congr (npt, variable, bc_symm, ptd, ptd, ztd, res
     &,           cf
     &,           mxscan, mscan, tolrsf
     &,           imask, iperm, jperm, iofs, nisle, nippts
     &,           converged, esterr)

!     correct for drifting dpsi on land mass "imain"

      if (imain .gt. 0) then
        dpsi1 = ptd(iperm(iofs(imain)+1), jperm(iofs(imain)+1))
        call con_adjust (ptd, dpsi1, map)
      endif

!-----------------------------------------------------------------------
!     test accuracy of solving for change in stream function
!-----------------------------------------------------------------------

      if (.not.converged) then
        write (stdout,'(a,i5,3(a,1pe10.3))')
     &  ' WARNING: SOLVER DID NOT CONVERGE in ',mscan
     &, ' scans. max(psi)='
     &, absmax(psi(1,1,2)), ' max(dpsi)=',absmax(ptd)
     &, ' estimated max(err)=', esterr
        nconv = nconv + 1
        if (nconv .gt. 50) stop 'nconv > 50 in tropic.f'
      endif

!-----------------------------------------------------------------------
!     update the stream function based upon the solution
!-----------------------------------------------------------------------

      if (euler2) then
        do jrow=1,jmt
          do i=1,imt
            psi(i,jrow,1) = psi(i,jrow,2) + ptd(i,jrow)
          enddo
        enddo
      else
        do jrow=1,jmt
          do i=1,imt
            res(i,jrow)    = psi(i,jrow,2) + ptd(i,jrow)
            psi(i,jrow,2)  = psi(i,jrow,1)
            psi(i,jrow,1)  = res(i,jrow)
          enddo
        enddo
      endif

!-----------------------------------------------------------------------
!     save ptd to compute 1st guess for relaxation next timestep
!     (..note.. on 1st pass of euler backward timestep, bypass this
!            save, since it will be done on the 2nd pass)
!     (..note.. on a mixing timestep, alter ptd to be consistent with
!            normal, leap-frog stepping)
!-----------------------------------------------------------------------

      if (.not. euler1) then

        if (.not. leapfrog) then
          do jrow=1,jmt
            do i=1,imt
              ptd(i,jrow)=c2*ptd(i,jrow)
            enddo
          enddo
        endif

        call oput (kflds, nwds, luptdb, ptd)

      endif

      return
      end

      subroutine sfforc (zu, dxu, dyu, csu, forc)

!=======================================================================

!           S T R E A M   F U N C T I O N   F O R C I N G
!=======================================================================

      implicit none

      real cddxu(0:1,0:1), cddyu(0:1,0:1)
      real cddxt(-1:0,-1:0), cddyt(-1:0,-1:0), p5, c0
      parameter (p5=0.5, c0=0.0)

      include "size.h"

      integer i, jrow, i1, j1

      real zu(imt,jmt,2), dxu(imt), dyu(jmt), csu(jmt)
      real forc(imt,jmt), ustuff(imt,jmt), vstuff(imt,jmt)

!-----------------------------------------------------------------------
!     initialize the forcing
!-----------------------------------------------------------------------

      do i=1,imt
        do jrow=1,jmt
          forc(i,jrow) = c0
        enddo
      enddo

!-----------------------------------------------------------------------
!     construct coefficients for partial differences. a partial
!     difference in "x" is defined as an "x" difference of a quantity
!     which is averaged in "y". (and symmetrically for "y" differences).
!     Note that this is an x difference and NOT an x derivitive.
!     partial differences of quantities on the "t" grid are defined on
!     the "u" grid and visa versa.
!     therefore partial differences at:
!     u/v points (i,j), involve nearby t/s points with subscripts:
!        (i  ,j+1)    (i+1,j+1)
!        (i  ,j  )    (i+1,j  )
!     t/s points (i,j), involve nearby u/v points with subscripts:
!        (i-1,j  )    (i  ,j  )
!        (i-1,j-1)    (i  ,j-1)
!     thus if qu(i,j) is defined on u/v points, its partial
!     difference ddxqt = ddxt(qu) is defined on t/s points and has the
!     value
!     ddxqt(i,j) = cddxt(-1,-1)*qu(i-1,j-1) + cddxt(-1,0)*qu(i-1,j+0)
!                + cddxt( 0,-1)*qu(i+0,j-1) + cddxt( 0,0)*qu(i+0,j+0)
!-----------------------------------------------------------------------

      cddxu( 0, 0) = -p5
      cddxu( 0, 1) = -p5
      cddxu( 1, 0) =  p5
      cddxu( 1, 1) =  p5

      cddxt(-1,-1) = -p5
      cddxt(-1, 0) = -p5
      cddxt( 0,-1) =  p5
      cddxt( 0, 0) =  p5

      cddyu( 0, 0) = -p5
      cddyu( 0, 1) =  p5
      cddyu( 1, 0) = -p5
      cddyu( 1, 1) =  p5

      cddyt(-1,-1) = -p5
      cddyt(-1, 0) =  p5
      cddyt( 0,-1) = -p5
      cddyt( 0, 0) =  p5

!-----------------------------------------------------------------------
!     multiply the u eqn by dx*cos, the v eqn by dy, then subtract their
!     partial differences to eliminate the unknown surface pressure from
!     the resulting equation
!-----------------------------------------------------------------------

      do i=1,imt-1
        do jrow=1,jmt-1
          ustuff(i,jrow) = zu(i,jrow,1)*dxu(i)*csu(jrow)
          vstuff(i,jrow) = zu(i,jrow,2)*dyu(jrow)
        enddo
      enddo

      do i1=-1,0
        do j1=-1,0
          do jrow=2,jmt-1
            do i=2,imt-1
              forc(i,jrow) = forc(i,jrow)
     &                     - cddyt(i1,j1)*ustuff(i+i1,jrow+j1)
     &                     + cddxt(i1,j1)*vstuff(i+i1,jrow+j1)
            enddo
          enddo
        enddo
      enddo

      return
      end

      subroutine sfc5pt (acor, f, c2dtsf, dxu, dyu, csu, hr, coef)

!=======================================================================

!     5  P T    C O E F F I C I E N T   I N I T I A L I A Z A T I O N

!     coefficient initialization for 5 point elliptic solvers

!     inputs:

!     acor   = implicit coriolis factor (0.0 => 1.0)
!     f      = 2*omega*sin(phi(j))
!     c2dtsf = twice the time step (seconds)
!     dxu    = width of "u" grid cell (cm)
!     dyu    = height of "u" grid cell (cm)
!     csu    = cosine of "u" grid cell
!     hr     = 1/depth at "u" cells (cm)

!     outputs:

!     coeff   = 3 x 3 array of coefficients at each (i,j) point
!=======================================================================

      implicit none

      integer jj, ii, j, i, i1, j1, i2, j2

      real p5, c0, c2dtsf, acor
      parameter (p5=0.5, c0=0.0)
      real cddxu(0:1,0:1), cddyu(0:1,0:1)
      real cddxt(-1:0,-1:0), cddyt(-1:0,-1:0)

      include "size.h"

      real csu(jmt), dxu(imt), dyu(jmt), hr(imt,jmt)
      real f(imt,jmt), coef(imt,jmt,-1:1,-1:1)
      real ustuff(imt,jmt), vstuff(imt,jmt)

!-----------------------------------------------------------------------
!     initialize the coefficients
!-----------------------------------------------------------------------

      do jj=-1,1
        do ii=-1,1
          do j=1,jmt
            do i=1,imt
              coef(i,j,ii,jj) = c0
            enddo
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     construct coefficients for partial differences. a partial
!     difference in "x" is defined as an "x" difference of a quantity
!     which is averaged in "y". (and symmetrically for "y" differences).
!     Note that this is an x difference and NOT an x derivitive.
!     partial differences of quantities on the "t" grid are defined on
!     the "u" grid and visa versa.
!     therefore partial differences at:
!     u/v points (i,j), involve nearby t/s points with subscripts:
!        (i  ,j+1)    (i+1,j+1)
!        (i  ,j  )    (i+1,j  )
!     t/s points (i,j), involve nearby u/v points with subscripts:
!        (i-1,j  )    (i  ,j  )
!        (i-1,j-1)    (i  ,j-1)
!     thus if qu(i,j) is defined on u/v points, its partial
!     difference ddxqt = ddxt(qu) is defined on t/s points and has the
!     value
!     ddxqt(i,j) = cddxt(-1,-1)*qu(i-1,j-1) + cddxt(-1,0)*qu(i-1,j+0)
!                + cddxt( 0,-1)*qu(i+0,j-1) + cddxt( 0,0)*qu(i+0,j+0)
!-----------------------------------------------------------------------

      cddxu( 0, 0) = -p5
      cddxu( 0, 1) = -p5
      cddxu( 1, 0) =  p5
      cddxu( 1, 1) =  p5

      cddxt(-1,-1) = -p5
      cddxt(-1, 0) = -p5
      cddxt( 0,-1) =  p5
      cddxt( 0, 0) =  p5

      cddyu( 0, 0) = -p5
      cddyu( 0, 1) =  p5
      cddyu( 1, 0) = -p5
      cddyu( 1, 1) =  p5

      cddyt(-1,-1) = -p5
      cddyt(-1, 0) =  p5
      cddyt( 0,-1) = -p5
      cddyt( 0, 0) =  p5

!-----------------------------------------------------------------------
!     compute coefficients for all points
!-----------------------------------------------------------------------

      do i=1,imt-1
        do j=1,jmt-1
          ustuff(i,j) = dxu(i)*csu(j)*hr(i,j) / (c2dtsf*dyu(j))
          vstuff(i,j) = dyu(j)*hr(i,j) / (c2dtsf*dxu(i)*csu(j))
        enddo
      enddo

!-----------------------------------------------------------------------
!     calculate 5 point coefficients

!     note that ne [and nw] coefficient adds to n coefficient in
!     ustuff term, but ne [and se] coefficient adds to e coefficient in
!     vstuff term for the 5 point operator.
!-----------------------------------------------------------------------

      do i1=0,1
        do j1=0,1
          do i2=-1,0
            do j2=-1,0
              do j=2,jmt-1
                do  i=2,imt-1
                  coef(i,j,0,j1+j2) = coef(i,j,0,j1+j2) +
     &               cddyu(i1,j1)*cddyt(i2,j2)*ustuff(i+i2,j+j2)
                  coef(i,j,i1+i2,0) = coef(i,j,i1+i2,0) +
     &               cddxu(i1,j1)*cddxt(i2,j2)*vstuff(i+i2,j+j2)
                enddo
              enddo
            enddo
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     augment coefficients for implicit treatment of coriolis term
!     all coefficients are calculated, but corner ones are zero.
!-----------------------------------------------------------------------

      if (acor .ne. 0.0) then
        do i=1,imt-1
          do j=1,jmt-1
            ustuff(i,j) = acor*hr(i,j)*(-f(i,j))
            vstuff(i,j) = acor*hr(i,j)*( f(i,j))
          enddo
        enddo
        do i1=0,1
          do j1=0,1
            do i2=-1,0
              do j2=-1,0
                do j=2,jmt-1
                  do  i=2,imt-1
                    coef(i,j,i1+i2,j1+j2) = coef(i,j,i1+i2,j1+j2)
     &                 - cddxu(i1,j1)*cddyt(i2,j2)*ustuff(i+i2,j+j2)
                    coef(i,j,i1+i2,j1+j2) = coef(i,j,i1+i2,j1+j2)
     &                 - cddyu(i1,j1)*cddxt(i2,j2)*vstuff(i+i2,j+j2)
                  enddo
                enddo
              enddo
            enddo
          enddo
        enddo
      endif

      return
      end

      subroutine sfc9pt  (acor, f, c2dtsf, dxu, dyu, csu, hr, coef)

!=======================================================================

!     9  P T    C O E F F I C I E N T   I N I T I A L I A Z A T I O N

!     coefficient initialization for 9 point elliptic solvers

!     inputs:

!     acor   = implicit coriolis factor (0.0 => 1.0)
!     f      = 2*omega*sin(phi(j))
!     c2dtsf = twice the time step (seconds)
!     dxu    = width of "u" grid cell (cm)
!     dyu    = height of "u" grid cell (cm)
!     csu    = cosine of "u" grid cell
!     hr     = 1/depth at "u" cells (cm)

!     outputs:

!     coeff   = 3 x 3 array of coefficients at each (i,j) point
!=======================================================================

      implicit none

      integer jj, ii, j, i, i1, j1, i2, j2

      real c0, p5, c2dtsf, acor
      parameter (c0=0.0, p5=0.5)
      real cddxu(0:1,0:1), cddyu(0:1,0:1)
      real cddxt(-1:0,-1:0), cddyt(-1:0,-1:0)

      include "size.h"

      real csu(jmt), dxu(imt), dyu(jmt), hr(imt,jmt)
      real f(imt,jmt)
      real coef(imt,jmt,-1:1,-1:1)
      real ustuff(imt,jmt), vstuff(imt,jmt)

!-----------------------------------------------------------------------
!     initialize the work area
!-----------------------------------------------------------------------

      do jj=-1,1
        do ii=-1,1
          do j=1,jmt
            do i=1,imt
              coef(i,j,ii,jj) = c0
            enddo
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     generate arrays of coefficients

!     construct coefficients for partial differences. a partial
!     difference in "x" is defined as an "x" difference of a quantity
!     which is averaged in "y". (and symmetrically for "y" differences).
!     Note that this is an x difference and NOT an x derivitive.
!     partial differences of quantities on the "t" grid are defined on
!     the "u" grid and visa versa.
!     therefore partial differences at:
!     u/v points (i,j), involve nearby t/s points with subscripts:
!        (i  ,j+1)    (i+1,j+1)
!        (i  ,j  )    (i+1,j  )
!     t/s points (i,j), involve nearby u/v points with subscripts:
!        (i-1,j  )    (i  ,j  )
!        (i-1,j-1)    (i  ,j-1)
!     thus if qu(i,j) is defined on u/v points, its partial
!     difference ddxqt = ddxt(qu) is defined on t/s points and has the
!     value
!     ddxqt(i,j) = cddxt(-1,-1)*qu(i-1,j-1) + cddxt(-1,0)*qu(i-1,j+0)
!                + cddxt( 0,-1)*qu(i+0,j-1) + cddxt( 0,0)*qu(i+0,j+0)
!-----------------------------------------------------------------------

      cddxu( 0, 0) = -p5
      cddxu( 0, 1) = -p5
      cddxu( 1, 0) =  p5
      cddxu( 1, 1) =  p5

      cddxt(-1,-1) = -p5
      cddxt(-1, 0) = -p5
      cddxt( 0,-1) =  p5
      cddxt( 0, 0) =  p5

      cddyu( 0, 0) = -p5
      cddyu( 0, 1) =  p5
      cddyu( 1, 0) = -p5
      cddyu( 1, 1) =  p5

      cddyt(-1,-1) = -p5
      cddyt(-1, 0) =  p5
      cddyt( 0,-1) = -p5
      cddyt( 0, 0) =  p5

!-----------------------------------------------------------------------
!     compute coefficients for all points
!-----------------------------------------------------------------------

      do i=1,imt-1
        do j=1,jmt-1
          ustuff(i,j) = dxu(i)*csu(j)*hr(i,j) / (c2dtsf*dyu(j))
          vstuff(i,j) = dyu(j)*hr(i,j) / (c2dtsf*dxu(i)*csu(j))
        enddo
      enddo

!---------------------------------------------------------------------
!     calculate 9 point coefficients
!---------------------------------------------------------------------

      do i1=0,1
        do j1=0,1
          do i2=-1,0
            do j2=-1,0
              do j=2,jmt-1
                do  i=2,imt-1
                  coef(i,j,i1+i2,j1+j2) = coef(i,j,i1+i2,j1+j2) +
     &               cddyu(i1,j1)*cddyt(i2,j2)*ustuff(i+i2,j+j2)  +
     &               cddxu(i1,j1)*cddxt(i2,j2)*vstuff(i+i2,j+j2)
                enddo
              enddo
            enddo
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     augment coefficients for implicit treatment of coriolis term
!     all coefficients are calculated, but corner ones are zero.
!-----------------------------------------------------------------------

      if (acor .ne. 0.0) then
        do i=1,imt-1
          do j=1,jmt-1
            ustuff(i,j) = acor*hr(i,j)*(-f(i,j))
            vstuff(i,j) = acor*hr(i,j)*( f(i,j))
          enddo
        enddo
        do i1=0,1
          do j1=0,1
            do i2=-1,0
              do j2=-1,0
                do j=2,jmt-1
                  do  i=2,imt-1
                    coef(i,j,i1+i2,j1+j2) = coef(i,j,i1+i2,j1+j2)
     &               - cddxu(i1,j1)*cddyt(i2,j2)*ustuff(i+i2,j+j2)
     &               - cddyu(i1,j1)*cddxt(i2,j2)*vstuff(i+i2,j+j2)
                  enddo
                enddo
              enddo
            enddo
          enddo
        enddo
      endif

      return
      end

      subroutine spforc (zu, dxu, dyu, csu, h, forc)

!=======================================================================
!           S U R F A C E   P R E S S U R E   F O R C I N G
!=======================================================================

      implicit none

      integer i, j, i1, j1

      real p5
      parameter (p5=0.5)
      real cddxu(0:1,0:1), cddyu(0:1,0:1)
      real cddxt(-1:0,-1:0), cddyt(-1:0,-1:0)

      include "size.h"

      real zu(imt,jmt,2), dxu(imt), dyu(jmt), csu(jmt), h(imt,jmt)
      real forc(imt,jmt), ustuff(imt,jmt), vstuff(imt,jmt)

!-----------------------------------------------------------------------
!     generate arrays of coefficients

!     construct coefficients for partial differences. a partial
!     difference in "x" is defined as an "x" difference of a quantity
!     which is averaged in "y". (and symmetrically for "y" differences).
!     Note that this is an x difference and NOT an x derivitive.
!     partial differences of quantities on the "t" grid are defined on
!     the "u" grid and visa versa.
!     therefore partial differences at:
!     u/v points (i,j), involve nearby t/s points with subscripts:
!        (i  ,j+1)    (i+1,j+1)
!        (i  ,j  )    (i+1,j  )
!     t/s points (i,j), involve nearby u/v points with subscripts:
!        (i-1,j  )    (i  ,j  )
!        (i-1,j-1)    (i  ,j-1)
!     thus if qu(i,j) is defined on u/v points, its partial
!     difference ddxqt = ddxt(qu) is defined on t/s points and has the
!     value
!     ddxqt(i,j) = cddxt(-1,-1)*qu(i-1,j-1) + cddxt(-1,0)*qu(i-1,j+0)
!                + cddxt( 0,-1)*qu(i+0,j-1) + cddxt( 0,0)*qu(i+0,j+0)
!-----------------------------------------------------------------------

      cddxu( 0, 0) = -p5
      cddxu( 0, 1) = -p5
      cddxu( 1, 0) =  p5
      cddxu( 1, 1) =  p5

      cddxt(-1,-1) = -p5
      cddxt(-1, 0) = -p5
      cddxt( 0,-1) =  p5
      cddxt( 0, 0) =  p5

      cddyu( 0, 0) = -p5
      cddyu( 0, 1) =  p5
      cddyu( 1, 0) = -p5
      cddyu( 1, 1) =  p5

      cddyt(-1,-1) = -p5
      cddyt(-1, 0) =  p5
      cddyt( 0,-1) = -p5
      cddyt( 0, 0) =  p5

!     weight "zu" and "zv" by the cell area and take the divergence

      do i=1,imt-1
        do j=1,jmt-1
          ustuff(i,j) = h(i,j)*zu(i,j,1)*dyu(j)
          vstuff(i,j) = h(i,j)*zu(i,j,2)*dxu(i)*csu(j)
        enddo
      enddo

      do i=1,imt
        do j=1,jmt
          forc(i,j) = 0.0
        enddo
      enddo

      do i1=-1,0
        do j1=-1,0
          do i=2,imt-1
            do j=2,jmt-1
              forc(i,j) = forc(i,j) + cddxt(i1,j1)*ustuff(i+i1,j+j1)
     &                              + cddyt(i1,j1)*vstuff(i+i1,j+j1)
            enddo
          enddo
        enddo
      enddo

      return
      end

      subroutine spc9pt (dxu, dyu, csu, h, coef)

!=======================================================================

!     S U R F A C E   P R E S S U R E    C O E F F I C I E N T

!                I N I T I A L I A Z A T I O N

!     inputs:

!     dxu    = width of "u" grid cell (cm)
!     dyu    = height of "u" grid cell (cm)
!     csu    = cosine of "u" grid cell
!     h      = depth at "u,v" cells (cm)

!     outputs:

!     coeff   = 3 x 3 array of coefficients at each (i,j) point
!=======================================================================

      implicit none

      integer i1, j1, i, j, i2, j2

      real c0, p5
      parameter (c0=0.0, p5=0.5)
      real cddxu(0:1,0:1), cddyu(0:1,0:1)
      real cddxt(-1:0,-1:0), cddyt(-1:0,-1:0)

      include "size.h"

      real csu(jmt), dxu(imt), dyu(jmt), h(imt,jmt)
      real coef(imt,jmt,-1:1,-1:1)
      real ustuff(imt,jmt), vstuff(imt,jmt)

!-----------------------------------------------------------------------
!     generate arrays of coefficients

!     construct coefficients for partial differences. a partial
!     difference in "x" is defined as an "x" difference of a quantity
!     which is averaged in "y". (and symmetrically for "y" differences).
!     Note that this is an x difference and NOT an x derivitive.
!     partial differences of quantities on the "t" grid are defined on
!     the "u" grid and visa versa.
!     therefore partial differences at:
!     u/v points (i,j), involve nearby t/s points with subscripts:
!        (i  ,j+1)    (i+1,j+1)
!        (i  ,j  )    (i+1,j  )
!     t/s points (i,j), involve nearby u/v points with subscripts:
!        (i-1,j  )    (i  ,j  )
!        (i-1,j-1)    (i  ,j-1)
!     thus if qu(i,j) is defined on u/v points, its partial
!     difference ddxqt = ddxt(qu) is defined on t/s points and has the
!     value
!     ddxqt(i,j) = cddxt(-1,-1)*qu(i-1,j-1) + cddxt(-1,0)*qu(i-1,j+0)
!                + cddxt( 0,-1)*qu(i+0,j-1) + cddxt( 0,0)*qu(i+0,j+0)
!-----------------------------------------------------------------------

      cddxu( 0, 0) = -p5
      cddxu( 0, 1) = -p5
      cddxu( 1, 0) =  p5
      cddxu( 1, 1) =  p5

      cddxt(-1,-1) = -p5
      cddxt(-1, 0) = -p5
      cddxt( 0,-1) =  p5
      cddxt( 0, 0) =  p5

      cddyu( 0, 0) = -p5
      cddyu( 0, 1) =  p5
      cddyu( 1, 0) = -p5
      cddyu( 1, 1) =  p5

      cddyt(-1,-1) = -p5
      cddyt(-1, 0) =  p5
      cddyt( 0,-1) = -p5
      cddyt( 0, 0) =  p5

!-----------------------------------------------------------------------
!     compute coefficients for all points
!-----------------------------------------------------------------------

!     initialize all 9 coefficients to zero

      do i1=-1,1
        do j1=-1,1
          do i=1,imt
            do j=1,jmt
              coef(i,j,i1,j1) = c0
            enddo
          enddo
        enddo
      enddo

      do j=1,jmt
        do i=1,imt
          ustuff(i,j) = 0.0
          vstuff(i,j) = 0.0
        enddo
      enddo
      do i=1,imt-1
        do j=1,jmt-1
          ustuff(i,j) = h(i,j)*dyu(j)/(dxu(i)*csu(j))
          vstuff(i,j) = h(i,j)*dxu(i)*csu(j)/dyu(j)
        enddo
      enddo

!     calculate divergence = ddx (ddx (ustuff)) + ddy( ddy (vstuff))

      do i1=0,1
        do j1=0,1
          do i2=-1,0
            do j2=-1,0
              do i=2,imt-1
                do j=2,jmt-1
                  coef(i,j,i1+i2,j1+j2) = coef(i,j,i1+i2,j1+j2)
     &                 + cddxu(i1,j1) * cddxt(i2,j2) * ustuff(i+i2,j+j2)
     &                 + cddyu(i1,j1) * cddyt(i2,j2) * vstuff(i+i2,j+j2)
                enddo
              enddo
            enddo
          enddo
        enddo
      enddo

      return
      end

      subroutine filz (fext, cf)

!=======================================================================
!     subroutine filz sets up input needed for fourier filtering
!     (when the "fourfil" ifdef is defined) -or- symmetric finite
!     impulse response filtering (when the "firfil" ifdef is defined)
!     of "fext" at the specified high latitudes. "fext" is forcing for
!     the external mode.
!=======================================================================

      implicit none

      integer jrow, jj, l, is, ie, ii, i, im, m, n, isv, iev, j, num

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "emode.h"
      include "grdvar.h"
      include "index.h"
      include "levind.h"

      real fext(imt,jmt)
      real cf(imt,jmt,3)
      real temp(imt)

!=======================================================================

      do jrow=jfrst,jmtm1
        if (jrow.le.jft1 .or. jrow.ge.jft2) then

!---------------------------------------------------------------------
!         fourier filter fext at high latitudes
!---------------------------------------------------------------------

            jj = jrow - jfrst + 1
            if (jrow .ge. jft2) jj = jj - jskpt + 1
            do l=1,lsegf
              is = iszf(jj,l)
              if (is .ne. 0) then
                ie = iezf(jj,l)
                do ii=is,ie
                  i = mod(ii-2,imtm2) + 2
                  temp(ii+1-is) = fext(i,jrow)
                enddo
                im = ie-is+1
                if (im .ne. imtm2) then
                   m = 1
                   n = nint(im*cst(jrow)*cstr(jft0))
                else
                   m = 3
                   n = nint(im*cst(jrow)*cstr(jft0)*p5)
                endif

                call filtr (temp(1), im, m ,n, 0)

                do ii=is,ie
                  i = mod(ii-2,imtm2)+2
                  fext(i,jrow) = temp(ii+1-is)
                enddo
              endif
            enddo
        endif
      enddo

      return
      end