#include "RADTRANS_OPTIONS.h"

CBOP
C     !ROUTINE: RADTRANS_SOLVE_TRIDIAG

C     !INTERFACE:
      SUBROUTINE RADTRANS_SOLVE_TRIDIAG(
     U                     a3d, b3d, c3d,
     U                     y3d,
     I                     n, myThid )

C     !DESCRIPTION:
C     Solve a tri-diagonal system A*X=Y (dimension n) by Gaussian
C     elimination with pivoting
C     - input arrays are overwritten
C     - result is returned in y3d

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

C     !INPUT PARAMETERS:
      INTEGER n,myThid

C     !INPUT/OUTPUT PARAMETERS:
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)
CEOP

#ifdef ALLOW_RADTRANS

C     !LOCAL VARIABLES:
C     d3d(1:Nr-2) :: second upper diagnonal, generated by pivoting
      INTEGER k
      _RL tmpc, tmpy
      _RL recVar
      _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

#endif /*ALLOW_RADTRANS */

      RETURN
      END