trid3d.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE trid3d
!                   *****************
!
     &(xaux,x,b,npoin,npoin2)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    SOLVES TRIDIAGONAL SYSTEMS FOR EVERY VERTICAL
!+                IN A MESH OF PRISMS.
!
!history  J-M HERVOUET (LNHE)
!+        30/09/05
!+        V5P6
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| B              |-->| RIGHT-HAND SIDE OF SYSTEM
!| NPOIN          |-->| NUMBER OF POINTS
!| NPOIN2         |-->| NUMBER OF POINTS IN 2D
!| X              |<--| SOLUTION OF SYSTEM
!| XAUX           |-->| TRI-DIAGONAL MATRIX
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_trid3d => trid3d
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NPOIN,NPOIN2
!
      DOUBLE PRECISION, INTENT(IN)    :: B(npoin2,*)
      DOUBLE PRECISION, INTENT(INOUT) :: XAUX(npoin,*),X(npoin2,*)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,IPLAN,I3D,NPLAN
      DOUBLE PRECISION EPS
!
      INTRINSIC abs
!
!-----------------------------------------------------------------------
!
!     XAUX(I,1) IS COEFFICIENT OF POINT BELOW I IN EQUATION OF POINT I
!     XAUX(I,2) IS THE DIAGONAL
!     XAUX(I,3) IS COEFFICIENT OF POINT ABOVE I IN EQUATION OF POINT I
!
!     XAUX(I,4) HERE USED AS WORKING ARRAY
!     XAUX(I,5) HERE USED AS WORKING ARRAY
!
!-----------------------------------------------------------------------
!
      nplan=npoin/npoin2
      eps=1.d-8
!
!-----------------------------------------------------------------------
!
!     BASIC ALGORITHM TAKEN IN "NUMERICAL RECIPES" PAGE 40 AND ADAPTED
!     TO STORAGE
!
!     XAUX(*,4) : WORKING ARRAY (SIZE NPOIN2)
!     XAUX(*,5) : WORKING ARRAY (SIZE NPOIN)
!
      DO i=1,npoin2
        xaux(i,4)=xaux(i,2)
        IF(abs(xaux(i,4)).LT.eps) THEN
          WRITE(lu,*) 'TRID3D: SYSTEM ILL-DEFINED'
          CALL plante(1)
          stop
        ENDIF
        x(i,1)=b(i,1)/xaux(i,4)
      ENDDO
!
      DO iplan=2,nplan
      DO i=1,npoin2
        i3d=i+npoin2*(iplan-1)
        xaux(i3d,5)=xaux(i3d-npoin2,3)/xaux(i,4)
        xaux(i,4)=xaux(i3d,2)-xaux(i3d,1)*xaux(i3d,5)
        IF(abs(xaux(i,4)).LT.eps) THEN
          WRITE(lu,*) 'TRID3D: SYSTEM ILL-DEFINED'
          WRITE(lu,*) '        PRECONDITIONING 17 IMPOSSIBLE'
          CALL plante(1)
          stop
        ENDIF
        x(i,iplan)=(b(i,iplan)-xaux(i3d,1)*x(i,iplan-1))/xaux(i,4)
      ENDDO
      ENDDO
!
      DO iplan=nplan-1,1,-1
      DO i=1,npoin2
        i3d=i+npoin2*iplan   ! PLAN TO THE TOP
        x(i,iplan)=x(i,iplan)-xaux(i3d,5)*x(i,iplan+1)
      ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END