! source file: /Users/oschlies/UVIC/master/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
!-----------------------------------------------------------------------

      include "param.h"
      include "emode.h"
      include "grdvar.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"
      include "tmngr.h"
      dimension 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