#include "STREAMICE_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP

      SUBROUTINE STREAMICE_TRIDIAG_SOLVE( 
     U                               cg_Uin,     ! x-velocities
     U                               cg_Vin,     ! x-velocities
     U                               cg_Bu,      ! force in x dir
     I                               A_uu,       ! section of matrix that multiplies u and projects on u
     I                               umask,
     I                               myThid )
C     /============================================================\
C     | SUBROUTINE                                                 |   
C     | o                                                          |
C     |============================================================|
C     |                                                            |
C     \============================================================/
      IMPLICIT NONE

#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "STREAMICE.h"
#include "STREAMICE_CG.h"

      _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
      _RS umask (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      INTEGER myThid

      INTEGER iMin,iMax,i,j
      _RL aMat(1:Nx)
      _RL bMat(1:Nx)
      _RL cMat(1:Nx)
      _RL yMat(1:Nx)
      _RL bet(1:Nx)
      _RL tmpvar
      INTEGER errCode
 
      
      IF (nPx.gt.1 .or. nSx.gt.1) THEN
       STOP 'must be serial for tridiag solve'
      ENDIF

      errCode = 0

      imax = 0
      iMin = 2
      do i=imin,Nx
       if (umask(i,1,1,1).eq.1.0) THEN

        aMat(i)=0.0
        bmat(i)=0.0
        cmat(i)=0.0
        ymat(i)=0.0
        do j=-1,1
        aMat(i) = amat(i)+A_uu(i,1,1,1,-1,j)
        bMat(i) = bmat(i)+A_uu(i,1,1,1,0,j)
        cMat(i) = cmat(i)+A_uu(i,1,1,1,1,j)
        enddo
        yMat(i) = ymat(i)+cg_Bu(i,1,1,1)
       else
        iMax = i-1
        exit
       endif
      enddo
      
      IF(imax.eq.0) THEN
       imax=Nx
      ENDIF
 

      IF ( bMat(imin).NE.0. _d 0 ) THEN 
       bet(imin) = 1. _d 0 / bMat(imin)
      ELSE
       bet(imin) = 0. _d 0
       errCode = 1
      ENDIF
       
      DO i=imin+1,imax
       tmpvar = bmat(i) - amat(i)*cmat(i-1)*bet(i-1)
       IF ( tmpvar .NE. 0. _d 0 ) THEN
        bet(i) = 1. _d 0 / tmpvar
       ELSE
        bet(i) = 0. _d 0
        errCode = 1
       ENDIF
      ENDDO


      ymat(imin) = ymat(imin)*bet(imin)

      DO i=imin+1,imax
       ymat(i) = ( ymat(i) 
     &            - amat(i)*ymat(i-1)
     &            )*bet(i) 
      ENDDO


      DO i=imax-1,imin,-1
          ymat(i) = ymat(i)
     &     - cmat(i)*bet(i)*ymat(i+1)
      ENDDO

      DO j=1,sNy
       DO i=imin,imax
        cg_Uin (i,j,1,1) = ymat(i)
       ENDDO
      ENDDO

      DO j=1,sNy
       DO i=1,sNx
        cg_Vin (i,j,1,1) = 0. _d 0
       ENDDO
      ENDDO

      print *, "ERRORCODE", errcode

      RETURN
      END