
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppDecomp( JVS, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse LU factorization
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE KPP_ROOT_Parameters
  USE KPP_ROOT_JacobianSP

      INTEGER  :: IER
      KPP_REAL :: JVS(LU_NONZERO), W(NVAR), a
      INTEGER  :: k, kk, j, jj

      a = 0. ! mz_rs_20050606
      IER = 0
      DO k=1,NVAR
        ! mz_rs_20050606: don't check if real value == 0
        ! IF ( JVS( LU_DIAG(k) ) .EQ. 0. ) THEN
        IF ( ABS(JVS(LU_DIAG(k))) < TINY(a) ) THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              W( LU_ICOL(kk) ) = JVS(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            a = -W(j) / JVS( LU_DIAG(j) )
            W(j) = -a
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVS(kk) = W( LU_ICOL(kk) )
         END DO
      END DO
      
END SUBROUTINE KppDecomp


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppDecompCmplx( JVS, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse LU factorization, complex
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE KPP_ROOT_Parameters
  USE KPP_ROOT_JacobianSP

      INTEGER        :: IER
      DOUBLE COMPLEX :: JVS(LU_NONZERO), W(NVAR), a
      KPP_REAL  :: b = 0.0
      INTEGER        :: k, kk, j, jj

      IER = 0
      DO k=1,NVAR
        IF ( ABS(JVS(LU_DIAG(k))) < TINY(b) ) THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              W( LU_ICOL(kk) ) = JVS(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            a = -W(j) / JVS( LU_DIAG(j) )
            W(j) = -a
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVS(kk) = W( LU_ICOL(kk) )
         END DO
      END DO
      
END SUBROUTINE KppDecompCmplx


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppDecompCmplxR( JVSR, JVSI, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Sparse LU factorization, complex
!   (Real and Imaginary parts are used instead of complex data type)     
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE KPP_ROOT_Parameters
  USE KPP_ROOT_JacobianSP

      INTEGER       :: IER
      KPP_REAL :: JVSR(LU_NONZERO), JVSI(LU_NONZERO) 
      KPP_REAL :: WR(NVAR), WI(NVAR), ar, ai, den
      INTEGER       :: k, kk, j, jj

      IER = 0
      ar  = 0.0
      DO k=1,NVAR
        IF (  ( ABS(JVSR(LU_DIAG(k))) < TINY(ar) ) .AND. &
              ( ABS(JVSI(LU_DIAG(k))) < TINY(ar) ) )  THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              WR( LU_ICOL(kk) ) = JVSR(kk)
              WI( LU_ICOL(kk) ) = JVSI(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            den = JVSR(LU_DIAG(j))**2 + JVSI(LU_DIAG(j))**2
            ar = -(WR(j)*JVSR(LU_DIAG(j)) + WI(j)*JVSI(LU_DIAG(j)))/den
            ai = -(WI(j)*JVSR(LU_DIAG(j)) - WR(j)*JVSI(LU_DIAG(j)))/den
            WR(j) = -ar
            WI(j) = -ai
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               WR( LU_ICOL(jj) ) = WR( LU_ICOL(jj) ) + ar*JVSR(jj) - ai*JVSI(jj)
               WI( LU_ICOL(jj) ) = WI( LU_ICOL(jj) ) + ar*JVSI(jj) + ai*JVSR(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVSR(kk) = WR( LU_ICOL(kk) )
            JVSI(kk) = WI( LU_ICOL(kk) )
         END DO
      END DO

END SUBROUTINE KppDecompCmplxR


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveIndirect( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse solve subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE KPP_ROOT_Parameters
  USE KPP_ROOT_JacobianSP

      INTEGER  :: i, j
      KPP_REAL :: JVS(LU_NONZERO), X(NVAR), sum

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1 
             X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
         END DO  
      END DO

      DO i=NVAR,1,-1
        sum = X(i);
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
          sum = sum - JVS(j)*X(LU_ICOL(j));
        END DO
        X(i) = sum/JVS(LU_DIAG(i));
      END DO
      
END SUBROUTINE KppSolveIndirect


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveTRIndirect( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Complex sparse solve transpose subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE KPP_ROOT_Parameters
  USE KPP_ROOT_JacobianSP

      INTEGER       :: i, j
      KPP_REAL :: JVS(LU_NONZERO), X(NVAR)

      DO i=1,NVAR
        X(i) = X(i)/JVS(LU_DIAG(i))
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO

      DO i=NVAR, 1, -1
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_CROW(i),LU_DIAG(i)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO
      
END SUBROUTINE KppSolveTRIndirect


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveCmplx( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Complex sparse solve subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE KPP_ROOT_Parameters
  USE KPP_ROOT_JacobianSP

      INTEGER        :: i, j
      DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR), sum

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1 
             X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
         END DO  
      END DO

      DO i=NVAR,1,-1
        sum = X(i);
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
          sum = sum - JVS(j)*X(LU_ICOL(j));
        END DO
        X(i) = sum/JVS(LU_DIAG(i));
      END DO
      
END SUBROUTINE KppSolveCmplx

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveCmplxR( JVSR, JVSI, XR, XI )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   Complex sparse solve subroutine using indirect addressing
!   (Real and Imaginary parts are used instead of complex data type)     
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE KPP_ROOT_Parameters
  USE KPP_ROOT_JacobianSP

      INTEGER       ::  i, j
      KPP_REAL ::  JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), sumr, sumi, den

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1 
             XR(i) = XR(i) - (JVSR(j)*XR(LU_ICOL(j)) - JVSI(j)*XI(LU_ICOL(j)))
             XI(i) = XI(i) - (JVSR(j)*XI(LU_ICOL(j)) + JVSI(j)*XR(LU_ICOL(j)))
         END DO  
      END DO

      DO i=NVAR,1,-1
        sumr = XR(i); sumi = XI(i)
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
            sumr = sumr - (JVSR(j)*XR(LU_ICOL(j)) - JVSI(j)*XI(LU_ICOL(j)))
            sumi = sumi - (JVSR(j)*XI(LU_ICOL(j)) + JVSI(j)*XR(LU_ICOL(j)))
        END DO
        den   = JVSR(LU_DIAG(i))**2 + JVSI(LU_DIAG(i))**2
        XR(i) = (sumr*JVSR(LU_DIAG(i)) + sumi*JVSI(LU_DIAG(i)))/den
        XI(i) = (sumi*JVSR(LU_DIAG(i)) - sumr*JVSI(LU_DIAG(i)))/den
      END DO
      
END SUBROUTINE KppSolveCmplxR


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveTRCmplx( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Complex sparse solve transpose subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE KPP_ROOT_Parameters
  USE KPP_ROOT_JacobianSP

      INTEGER        :: i, j
      DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR)

      DO i=1,NVAR
        X(i) = X(i)/JVS(LU_DIAG(i))
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO

      DO i=NVAR, 1, -1
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_CROW(i),LU_DIAG(i)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO
      
END SUBROUTINE KppSolveTRCmplx


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveTRCmplxR( JVSR, JVSI, XR, XI )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   Complex sparse solve transpose subroutine using indirect addressing
!   (Real and Imaginary parts are used instead of complex data type)     
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE KPP_ROOT_Parameters
  USE KPP_ROOT_JacobianSP

      INTEGER       ::  i, j
      KPP_REAL ::  JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), den

      DO i=1,NVAR
        den   = JVSR(LU_DIAG(i))**2 + JVSI(LU_DIAG(i))**2
        XR(i) = (XR(i)*JVSR(LU_DIAG(i)) + XI(i)*JVSI(LU_DIAG(i)))/den
        XI(i) = (XI(i)*JVSR(LU_DIAG(i)) - XR(i)*JVSI(LU_DIAG(i)))/den
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
          XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)*XR(i) - JVSI(j)*XI(i))
          XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)*XR(i) + JVSR(j)*XI(i))
        END DO
      END DO

      DO i=NVAR, 1, -1
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_CROW(i),LU_DIAG(i)-1
          XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)*XR(i) - JVSI(j)*XI(i))
          XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)*XR(i) + JVSR(j)*XI(i))
        END DO
      END DO
      
