! source file: /usr/local/models/UVic_ESCM/2.9/source/mom/congrad.F
      subroutine congr  (npt, variable, bc_symm
     &,                  guess, dpsi, forc, res
     &,                  cf
     &,                  max_iterations, iterations, epsilon
     &,                  imask, iperm, jperm, iofs, nisle, nippts
     &,                  converged
     &,                  estimated_error
     &                  )

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

!                            C O N G R A D

!      solve:

!             A * dpsi = forc

!      for "dpsi" with dirichlet boundary conditions (dpsi=const on
!      each component of the boundary) by a preconditioned conjugate
!      gradient algorithm.

!      inputs:
!              npt   = 5 or 9 (active coefficients)
!              variable = character string identifying solution variable
!              bc_symm = equatorial symmetry type (used only when the
!                        symmetry option is on. otherwise ignore it)
!              guess = initial approximation to solution
!              A     = linear operator (assumed symmetric)
!                      typically A is  grad{(1/h)*grad(dpsi)} -
!                      2dt*acor*{grad(f/h) x grad(dpsi)}
!                      using 5 or 9 pt discretizations
!              cf    = imt x jmt x 3 x 3 array of coefficients of A
!              forc  = the sum of all terms evaluated at times tau
!                      or tau-1
!              epsilon = convergence criterion
!              max_iterations = maximum number of iterations
!              imask = shows which land masses have perimeter equations
!              iperm = i coordinate of island perimeter points
!              jperm = j coordinate of island perimeter points
!              iofs  = offset in iperm, jperm for start of perimeter
!                      of land_mass(isle)
!              nisle = actual number of land_masses
!              nippts = number of perimeter ocean points for a land_mass
!      output:
!              dpsi   = answer
!              iterations = actual number of iterations performed
!              converged = logical value
!              estimated_error = estimated maximum error in solution
!                          based on step sizes and convergence rate

!      based on the preconditioned conjugate gradient algorithm given
!      in:

!     A Reformulation and Implementation of the Bryan-Cox-Semtner
!     Ocean Model on the Connection Machine
!     J.K. Dukowicz, R.D. Smith, and R.C. Malone
!     Journal of Atmospheric and Oceanic Technology
!     Vol 10. No. 2 April 1993

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

!      more specifically, the equations to be solved are

!             sum (A(ij,i'j') * dpsi(i'j')) = forc(ij)

!      where the subscripts ij and i'j' range over all "free ocean"
!      T cells ij=(i,j) that are not adjacent to land T cells,
!      and one ij=isle for each boundary component of the ocean.

!      with this choice of variables, in the absence of coriolis terms
!      (acor=0), the operator A is symmetric, i.e.,

!             A(ij,i'j') = A(i'j',ij)

!      the algorithm (essentially executable in Fortran 90) is...

!      subroutine congrad (A, guess, forc, dpsi, iterations)

!      use matrix_module

!      intent (in)     :: A, guess, forc
!      intent (out)    :: dpsi, iterations

!      type(dpsi_type) :: guess, dpsi, Zres, s
!      type(res_type)  :: res, As, forc
!      type(operator)  :: A
!      type(inv_op)    :: Z
!      dimension (0:max_iterations) :: dpsi, res, s, As, beta, alpha

!      dpsi(0) = guess
!      res(0)  = forc - A * dpsi(0)
!      beta(0) = 1
!      s(0)    = zerovector()
!      do k = 1 to max_iterations
!        Zres(k-1) = Z * res(k-1)
!        beta(k)   = res(k-1) * Zres(k-1)
!        s(k)      = Zres(k-1) + (beta(k)/beta(k-1)) * s(k-1)
!        As(k)     = A * s(k)
!        alpha(k)  = beta(k) / (s(k) * As(k))
!        dpsi(k)   = dpsi(k-1) + alpha(k) * s(k)
!        res(k)    = res(k-1) - alpha(k) * As(k)
!        estimated_error = err_est(k, alpha(k), s(k))
!        if (estimated_error) < epsilon) exit
!      enddo
!      if (k > max_iterations) then
!        print *, 'did not converge in ',k,' iterations'
!        stop '=>congrad'
!      endif

!      iterations = k
!      dpsi = dpsi(k)

!      end

!      where...

!      the "vector" and "operator" types used in conjugate gradient
!      are mapped to ordinary 2-dimensional fortran arrays as follows:

