pkg/radtrans/radtrans_solve_tridiag.F

b55e95f1ff Oliv*0001 #include "RADTRANS_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: RADTRANS_SOLVE_TRIDIAG
                0005
                0006 C     !INTERFACE:
                0007       SUBROUTINE RADTRANS_SOLVE_TRIDIAG(
                0008      U                     a3d, b3d, c3d,
                0009      U                     y3d,
                0010      I                     n, myThid )
                0011
                0012 C     !DESCRIPTION:
                0013 C     Solve a tri-diagonal system A*X=Y (dimension n) by Gaussian
                0014 C     elimination with pivoting
                0015 C     - input arrays are overwritten
                0016 C     - result is returned in y3d
                0017
                0018 C     !USES:
                0019       IMPLICIT NONE
                0020 #include "SIZE.h"
                0021
                0022 C     !INPUT PARAMETERS:
                0023       INTEGER n,myThid
                0024
                0025 C     !INPUT/OUTPUT PARAMETERS:
                0026 C     a3d(2:n)   :: matrix lower diagnonal
                0027 C     b3d(1:n)   :: matrix main  diagnonal
                0028 C     c3d(1:n-1) :: matrix upper diagnonal
                0029 C     y3d(1:n)   :: Input = Y vector ; Output = X = solution of A*X=Y
                0030       _RL a3d(2*Nr)
                0031       _RL b3d(2*Nr)
                0032       _RL c3d(2*Nr)
                0033       _RL y3d(2*Nr)
                0034 CEOP
                0035
                0036 #ifdef ALLOW_RADTRANS
                0037
                0038 C     !LOCAL VARIABLES:
                0039 C     d3d(1:Nr-2) :: second upper diagnonal, generated by pivoting
                0040       INTEGER k
                0041       _RL tmpc, tmpy
                0042       _RL recVar
                0043       _RL d3d(2*Nr)
                0044
                0045       c3d(n) = 0.
                0046
                0047 c eliminate lower diagonal
                0048 c     only c3d, d3d and y3d are modified and will be used later
                0049       DO k=1,n-1
                0050         IF(ABS(b3d(k)).GE.ABS(a3d(k+1)))THEN
                0051 c         scale current row, subtract from next
                0052           recVar = 1. / b3d(k)
                0053           c3d(k) = c3d(k)*recVar
                0054           y3d(k) = y3d(k)*recVar
                0055           b3d(k+1) = b3d(k+1) - a3d(k+1)*c3d(k)
                0056           y3d(k+1) = y3d(k+1) - a3d(k+1)*y3d(k)
                0057           d3d(k) = 0.
                0058         ELSE
                0059 c         swap rows, then scale current and subtract from next
                0060           recVar = 1. / a3d(k+1)
                0061           tmpc = c3d(k)
                0062           tmpy = y3d(k)
                0063           c3d(k) = b3d(k+1)*recVar
                0064           d3d(k) = c3d(k+1)*recVar
                0065           y3d(k) = y3d(k+1)*recVar
                0066           b3d(k+1) = tmpc - b3d(k)*c3d(k)
                0067           c3d(k+1) =      - b3d(k)*d3d(k)
                0068           y3d(k+1) = tmpy - b3d(k)*y3d(k)
                0069         ENDIF
                0070       ENDDO
                0071       recVar = 1. / b3d(n)
                0072       y3d(n) = y3d(n)*recVar
                0073
                0074 c backsubstitute
                0075 c     y3d(n) is already good
                0076       y3d(n-1) = y3d(n-1) - c3d(n-1)*y3d(n)
                0077       DO k=n-2,1,-1
                0078         y3d(k) = y3d(k) - c3d(k)*y3d(k+1) - d3d(k)*y3d(k+2)
                0079       ENDDO
                0080
                0081 #endif /*ALLOW_RADTRANS */
                0082
                0083       RETURN
                0084       END
                0085