END SUBROUTINE KppSolveTRCmplxR


!
! Next few commented subroutines perform sparse big linear algebra
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE KppDecompBig( JVS, IP, IER )
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!        Sparse LU factorization
!!        for the Runge Kutta (3n)x(3n) linear system
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  USE KPP_ROOT_Parameters
!  USE KPP_ROOT_JacobianSP
!
!      INTEGER  :: IP3(3), IER, IP(3,NVAR)
!      KPP_REAL :: JVS(3,3,LU_NONZERO), W(3,3,NVAR), a(3,3), E(3,3)
!      INTEGER  :: k, kk, j, jj
!
!      a = 0.0d0
!      IER = 0
!      DO k=1,NVAR
!        DO kk = LU_CROW(k), LU_CROW(k+1)-1
!              W( 1:3,1:3,LU_ICOL(kk) ) = JVS(1:3,1:3,kk)
!        END DO
!        DO kk = LU_CROW(k), LU_DIAG(k)-1
!            j = LU_ICOL(kk)
!            E(1:3,1:3) = JVS( 1:3,1:3,LU_DIAG(j) )
!            ! CALL DGETRF(3,3,E,3,IP3,IER) 
!            CALL FAC3(E,IP3,IER)
!            IF ( IER /= 0 )  RETURN
!            ! a = W(j) / JVS( LU_DIAG(j) )
!            a(1:3,1:3) = W( 1:3,1:3,j )
!            ! CALL DGETRS ('N',3,3,E,3,IP3,a,3,IER) 
!            CALL SOL3('N',E,IP3,a(1,1))
!            CALL SOL3('N',E,IP3,a(1,2))
!            CALL SOL3('N',E,IP3,a(1,3))
!            W(1:3,1:3,j) = a(1:3,1:3)
!            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
!               W( 1:3,1:3,LU_ICOL(jj) ) = W( 1:3,1:3,LU_ICOL(jj) ) &
!                        - MATMUL( a(1:3,1:3) , JVS(1:3,1:3,jj) )
!            END DO
!         END DO
!         DO kk = LU_CROW(k), LU_CROW(k+1)-1
!            JVS(1:3,1:3,kk) = W( 1:3,1:3,LU_ICOL(kk) )
!         END DO
!      END DO
!
!      DO k=1,NVAR
!         ! CALL WGEFA(JVS(1,1,LU_DIAG(k)),3,3,IP(1,k),IER)
!         ! CALL DGETRF(3,3,JVS(1,1,LU_DIAG(k)),3,IP(1,k),IER)
!         CALL FAC3(JVS(1,1,LU_DIAG(k)),IP(1,k),IER)
!         IF ( IER /= 0 )  RETURN
!      END DO 
!      
!END SUBROUTINE KppDecompBig
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE KppSolveBig( JVS, IP, X )
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!        Sparse solve subroutine using indirect addressing
!!        for the Runge Kutta (3n)x(3n) linear system
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  USE KPP_ROOT_Parameters
!  USE KPP_ROOT_JacobianSP
!
!      INTEGER  :: i, j, k, m, IP3(3), IP(3,NVAR), IER
!      KPP_REAL :: JVS(3,3,LU_NONZERO), X(3,NVAR), sum(3)
!
!      DO i=1,NVAR
!        DO j = LU_CROW(i), LU_DIAG(i)-1 
!          !X(1:3,i) = X(1:3,i) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j)));
!          DO k=1,3
!            DO m=1,3
!              X(k,i) = X(k,i) - JVS(k,m,j)*X(m,LU_ICOL(j))
!            END DO
!          END DO
!        END DO  
!      END DO
!
!      DO i=NVAR,1,-1
!        sum(1:3) = X(1:3,i);
!        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
!          !sum(1:3) = sum(1:3) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j)));
!          DO k=1,3
!            DO m=1,3
!              sum(k) = sum(k) - JVS(k,m,j)*X(m,LU_ICOL(j))
!            END DO
!          END DO
!        END DO
!        ! X(i) = sum/JVS(LU_DIAG(i));
!        ! CALL DGETRS ('N',3,1,JVS(1:3,1:3,LU_DIAG(i)),3,IP(1,i),sum,3,0) 
!        ! CALL WGESL('N',JVS(1,1,LU_DIAG(i)),3,3,IP(1,i),sum)
!        CALL SOL3('N',JVS(1,1,LU_DIAG(i)),IP(1,i),sum)
!        X(1:3,i) = sum(1:3)
!      END DO
!      
!END SUBROUTINE KppSolveBig
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE KppSolveBigTR( JVS, IP, X )
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!        Big sparse transpose solve using indirect addressing
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  USE KPP_ROOT_Parameters
!  USE KPP_ROOT_JacobianSP
!
!      INTEGER       :: i, j, k, m, IP(3,NVAR)
!      KPP_REAL :: JVS(3,3,LU_NONZERO), X(3,NVAR)
!
!      DO i=1,NVAR
!        ! X(i) = X(i)/JVS(LU_DIAG(i))
!        CALL SOL3('T',JVS(1,1,LU_DIAG(i)),IP(1,i),X(1,i))
!        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
!         !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) &
!          !    - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) )
!          DO k=1,3
!            DO m=1,3
!              X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i)
!            END DO
!          END DO
!       END DO
!      END DO
!
!      DO i=NVAR, 1, -1
!        DO j=LU_CROW(i),LU_DIAG(i)-1
!         !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) &
!          !   - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) )
!          DO k=1,3
!            DO m=1,3
!              X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i)
!            END DO
!          END DO
!       END DO
!      END DO
!      
!END SUBROUTINE KppSolveBigTR
!
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE FAC3(A,IPVT,INFO)
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!     FAC3 FACTORS THE MATRIX A (3,3) BY
!!           GAUSS ELIMINATION WITH PARTIAL PIVOTING
!!     LINPACK - LIKE 
!!
!!     Remove comments to perform pivoting
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!
!      KPP_REAL :: A(3,3)
!      INTEGER       :: IPVT(3),INFO
!!      INTEGER       :: L
!!      KPP_REAL :: t, dmax, da, TMP(3)
!      KPP_REAL, PARAMETER :: ZERO = 0.0, ONE = 1.0
!
!      info = 0
!!      t = TINY(da)
!!      
!!      da = ABS(A(1,1)); L = 1
!!      IF ( ABS(A(2,1))>da ) THEN
!!        da = ABS(A(2,1)); L = 2
!!        IF ( ABS(A(3,1))>da ) THEN
!!          L = 3
!!        END IF  
!!      END IF  
!!      IPVT(1)  = L
!!      IF (L /=1 ) THEN
!!         TMP(1:3) = A(L,1:3)
!!         A(L,1:3) = A(1,1:3)
!!         A(1,1:3) = TMP(1:3)
!!      END IF
!!      IF (ABS(A(1,1)) < t) THEN
!!         info = 1
!!         return
!!      END IF   
!!
!      A(2,1) = A(2,1)/A(1,1)
!      A(2,2) = A(2,2) - A(2,1)*A(1,2)
!      A(2,3) = A(2,3) - A(2,1)*A(1,3)
!      A(3,1) = A(3,1)/A(1,1)
!      A(3,2) = A(3,2) - A(3,1)*A(1,2)
!      A(3,3) = A(3,3) - A(3,1)*A(1,3)
!      
!!      IPVT(2)  = 2
!!      IF (ABS(A(3,2))>ABS(A(2,2))) THEN
!!         IPVT(2)  = 3
!!         TMP(2:3) = A(3,2:3)
!!         A(3,2:3) = A(2,2:3)
!!         A(2,2:3) = TMP(2:3)
!!      END IF
!!      IF (ABS(A(2,2)) < t) THEN
!!         info = 1
!!         return
!!      END IF   
!!      
!      A(3,2)   = A(3,2)/A(2,2)
!      A(3,3)   = A(3,3) - A(3,2)*A(2,3)
!      IPVT(3)  = 3
!      
!END SUBROUTINE FAC3
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE SOL3(Trans,A,IPVT,b)
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!     SOL3 solves the system 3x3
!!     A * x = b  or  trans(a) * x = b
!!     using the factors computed by WGEFA.
!!
!!     Trans      = 'N'   to solve  A*x = b ,
!!                = 'T'   to solve  transpose(A)*x = b
!!     LINPACK - LIKE 
!!
!!     Remove comments to use pivoting
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!      CHARACTER     :: Trans
!      KPP_REAL :: a(3,3),b(3)
!      INTEGER       :: IPVT(3)
!!      INTEGER       :: L
!!      KPP_REAL :: TMP
!      
!      SELECT CASE (Trans)
!
!      CASE ('n','N')  !  Solve  A * x = b
!
!!     Solve  L*y = b
!!         L = IPVT(1)
!!         IF (L /= 1) THEN
!!            TMP = B(1); B(1) = B(L); B(L) = TMP
!!         END IF
!         b(2) = b(2)-A(2,1)*b(1)
!         b(3) = b(3)-A(3,1)*b(1)
!         
!!         L = IPVT(2)
!!         IF (L /= 2) THEN
!!            TMP = B(2); B(2) = B(L); B(L) = TMP
!!         END IF
!         b(3) = b(3)-A(3,2)*b(2)
!
!!     Solve  U*x = y
!         b(3) = b(3)/A(3,3)
!         b(2) = (b(2)-A(2,3)*b(3))/A(2,2)
!         b(1) = (b(1)-A(1,3)*b(3)-A(1,2)*b(2))/A(1,1)
!      
!      
!      CASE ('t','T')  !  Solve transpose(A) * x = b
!
!!      Solve transpose(U)*y = b
!         b(1) = b(1)/A(1,1)
!         b(2) = (b(2)-A(1,2)*b(1))/A(2,2)
!         b(3) = (b(3)-A(1,3)*b(1)-A(2,3)*b(2))/A(3,3)
!
!!      Solve transpose(L)*x = y
!         b(2) = b(2)-A(3,2)*b(3)
!!         L = ipvt(2)
!!         IF (L /= 2) THEN
!!            TMP = B(2); B(2) = B(L); B(L) = TMP
!!         END IF
!         b(1) = b(1)-A(3,1)*b(3)-A(2,1)*b(2)
!!         L = ipvt(1)
!!         IF (L /= 1) THEN
!!            TMP = B(1); B(1) = B(L); B(L) = TMP
!!         END IF
!   
!      END SELECT
!
!END SUBROUTINE SOL3

