! source file: /Users/csomes/Research/Models/UVic_ESCM/2.9/source/mom/calc_psgrad.F
      subroutine calc_psgrad (psgrad, uext, vext, js, je, is, ie)

!-----------------------------------------------------------------------
!     compute the surface pressure gradients

!     inputs:

!     js   = index of starting row
!     je   = index of ending row
!     is   = index of starting longitude
!     ie   = index of ending longitude

!     outputs:

!     psgrad  = grad(surf press)
!     uext = external mode u (tau+1) for point (ie,je) only
!     vext = external mode v (tau+1) for point (ie,je) only
!-----------------------------------------------------------------------

      implicit none

      integer is, ie, js, je, js1, je1, is1, ie1, jrow, i, kz

      real fxa, r2dtuv, f3, atosp, f2, uext, vext, d1, d2, diag1, diag0
      real diag3, diag4, dubdt, dvbdt

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "emode.h"
      include "grdvar.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"
      include "tmngr.h"

      real psgrad(is:ie,js:je,2)

      js1 = max(js,2)
      je1 = min(je,jmt-1)
      is1 = max(is,2)
      ie1 = min(ie,imtm1)

!     on mixing time steps "ptd" has been multiplied by a factor of
!     two and the time step has to be adjusted also.

      if (mod(itt,nmix) .eq. 1) then
        fxa = p5
      else
        fxa = c1
      endif
      r2dtuv = c1/c2dtuv
!     f3    = c2dtuv/c2dtsf
      f3    = c2dtuv
      do jrow=js1,je1
        do i=is1,ie1
          atosp = acor*cori(i,jrow,1)
          f2    = atosp*c2dtuv
          kz = kmu(i,jrow)
          if (kz .ne. 0) then
            diag1        = psi(i+1,jrow+1,1)-psi(i  ,jrow,1)
            diag0        = psi(i  ,jrow+1,1)-psi(i+1,jrow,1)
            uext         = -(diag1+diag0)*dyu2r(jrow)*hr(i,jrow)
            vext         =  (diag1-diag0)*dxu2r(i)*hr(i,jrow)*csur(jrow)
            diag3        = fxa*(ptd(i+1,jrow+1)-ptd(i  ,jrow))
            diag4        = fxa*(ptd(i  ,jrow+1)-ptd(i+1,jrow))
            dubdt        = (diag3+diag4)*dyu2r(jrow)*hr(i,jrow)
            dvbdt        = (diag3-diag4)*dxu2r(i)*hr(i,jrow)*csur(jrow)
            psgrad(i,jrow,1)=r2dtuv*(dubdt + f3*zu(i,jrow,1) + f2*dvbdt)
            psgrad(i,jrow,2)=r2dtuv*(-dvbdt+ f3*zu(i,jrow,2) + f2*dubdt)
          else
            psgrad(i,jrow,1) = c0
            psgrad(i,jrow,2) = c0
          endif
        enddo
      enddo

      return
      end