!      type(dpsi_type) :: guess, dpsi, Zres, s
!          if ij=(i,j) is a mid-ocean point, map dpsi(ij)-->dpsi(i,j)
!          if ij=isle is an ocean boundary subscript, replicate the
!          value dpsi(isle) in dpsi(i,j) for each (i,j) in the ocean
!          perimeter of land_mass(isle).  the arrays iperm(isle) and
!          jperm(isle), along with iofs(isle) locate these ocean
!          perimeter T cells.
!      type(res_type)  :: res, As, forc
!          if ij=(i,j) is a mid-ocean point, res(ij)-->res(i,j)
!          if ij=isle is an ocean boundary subscript, the value of
!          res(isle) = sum (res(i,j))
!          where the sum is taken over all (i,j) in the ocean perimeter
!          of land_mass(isle).  sometimes, the computed values
!          res(i,j) represent contributions of T cell (i,j) to the
!          component res(isle), and sometimes the values are balanced
!          so that res(i,j)=res(isle)/nippts(isle).  note that, even
!          when balanced, the relation between type(res_type) variables
!          res(isle) and res(i,j) differs from that of type(dpsi_type)
!          variables dpsi(isle) and dpsi(i,j) on T cells in the ocean
!          perimeter.
!      type(operator)  :: A
!          the nearly diagonal quality of the operators used
!          permits a representation as a small collection of
!          2-dimensional arrays.
!          the diagonal, A(ij,ij), is stored in an array cfdiag(i,j)
!          as follows:
!          if ij=(i,j) is a mid-ocean point, A(ij,ij) = cfdiag(i,j)
!          if ij=isle is an ocean boundary subscript,
!          A(isle, isle) = sum (cfdiag(i,j))
!          where the sum is taken over all (i,j) in the ocean perimeter
!          of land_mass(isle).  each cfdiag(i,j) represents the contribution
!          of T cell (i,j) to the island variable diagonal coefficient.
!          the off-diagonal terms A(ij,i`j`) are stored in 4 arrays
!          cfn, cfs, cfe, and cfw if A is a 5-point operator, and in
!          these and 4 additional arrays, cfne, cfnw, cfse, cfsw, if
!          A is a 9-point operator.  For example, if i`=i and j`=j+1,
!          then A(ij,i`j`) is stored in cfn(i,j).
!          if ij=(i,j) is a mid-ocean point and i`j`=isle` is and ocean
!          perimeter subscript, with i`=i and j`=j+1, then
!          cfn(i,j)=A(ij,isle`) is the coefficient of the island
!          variable dpsi(isle`) in the equation for mid-ocean point
!          dpsi(ij)=dpsi(i,j).
!          if ij=isle is an ocean perimeter point and i`j`=(i`,j`) is
!          a mid-ocean point, with i`=i and j`=j-1, then
!          cfs(i,j)=A(isle,i`j`) is the coefficient of the mid-ocean
!          variable dpsi(i`j`)=dpsi(i,j) in the equation for the island
!          variable dpsi(isle).  note that equations for island
!          variables dpsi(isle) are "non-local" in the sense that
!          they usually contain more than 5 or 9 terms, some of which
!          involve values dpsi(i`j`) outside of a compact 5-point
!          or 9-point neighborhood.
!      type(inv_op)    :: Z
!          the approximate inverse operator Z used at present is a
!          diagonal operator Z(ij,ij) = 1/A(ij,ij).
!          if ij=(i,j) is a mid-ocean point,
!            then Z(i,j)=Z(ij)=1/A(ij)=1/cfdiag(i,j)
!          if ij=isle is an ocean perimeter point, then
!          Z(isle) is replicated at each ocean perimeter T cell
!          bordering land_mass(isle).
!            Z(i,j)=Z(isle)=1/A(isle)=1/sum(A(i,j))

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

      implicit none

      character(16) :: variable
      character(*) :: bc_symm

      integer nerror, j, nisle, npt, k, imax, jmax, max_iterations
      integer is, ie, js, je, js1, je1, is1, ie1, jrow, i, kz

      logical converged, diverging

      real zresmax, absvecmax, absmax, epsilon, estimated_error
      real betakm1, betak, dot2, betak_min, smax, step, alpha
      real convergence_rate, betaquot, s_dot_as, step1, cfactor
      real fxa, r2dtuv, f3, atosp, f2, uext, vext, d1, d2
      real diag1, diag0, diag3, diag4, dubdt, dvbdt

      include "size.h"

      integer iperm(maxipp), jperm(maxipp), iofs(mnisle)
      integer iterations, nippts(mnisle)

      logical imask(-mnisle:mnisle)

      real guess(imt,jmt), dpsi(imt,jmt), Zres(imt,jmt)
      real s(imt,jmt), res(imt,jmt), As(imt,jmt), forc(imt,jmt)
      real cf(imt,jmt,-1:1,-1:1), Z(imt,jmt)

!-----------------------------------------------------------------------
!     impose boundary conditions on guess
!     dpsi(0) = guess
!-----------------------------------------------------------------------

      call border(guess, bc_symm)

      do i=1,imt
        do j=1,jmt
          dpsi(i,j) = guess(i,j)
        enddo
      enddo

!-----------------------------------------------------------------------
!     make approximate inverse operator Z (always even symmetry)
!-----------------------------------------------------------------------

      call make_inv (cf, Z,
     &               imask, iperm, jperm, iofs, nisle, nippts)
      call border(Z, 't even')

!-----------------------------------------------------------------------
!     res(0)  = forc - A * dpsi(0)
!     impose cyclic and/or symmetry conditions on res(i,j)
!-----------------------------------------------------------------------

      if (npt .eq. 5) then
        call op5_vec(cf, dpsi, res)
      else
        call op9_vec(cf, dpsi, res)
      endif
      do i=2,imt-1
        do j=2,jmt-1
          res(i,j) = forc(i,j) - res(i,j)
        enddo
      enddo

      call border(res, bc_symm)

!-----------------------------------------------------------------------
!     Zres(k-1) = Z * res(k-1)
!     see if guess is a solution, bail out to avoid division by zero
!-----------------------------------------------------------------------

      k = 0
      call inv_op(Z, res, Zres,
     &            imask, iperm, jperm, iofs, nisle, nippts)

!     set borders of Zres using cyclic/symmetry, if defined.

      call border(Zres, bc_symm)

      Zresmax = absmax(Zres)

      diverging = .false.

!     Assume convergence rate of 0.99 to extrapolate error

      if (100.0 * Zresmax .lt. epsilon) then
        estimated_error = 100.0 * Zresmax
        goto 101
      endif

!-----------------------------------------------------------------------
!     beta(0) = 1
!     s(0)    = zerovector()
!-----------------------------------------------------------------------

      betakm1 = 1.0
      call zero_vec(s)

!-----------------------------------------------------------------------
!     begin iteration loop
!-----------------------------------------------------------------------

      do k = 1,max_iterations

!-----------------------------------------------------------------------
!       Zres(k-1) = Z * res(k-1)
!-----------------------------------------------------------------------

        call inv_op(Z, res, Zres,
     &              imask, iperm, jperm, iofs, nisle, nippts)

!       set borders of Zres using cyclic/symmetry, if defined.

        call border(Zres, bc_symm)

!-----------------------------------------------------------------------
!       beta(k)   = res(k-1) * Zres(k-1)
!-----------------------------------------------------------------------

        betak = dot2(Zres, res)
        if (k .eq. 1) then
          betak_min = abs(betak)
        elseif (k .gt. 2) then
          betak_min = min(betak_min, abs(betak))
          if (abs(betak) .gt. 100.0*betak_min) then
            write (*,'(/2(a/))')
     &      'WARNING: conjugate gradient solver terminated because'
     &,     '         correction steps are diverging.'
            write (*,'(/7(a/))')
     &      'PROBABLE CAUSES:'
     &,     '         1. convergence criterion is too tight...'
     &,     '            roundoff error prevents convergence'
     &,     '     or  2. the solution is beginning to blow up...'
     &,     '            if so, it is extremely unlikely that usable'
     &,     '            results can be obtained in subsequent time'
     &,     '            steps.'
            write (*,'(/3(a/))')
     &      'ERROR:   It is assumed that the solution is blowing up.'
     &,     '         It is extremely unlikely that usable results can'
     &,     '         be obtained in subsequent time steps.'
            if (variable .ne. 'surfpres') then
              stop '==>congrad'
            endif
            diverging = .true.
            smax = absmax(s)
            step = abs(alpha) * smax
            estimated_error=step*convergence_rate/(1.0-convergence_rate)
            go to 101
          endif
        endif

!-----------------------------------------------------------------------
!       s(k)      = Zres(k-1) + (beta(k)/beta(k-1)) * s(k-1)
!-----------------------------------------------------------------------

        betaquot = betak/betakm1
        do i=1,imt
          do j=1,jmt
            s(i,j) = Zres(i,j) + betaquot * s(i,j)
          enddo
        enddo

!-----------------------------------------------------------------------
!       As(k)     = A * s(k)
!-----------------------------------------------------------------------

        if (npt .eq. 5) then
          call op5_vec(cf, s, As)
        else
          call op9_vec(cf, s, As)
        endif

        call border(As, bc_symm)

!-----------------------------------------------------------------------
!       If s=0 then the division for alpha(k) gives a float exception.
!       Assume convergence rate of 0.99 to extrapolate error.
!       Also assume alpha(k) ~ 1.
!-----------------------------------------------------------------------

        s_dot_As = dot2(s, As)
        if (abs(s_dot_As) .lt. abs(betak)*1.e-10) then
          smax = absmax(s)
          estimated_error = 100.0 * smax
          goto 101
        endif

!-----------------------------------------------------------------------
!       alpha(k)  = beta(k) / (s(k) * As(k))
!-----------------------------------------------------------------------

        alpha = betak / s_dot_As

!-----------------------------------------------------------------------
!       update values:
!       dpsi(k)   = dpsi(k-1) + alpha(k) * s(k)
!       res(k)    = res(k-1) - alpha(k) * As(k)
!-----------------------------------------------------------------------

        do i=1,imt
          do j=1,jmt
            dpsi (i,j) = dpsi(i,j) + alpha * s(i,j)
            res  (i,j) = res (i,j) - alpha * As(i,j)
          enddo
        enddo
        call avg_dist (res,
     &          imask, iperm, jperm, iofs, nisle, nippts)

        call border(res, bc_symm)

        smax = absmax(s)

!-----------------------------------------------------------------------
!       test for convergence
!       if (estimated_error) < epsilon) exit
!-----------------------------------------------------------------------

        step = abs(alpha) * smax
        if (k .eq. 1) then
          step1 = step
          estimated_error = step
          if (step .lt. epsilon) goto 101
        elseif (step .lt. epsilon) then
          cfactor = log(step/step1)
          convergence_rate = exp(cfactor/(k-1))
          estimated_error = step*convergence_rate/(1.0-convergence_rate)
          if (estimated_error .lt. epsilon) goto 101
        endif

        betakm1 = betak

      enddo

!-----------------------------------------------------------------------
!     end of iteration loop
!-----------------------------------------------------------------------

  101 continue
      if (k .gt. max_iterations) then
        cfactor = log(step/step1)
        convergence_rate = exp(cfactor/(k-1))
        estimated_error = step*convergence_rate/(1.0-convergence_rate)
        converged = .false.
      else
        if (diverging) then
          converged = .false.
        else
          converged = .true.
        endif
      endif

      iterations = k

!-----------------------------------------------------------------------
!     return the last increment of dpsi in the argument res
!-----------------------------------------------------------------------

      if (iterations .eq. 0) then
        do i=1,imt
          do j=1,jmt
            res(i,j) = Zres(i,j)
          enddo
        enddo
      else
        do i=1,imt
          do j=1,jmt
            res(i,j) = alpha * s(i,j)
          enddo
        enddo
      endif

      return
      end

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

!     M A T R I X   M O D U L E   F O R   C O N G R A D

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

      subroutine zero_vec (v)

      implicit none

      integer i, j

      include "size.h"

      real v(imt,jmt)

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

      return
      end

      subroutine add_vec (v,w,vpw)

      implicit none

      integer i, j

      include "size.h"

      real v(imt,jmt), w(imt,jmt), vpw(imt,jmt)

      do i=1,imt
        do j=1,jmt
          vpw(i,j) = v(i,j) + w(i,j)
        enddo
      enddo

      return
      end

      subroutine sub_vec (v,w,vmw)

      implicit none

      integer i, j

      include "size.h"

      real v(imt,jmt), w(imt,jmt), vmw(imt,jmt)

      do i=1,imt
        do j=1,jmt
          vmw(i,j) = v(i,j) - w(i,j)
        enddo
      enddo

      return
      end

      subroutine mult_vec(v,w,vtw)

      implicit none

      integer i, j

      include "size.h"

      real v(imt,jmt), w(imt,jmt), vtw(imt,jmt)

      do i=1,imt
        do j=1,jmt
          vtw(i,j) = v(i,j) * w(i,j)
        enddo
      enddo

      return
      end

      subroutine div_vec(v,w,vdw)

      implicit none

      integer i, j

      include "size.h"

      real v(imt,jmt), w(imt,jmt), vdw(imt,jmt)

      do i=1,imt
        do j=1,jmt
          if (w(i,j) .ne. 0) then
            vdw(i,j) = v(i,j) / w(i,j)
          else
            vdw(i,j) = 0.0
          endif
        enddo
      enddo

      return
      end

      subroutine scalar_vec (scalar,w,sw)

      implicit none

      integer i, j

      include "size.h"

      real w(imt,jmt), sw(imt,jmt), scalar

      do i=1,imt
        do j=1,jmt
          sw(i,j) = scalar * w(i,j)
        enddo
      enddo

      return
      end

      subroutine neg_vec (v)

      implicit none

      integer i, j

      include "size.h"

      real v(imt,jmt)

      do i=1,imt
        do j=1,jmt
          v(i,j) = -v(i,j)
        enddo
      enddo

      return
      end

      function dot2 (dp_vec, res_vec)

!     this dot product produces the correct answers because for
!     ocean perimeter subscripts, ij=isle, the value on a
!     type(dpsi_type) vector, dp_vec(isle)=dp_vec(i,j), i.e., the true
!     value is replicated, and for a type(res_type) vector,
!     res_vec(isle) = sum (res_vec(i,j)), i.e., the true value is the
!     accumulation of the distributed values.

      implicit none

      integer j, i

      real dot2

      include "size.h"

      real dp_vec(imt,jmt), res_vec(imt,jmt), rowsum (jmt)

      do j=2,jmt-1
        rowsum(j) = 0.0
        do i=2,imt-1
          rowsum(j) = rowsum(j) + dp_vec(i,j) * res_vec(i,j)
        enddo
      enddo

      dot2 = 0.0
      do j=2,jmt-1
        dot2 = dot2 + rowsum(j)
      enddo

      return
      end

      subroutine op5_vec(cf, dpsi, res)

!                       res = A * dpsi

!     this subroutine does not collect the terms of the true value
!     of res(isle) = sum (res(i,j)).  the contributions to the sum
!     remain distributed among the T cells (i,j) that form the
!     ocean perimeter of land_mass(isle).

!     at present, borders are not computed [i=1 or imt] [j=1 or jmt]

      implicit none

      integer i, j

      include "size.h"

      real cf(imt,jmt,-1:1,-1:1), dpsi(imt,jmt), res(imt,jmt)

      do j=2,jmt-1
        do i=2,imt-1
          res(i,j) = cf(i,j, 0, 0) * dpsi(i,j)   +
     &               cf(i,j, 0, 1) * dpsi(i,j+1) +
     &               cf(i,j, 0,-1) * dpsi(i,j-1) +
     &               cf(i,j, 1, 0) * dpsi(i+1,j) +
     &               cf(i,j,-1, 0) * dpsi(i-1,j)
        enddo
      enddo

      return
      end

      subroutine op9_vec(cf, dpsi, res)

!                       res = A * dpsi

!     this subroutine does not collect the terms of the true value
!     of res(isle) = sum (res(i,j)).  the contributions to the sum
!     remain distributed among the T cells (i,j) that form the
!     ocean perimeter of land_mass(isle).

!     at present, borders are not computed [i=1 or imt] [j=1 or jmt]

      implicit none

      integer i, j

      include "size.h"

      real cf(imt,jmt,-1:1,-1:1), dpsi(imt,jmt), res(imt,jmt)

      do j=2,jmt-1
        do i=2,imt-1
          res(i,j) = cf(i,j, 0, 0) * dpsi(i  ,j  ) +
     &               cf(i,j, 0, 1) * dpsi(i  ,j+1) +
     &               cf(i,j, 0,-1) * dpsi(i  ,j-1) +
     &               cf(i,j, 1, 0) * dpsi(i+1,j  ) +
     &               cf(i,j,-1, 0) * dpsi(i-1,j  ) +
     &               cf(i,j, 1, 1) * dpsi(i+1,j+1) +
     &               cf(i,j,-1, 1) * dpsi(i-1,j+1) +
     &               cf(i,j, 1,-1) * dpsi(i+1,j-1) +
     &               cf(i,j,-1,-1) * dpsi(i-1,j-1)
        enddo
      enddo

      return
      end

      subroutine subset (a, b, nerror)

!     verifies that the set of subscripts for which a(i,j) .ne. 0.0
!     is a subset of the set of subscripts for which b(i,j) .ne. 0.0

      implicit none

      integer i, j, nerror

      include "size.h"

      real a(imt,jmt), b(imt,jmt)

      nerror = 0
      do i=2,imt-1
        do j=2,jmt-1
          if (a(i,j) .ne. 0.0 .and. b(i,j) .eq. 0.0) then
            nerror = nerror + 1
            print '(a,i3,a,i3,a,a)', '(',i,',',j,')'
     &                           ,' forcing is reset to zero'
!           set forcing (i.e., a(i,j)) to zero
            a(i,j) = 0.0
          endif
        enddo
      enddo

      return
      end

      subroutine inv_op(Z, res, Zres,
     &              imask, iperm, jperm, iofs, nisle, nippts)

!     apply and approximate inverse Z or the operator A

!     res is type(res_type), i.e., perimeter values res(isle)
!         are the sum of the distributed contributions res(i,j)
!     Zres is type(dpsi_type), i.e., perimeter values Zres(isle)
!         must be replicated at each perimeter point Zres(i,j)

!     borders  of Zres [i=1 or imt] [j=1 or jmt] must be defined
!     and must satisfy cyclic and/or symmetry, if defined.

!     currently, Z is diagonal:  Z(ij) = 1/A(ij)
!     and is stored in type(dpsi_type) format, i.e., Z(isle) is
!     replicated and stored in each Z(i,j) in the perimeter of
!     land_mass(isle).

      implicit none

      integer i, j, nisle

      include "size.h"

      logical imask(-mnisle:mnisle)

      integer iperm(maxipp), jperm(maxipp), iofs(mnisle), nippts(mnisle)

      real Z(imt,jmt), res(imt,jmt), Zres(imt,jmt)

      do i=1,imt
        do j=1,jmt
          Zres(i,j) = Z(i,j) * res(i,j)
        enddo
      enddo

!     sum contributions to Zres(isle)
!     distribute Zres(isle) to all perimeter points

      call sum_dist (Zres,
     &        imask, iperm, jperm, iofs, nisle, nippts)

      return
      end

      function absvecmax(res, imax, jmax)

      implicit none

      integer i, j, imax, jmax

      real absvecmax

      include "size.h"

      real res(imt,jmt)

      absvecmax = 0.0
      do i=2,imt-1
        do j=2,jmt-1
          if (abs(res(i,j)) .gt. absvecmax) then
            absvecmax = abs(res(i,j))
            imax = i
            jmax = j
          endif
        enddo
      enddo
      return
      end

      function absmax (f)

      implicit none

      integer i, j

      real amax, absmax

      include "size.h"

      real f(imt,jmt)

      amax = 0.0
      do i=1,imt
        do j=1,jmt
          amax = max(amax, abs(f(i,j)))
        enddo
      enddo
      absmax = amax
      return
      end

      function absmin (f)

      implicit none

      integer i, j

      real amin, absmin

      include "size.h"

      real f(imt,jmt)

      amin = 1.0e37
      do i=1,imt
        do j=1,jmt
          if (f(i,j) .ne. 0 .and. abs(f(i,j)) .lt. amin) then
            amin = abs(f(i,j))
          endif
        enddo
      enddo
      absmin = amin
      return
      end

      subroutine make_inv (cf, Z,
     &              imask, iperm, jperm, iofs, nisle, nippts)

!     construct an approximate inverse Z to A

!     Z will be diagonal:  Z(ij) = 1/A(ij)
!     and values for ocean perimeter entries Z(isle) will be replicated
!     at all T cells Z(i,j) in the ocean perimeter of land_mass(isle).

!     T cells (i,j) for which there is no diagonal coefficient
!     i.e., A(ij)=A(i,j)=0, are masked off by assigning Z(i,j)=0.
!     there are effectively no equations and no variables dpsi(i,j)
!     at these points.

      implicit none

      integer i, j, nisle, isle, n

      include "size.h"

      integer iperm(maxipp),jperm(maxipp)
      integer iofs(mnisle), nippts(mnisle)

      logical imask(-mnisle:mnisle)

      real cf(imt,jmt,-1:1,-1:1), Z(imt,jmt)

!     copy diagonal coefficients of A to Z

      do i=2,imt-1
        do j=2,jmt-1
          Z(i,j) = cf(i,j,0,0)
        enddo
      enddo

!     for each land_mass(isle),
!     sum the contributions to cfdiag(isle)=A(isle,isle)
!     now stored in Z(i,j) at ocean perimeter T cells and replicate
!     the sum in all Z(i,j) for which (i,j) is in ocean perimeter
!     of land_mass(isle).

      call sum_dist (Z,
     &        imask, iperm, jperm, iofs, nisle, nippts)

!     now invert Z

      do i=2,imt-1
        do j=2,jmt-1
          if (Z(i,j) .ne. 0.0) then
            Z(i,j) = 1/Z(i,j)
          else
            Z(i,j) = 0.0
          endif
        enddo
      enddo

!     make inverse zero on island perimeters that are not integrated

      do isle=1,nisle
        if (.not. imask(isle)) then
          do n=1,nippts(isle)
            i = iperm(iofs(isle)+n)
            j = jperm(iofs(isle)+n)
            Z(i,j) = 0.0
          enddo
        endif
      enddo

      return
      end

      subroutine sum_dist (Zres,
     &              imask, iperm, jperm, iofs, nisle, nippts)

!     sum contributions to Zres(isle)
!     distribute Zres(isle) to all perimeter points

!     this subroutine converts a type(res_type) vector with
!     distributed contributions to perimeter values
!        Zres(isle) = sum (Zres(i,j))
!     into a type (dpsi_type) vector with replicated values
!     for land_mass perimeters
!        Zres(isle) = Zres(i,j)
!     for all (i,j) in the ocean perimeter of land_mass(isle).

      implicit none

      integer i, j, isle, nisle, n

      include "size.h"

      integer iperm(maxipp),jperm(maxipp)
      integer iofs(mnisle), nippts(mnisle)

      logical imask(-mnisle:mnisle)

      real Zres(imt,jmt), Zresisle(mnisle)

!     sum contributions to Zres(isle)

      do isle=1,nisle
        if (imask(isle)) then
          Zresisle(isle) = 0.0
          do n=1,nippts(isle)
            i = iperm(iofs(isle)+n)
            j = jperm(iofs(isle)+n)
            Zresisle(isle) = Zresisle(isle) + Zres(i,j)
          enddo
        endif
      enddo

!     distribute Zres(isle) to all perimeter points

      do isle=1,nisle
        if (imask(isle)) then
          do n=1,nippts(isle)
            i = iperm(iofs(isle)+n)
            j = jperm(iofs(isle)+n)
            Zres(i,j) = Zresisle(isle)
          enddo
        endif
      enddo

      return
      end

      subroutine avg_dist (Zres,
     &              imask, iperm, jperm, iofs, nisle, nippts)

!     avg contributions to Zres(isle)
!     distribute Zres(isle) to all perimeter points

!     this subroutine converts a type(res_type) vector with
!     distributed contributions to perimeter values
!        Zres(isle) = avg (Zres(i,j))
!     into a type (dpsi_type) vector with replicated values
!     for land_mass perimeters
!        Zres(isle) = Zres(i,j)
!     for all (i,j) in the ocean perimeter of land_mass(isle).

      implicit none

      integer i, j, isle, nisle, n

      include "size.h"

      integer iperm(maxipp),jperm(maxipp)
      integer iofs (mnisle), nippts(mnisle)

      logical imask(-mnisle:mnisle)

      real Zres(imt,jmt), Zresisle(mnisle)

!     avg contributions to Zres(isle)

      do isle=1,nisle
        if (imask(isle)) then
!         print *,'isle=',isle,' nisle=',nisle
          Zresisle(isle) = 0.0
          do n=1,nippts(isle)
            i = iperm(iofs(isle)+n)
            j = jperm(iofs(isle)+n)
            Zresisle(isle) = Zresisle(isle) + Zres(i,j)
          enddo
        endif
      enddo

!     distribute Zres(isle) to all perimeter points

      do isle=1,nisle
        if (imask(isle)) then
          do n=1,nippts(isle)
            i = iperm(iofs(isle)+n)
            j = jperm(iofs(isle)+n)
            Zres(i,j) = Zresisle(isle)/nippts(isle)
          enddo
        endif
      enddo

      return
      end