C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_tridiag_solve.F,v 1.2 2013/06/21 20:49:51 jmc Exp $
C $Name:  $

#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,k
      _RL aMat(1:Nx)
      _RL bMat(1:Nx)
      _RL cMat(1:Nx)
      _RL yMat(1:Nx)
      _RL bet(1:Nx)
      _RL tmpvar
      INTEGER errCode


!      CALL WRITE_FLD_XY_RL ("taud_tri","",cg_Bu,0,mythid)
!      CALL WRITE_FLD_XY_RL ("A_m1m1","",A_uu(:,:,:,:,-1,-1),0,mythid)
!      CALL WRITE_FLD_XY_RL ("A_m1_0","",A_uu(:,:,:,:,-1,0),0,mythid)
!      CALL WRITE_FLD_XY_RL ("A_m1p1","",A_uu(:,:,:,:,-1,1),0,mythid)
!      CALL WRITE_FLD_XY_RL ("A_0_m1","",A_uu(:,:,:,:,0,-1),0,mythid)
!      CALL WRITE_FLD_XY_RL ("A_0_0","",A_uu(:,:,:,:,0,0),0,mythid)
!      CALL WRITE_FLD_XY_RL ("A_0_p1","",A_uu(:,:,:,:,0,1),0,mythid)
!      CALL WRITE_FLD_XY_RL ("A_p1m1","",A_uu(:,:,:,:,1,-1),0,mythid)
!      CALL WRITE_FLD_XY_RL ("A_p1_0","",A_uu(:,:,:,:,1,0),0,mythid)
!      CALL WRITE_FLD_XY_RL ("A_p1p1","",A_uu(:,:,:,:,1,1),0,mythid)



      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
        do k=1,3
        aMat(i) = amat(i)+A_uu(i,k,1,1,-1,j)
        bMat(i) = bmat(i)+A_uu(i,k,1,1,0,j)
        cMat(i) = cmat(i)+A_uu(i,k,1,1,1,j)
        enddo
        yMat(i) = ymat(i)+cg_Bu(i,j+2,1,1)
        enddo
       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