C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/darwin2/pkg/monod/solve_tridiagonal_pivot.F,v 1.1 2012/08/23 21:53:00 jahn Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: SOLVE_TRIDIAGONAL_PIVOT
C     !INTERFACE:
      SUBROUTINE SOLVE_TRIDIAGONAL_PIVOT(
     U                     a3d, b3d, c3d,
     U                     y3d,
     I                     n, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R SOLVE_TRIDIAGONAL_PIVOT
C     | o Solve a tri-diagonal system A*X=Y (dimension n)
C     |   by Gaussian elimination with pivoting
C     *==========================================================*
C     | o input arrays are overwritten
C     | o result is returned in y3d
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global data ==
#include "SIZE.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
C     a3d(2:n)   :: matrix lower diagnonal
C     b3d(1:n)   :: matrix main  diagnonal
C     c3d(1:n-1) :: matrix upper diagnonal
C     y3d(1:n)   :: Input = Y vector ; Output = X = solution of A*X=Y
      _RL a3d(2*Nr)
      _RL b3d(2*Nr)
      _RL c3d(2*Nr)
      _RL y3d(2*Nr)
      INTEGER n,myThid
CEOP

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER k
      _RL tmpc, tmpy
      _RL recVar
C     d3d(1:Nr-2) :: second upper diagnonal, generated by pivoting
      _RL d3d(2*Nr)

      c3d(n) = 0.

c eliminate lower diagonal
c     only c3d, d3d and y3d are modified and will be used later
      DO k=1,n-1
        IF(ABS(b3d(k)).GE.ABS(a3d(k+1)))THEN
c         scale current row, subtract from next
          recVar = 1. / b3d(k)
          c3d(k) = c3d(k)*recVar
          y3d(k) = y3d(k)*recVar
          b3d(k+1) = b3d(k+1) - a3d(k+1)*c3d(k)
          y3d(k+1) = y3d(k+1) - a3d(k+1)*y3d(k)
          d3d(k) = 0.
        ELSE
c         swap rows, then scale current and subtract from next
          recVar = 1. / a3d(k+1)
          tmpc = c3d(k)
          tmpy = y3d(k)
          c3d(k) = b3d(k+1)*recVar
          d3d(k) = c3d(k+1)*recVar
          y3d(k) = y3d(k+1)*recVar
          b3d(k+1) = tmpc - b3d(k)*c3d(k)
          c3d(k+1) =      - b3d(k)*d3d(k)
          y3d(k+1) = tmpy - b3d(k)*y3d(k)
        ENDIF
      ENDDO
      recVar = 1. / b3d(n)
      y3d(n) = y3d(n)*recVar

c backsubstitute
c     y3d(n) is already good
      y3d(n-1) = y3d(n-1) - c3d(n-1)*y3d(n)
      DO k=n-2,1,-1
        y3d(k) = y3d(k) - c3d(k)*y3d(k+1) - d3d(k)*y3d(k+2)
      ENDDO

      RETURN
      END