戻る


PROGRAM MAIN
C ********************************************************************* C
C C
C CVDF1D - FDM C
C C
C 1-D - Convection-Diffusion Equation Solver C
C C
C by Finite Difference Method C
C C
C ( Implicit Method & 7 Explicit methods ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ********************************************************************* C
C
IMPLICIT REAL*8(A-H,O-Z)

DIMENSION R ( 101), T ( 101), TD ( 101), TX( 101), X( 101),
1 T0( 101), TA1( 101), TN1( 101), A ( 5,101)
C
COMMON / CS2 / Bx, Rx
COMMON / CST / U, DX, RKV, DT, LOP, NM, TMX
C
OPEN ( 7, FILE='CVDF1_OUT.DATA' )
OPEN ( 8, FILE='CVDF1_ANL.DATA' )
C
C ----------------------------
C
CALL TTLE
C
C ----------------------------
C
CALL SBIP ( NNP, *999 ) ! Ex. NNP= 11, U= 0.1, K= 0.01, DT= 0.4
C
C ----------------------------
C
NM = 0
C
C ----- Implicit FDM ----------------------------------------
C
CALL CNMT ( R, T, TD, TX, A, X, T0, TA1, TN1,NNP )
C
C -----------------------------------------------------------
C
999 NM = NM + 1
C
C ----- Explicit FDMs ---------------------------------------------
C
CALL CEFD ( R, T, TD, TX, A, X, T0, TA1, TN1,
1 NNP, NM, *999, *900 )
C
C -----------------------------------------------------------------
C
900 CONTINUE
C
CLOSE( 7)
CLOSE( 8)
C
STOP
END
C **********************************************************************
C
SUBROUTINE SBIP ( NNP, * )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON / CS2 / Bx, Rx
COMMON / CST / U, DX, RKV, DT, LOP, NM, TMX
C
C ----- Ex. NNP= 11, U= 0.1, K= 0.01, DT= 0.1 ---
C
WRITE(6,*) ' '
WRITE(6,*) ' *** Please Input "Node Numbers" '
READ(5,*) NNP
WRITE(6,*) ' *** Please Input "U" ! '
** READ(5,*) U
U = 0.1D0
C
WRITE(6,*) ' '
WRITE(6,*) ' Please Input "DT" ! '
READ(5,*) DT
WRITE(6,*) ' Please Input "KX" ! '
** READ(5,*) RKV
RKV = 0.01D0
C
C ***** Attention ***
C
LOP = 0.065D0/ ( ( RKV + 0.13D0* U )* DT) + 0.0001D0
C --------------------------------------------------------
C
IF ( LOP.EQ.0) WRITE(6,*) ' **** Loop= 0 & STOP ! ***'
IF ( LOP.EQ.0) RETURN 1
C
DX = 2.D0/ DFLOAT( NNP-1)
C
Bx = U* DT/ DX
C -------------------
C
EL = 20.D0
C ***************
C
IF ( RKV.LT.1.D-10) RKV = 1.D-10
C
Rx = RKV* DT/( DX* DX)
C ---------------------------
C
WRITE(6,*) ' '
WRITE(6,*) '*** Input Data ***'
WRITE(6,2000)NNP, U, DT, RKV, DX, Bx, Rx, LOP
C
C --------------------------------------------------
WRITE(7,2100)NNP, U, DT, RKV, DX, Bx, Rx, LOP
2100 FORMAT(I5,6F7.3,I3)
C --------------------------------------------------
C
2000 FORMAT(/,' * No. of Nodes =',I4,' U =',F8.4,' DT =',F8.4,
1 ' Kx =',F8.4,/,26X,' DX =',F8.4, ' Bx =',F8.4,' Rx =',
2 F8.4,/,' * LOOP =', I4 )
C
TMX = DT* DFLOAT( LOP)
C
WRITE(6,2200) TMX
2200 FORMAT(' * Final Solution : Time = ',F5.3)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CEFD ( R, T, TD, TX, A, X, T0, TA1, TN1, NNP, *, * )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON / CS2 / Bx, Rx
COMMON / CST / U, DX, RKV, DT, LOP, NM, TMX
C
C ----- Explicit FDM Methods ---
C
DIMENSION R( NNP), T ( NNP), TD ( NNP), TX( NNP), A( 5,NNP),
1 X( NNP), T0( NNP), TA1( NNP), TN1( NNP)
C
NM = NM
C **************
C
TMX = DT* DFLOAT( LOP)
C
C ----- To initialize 'T' & evaluate 'TX' ------------
C
CALL EXSL ( NNP, X, T, TX, DX, U, RKV, TMX )
C
C -----------------------------------------------------
C
DO 100 J = 1,NNP
DO 100 K = 1, 5
C
A(K,J) = 0.D0
100 CONTINUE
C
C ---------------------------------------------------------------------
DO 200 N = 1, LOP
C
C ------------------------------------------------------
C
CALL SEFDM ( T, TN1, U, DT, RKV, DX,NNP, NM )
C
C ------------------------------------------------------
DO 220 J = 1,NNP
C
IF ( T( J).GE.2.D0.OR.T( J).LE.-1.D0 ) GO TO 990
220 CONTINUE
C -----------------------------------------------------
IF ( DABS( DFLOAT( N-LOP)/2.D0 ).GT.0.55D0) GO TO 200
C
T1 = DT* DFLOAT( N)
C
C ------------------------------------------------------------------
C
CALL EXSL ( NNP, X, T0, TA1, DX, U, RKV, T1 )
C
C ------------------------------------------------------------------
C
DO 240 J = 1,NNP
C
TN1( J) = T( J)
240 CONTINUE
C
200 CONTINUE
C ---------------------------------------------------------------------
C
** WRITE(6,2000) TMX
** 2000 FORMAT(/,' ** Final Solution * TIME = ',F5.3)
** WRITE(6,2100) ( T( J),J = 1,NNP)
** 2100 FORMAT('* Num_R.=',1X,10F6.3,1X,10F6.3,1X,10F6.3,
** 1 1X,10F6.3,1X,10F6.3,1X,10F6.3,1X,10F6.3)
C
C ------------------------------------------------
SUM = 0.D0
DO 260 J = 2, NNP-1
C
SUM = SUM + ( T( J)-TX( J) )** 2
260 CONTINUE
C ------------------------------------------------
C
RMS = DSQRT( SUM/ DFLOAT( NNP-2) )
C
IF ( NM.EQ.1) WRITE(6,2001) RMS*100.D0
2001 FORMAT(/,' * FTCS Method * RMS*100=', F10.3)
IF ( NM.EQ.2) WRITE(6,2002) RMS*100.D0
2002 FORMAT(/,' * Lax-Wendroff Method * RMS*100=', F10.3)
IF ( NM.EQ.3) WRITE(6,2003) RMS*100.D0
2003 FORMAT(/,' * 1st Order Upwind Method * RMS*100=', F10.3)
IF ( NM.EQ.4) WRITE(6,2004) RMS*100.D0
2004 FORMAT(/,' * 4-Ppoint Upwind Method * RMS*100=', F10.3)
IF ( NM.EQ.5) WRITE(6,2005) RMS*100.D0
2005 FORMAT(/,' * Quick Method * RMS*100=', F10.3)
IF ( NM.EQ.6) WRITE(6,2006) RMS*100.D0
2006 FORMAT(/,' * 3rd ord. Upwind Method * RMS*100=', F10.3)
IF ( NM.EQ.7) WRITE(6,2007) RMS*100.D0
2007 FORMAT(/,' * 4th ord. Centered Method * RMS*100=', F10.3)
C
C -----------------------------------------------
WRITE(7,2300) NM, RMS
2300 FORMAT(I3,E12.5)
WRITE(7,2400) TMX, ( T( J), J = 1, NNP )
2400 FORMAT(10F7.3)
C -----------------------------------------------
C
IF ( NM.GE.7) THEN
WRITE(6,*) ' '
WRITE(6,*) '*** Simulation End ***'
RETURN 2
END IF
C
RETURN 1
C
990 CONTINUE
WRITE(6,*) ' *** Divergence ! ***'
WRITE(6,*) ' ** T( I) = ', ( T( J), J = 1, NNP )
C
RETURN 2
C
END
C **********************************************************************
C
SUBROUTINE CNMT ( R, T, TD, TX, A, X, T0, TA1, TN1, NNP )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON / CS2 / Bx, Rx
COMMON / CST / U, DX, RKV, DT, LOP, NM, TMX
C
C ----- Crank - Nicolson Method ---
C
DIMENSION R ( NNP), T( NNP), TD( NNP), TX ( NNP),
1 A ( 5,NNP), X( NNP), T0( NNP), TA1( NNP),
2 TN1( NNP)
C
TMX = DT* DFLOAT( LOP)
C
C ------------------------------------------------------------------
C
CALL EXSL ( NNP, X, T, TX, DX, U, RKV, TMX )
C
C ----------------------------------------------------------
DO 100 J = 1,NNP
DO 100 K = 1, 5
C
A( K,J) = 0.D0
100 CONTINUE
C ----------------------------------------------------------
C
WRITE(6,2000)
2000 FORMAT(/,' *** Initial Conditions *',/)
WRITE(6,2100) ( X( J),J = 1,NNP)
2100 FORMAT(' * Xi =',1X,10F6.3,1X,10F6.3,1X,10F6.3,1X,10F6.3,
1 1X,10F6.3,1X,10F6.3,1X,10F6.3)
WRITE(6,2200) ( T( J),J = 1,NNP)
2200 FORMAT(' * Ci =',1X,10F6.3,1X,10F6.3,1X,10F6.3,1X,
1 10F6.3,1X,10F6.3,1X,10F6.3,1X,10F6.3)
WRITE(6,2300) ( TX( J),J = 1,NNP)
2300 FORMAT(/,' * Anl_R.=',1X,10F6.3,1X,10F6.3,1X,10F6.3,1X,
1 10F6.3,1X,10F6.3,1X,10F6.3,1X,10F6.3)
WRITE(6,*) ' '
C
C ----------------------------------------------------------
DO 200 N = 1, LOP
C
C ----- Tridiagonal system for implicit schemes ---
C
C -------------------------------------------------
DO 220 J = 2, NNP-1
C
A( 1,J-1) = 0.D0
A( 2,J-1) = - 0.25D0*Bx - 0.5D0*Rx
A( 3,J-1) = 1.D0 + Rx
A( 4,J-1) = 0.25D0*Bx - 0.5D0*Rx
C
R( J-1) = ( 0.25D0* Bx + 0.5D0* Rx )* T( J-1)
1 + ( 1.D0 - Rx )* T( J)
2 + ( - 0.25D0* Bx + 0.5D0* Rx )* T( J+1)
C
220 CONTINUE
C -------------------------------------------------
C
R( 1) = R( 1) - A( 2,1)* T( 1)
C
A( 2, 1) = 0.D0
A( 4, NNP-2) = 0.D0
C
C --------------------------------------
C
CALL BNFC ( A, NNP-2, 1)
C
C --------------------------------------
C
CALL BNSL ( R, TD, A, NNP-2, 1)
C
C ------------------------------------------------
DO 250 J = 2, NNP-1
C
T( J) = TD( J-1)
250 CONTINUE
C -------------------------------------------------
DO 270 J = 1, NNP
C
IF ( T( J).GE.2.D0.OR.T( J).LE.-1.D0) GO TO 590
270 CONTINUE
C -------------------------------------------------
C
IF ( DABS( DFLOAT( N-LOP)/ 2.D0 ).GT.0.55D0) GO TO 200
C
T1 = DT* DFLOAT( N)
C
C ------------------------------------------------------------------
C
CALL EXSL ( NNP, X, T0, TA1, DX, U, RKV, T1 )
C
C ------------------------------------------------------------------
DO 290 J = 1,NNP
C
TN1( J) = T( J)
290 CONTINUE
C ------------------------------------------------
C
200 CONTINUE
C ----------------------------------------------------------
** WRITE(6,*) ' '
** WRITE(6,2400) TMX
** 2400 FORMAT(' ** Final Solution * TIME = ',F5.3)
** WRITE(6,2500) ( T( J),J = 1,NNP)
** 2500 FORMAT(' * Num_R.=',1X,10F6.3,1X,10F6.3,1X,10F6.3,
** 1 1X,10F6.3,1X,10F6.3,1X,10F6.3,1X,10F6.3)
C ----------------------------------------------------------
C
SUM = 0.D0
C
DO 300 J = 2, NNP-1
C
SUM = SUM + ( T( J)-TX( J))** 2
300 CONTINUE
C
RMS = DSQRT( SUM/ ( DFLOAT( NNP-2)) )
C
WRITE(6,*) '*** Implicit Method ***'
WRITE(6,2600) RMS*100.D0
2600 FORMAT(/,' * Crank-Nicolson Method',7X,' RMS*100=',F10.3)
C
C ------------------------------------------------
WRITE(7,2700) RMS
2700 FORMAT(E12.5)
WRITE(7,2800) ( T ( J), J = 1, NNP )
WRITE(8,2800) ( TX( J), J = 1, NNP )
2800 FORMAT(10F7.3)
C ------------------------------------------------
C
WRITE(6,*) ' '
WRITE(6,*) '*** Explicit Method ***'
C
RETURN
C
590 CONTINUE
WRITE(6,*) '*** Divergence ! ***'
WRITE(6,*) ' * T( I) =', ( T( J), J = 1, NNP )
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,100)
WRITE(6,101) 'C *********************************************** C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C C V D F (1D) - FDM C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C ( V.3.5 ) C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C 1D - Convection-Diffusion Equation Solver C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C by F.D.M. C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C Copyright 2011 : Yasuhiro MATSUDA C'
WRITE(6,101) 'C C'
WRITE(6,101) 'C *********************************************** C'
C
WRITE(6,*) ' '
WRITE(6,*) ' ** Implicit Method : Crank-Nicolson Method '
WRITE(6,*) ' '
WRITE(6,*) ' ** Explicit Methods :'
WRITE(6,*) ' '
WRITE(6,*) ' - FTCS M. - Lax-Wendroff M. - 1
1st ord. Upwind M.'
WRITE(6,*) ' - 4-Pts. Upwind M. - QUICK M. - 3
1rd ord. Upwind M.'
WRITE(6,*) ' - 4th ord. Centered M. '
C
100 FORMAT(/)
101 FORMAT(12X,A55)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SEFDM ( TN1, TN, U, DT, RKX, DX, NNP, NM )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON / CS2 / Bx, Rx
C
DIMENSION TN1( NNP), TN( NNP)
C
C ----- Explicit Methods ---
C
C -------------------------
DO 100 I = 1, NNP
C
TN( I) = TN1( I)
100 CONTINUE
C ----------------------------------------------------------
DO 200 I = 1, NNP
C
IF ( I.EQ.1 ) GO TO 123
IF ( I.EQ.NNP) GO TO 300
C
IF ( I.EQ.NNP-1.AND.NM.EQ.6 ) GO TO 145
IF ( I.EQ.NNP-1.AND.NM.EQ.7 ) GO TO 165
IF ( I.NE.2) GO TO 70
C
GO TO ( 70, 70, 70, 80, 120, 140, 160 ) NM
C
70 CONTINUE
C
IF ( NM.EQ.1) ! FTCS
C
1 TN1( I) = TN( I) - Bx* (TN( I+1) - TN( I-1)) / 2.D0
2 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
IF ( NM.EQ.2) ! Lax-Wendroff
C
1 TN1( I) = TN( I) - Bx* (TN( I+1) - TN( I-1) ) / 2.D0
2 + Bx** 2* ( TN( I+1) - 2.D0* TN( I)+ TN( I-1))/ 2.D0
3 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
IF ( NM.EQ.3) ! !st ordered Upwind
C
1 TN1( I) = TN( I) - Bx* ( TN( I) - TN( I-1) )
2 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
IF ( NM.EQ.4) ! 4 pts. Upwind
C
1 TN1( I) = TN( I) - Bx* ( 2.D0* TN( I+1) + 3.D0* TN( I)
2 - 6.D0* TN( I-1) + TN( I-2) ) / 6.D0
3 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
IF ( NM.EQ.5) ! QUICK
C
1 TN1( I) = TN( I) - Bx * ( 3.D0* TN( I+1) + 3.D0* TN( I)
2 - 7.D0* TN( I-1) + TN( I-2) ) / 8.D0
3 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
IF ( NM.EQ.6) ! #rd ordered Upwind
C
1 TN1( I) = TN( I) - Bx* ( - TN ( I+2) + 12.D0* TN( I+1)
2 + 6.D0* TN( I) - 20.D0* TN( I-1) + 3.D0* TN( I-2)) / 24.D0
3 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
IF ( NM.EQ.7) ! 4th ordered Center
C
1 TN1( I) = TN( I) - Bx* ( - TN( I+2) + 8.D0* TN( I+1)
2 - 8.D0* TN( I-1) + TN( I-2) ) / 12.D0
3 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
GO TO 200
C
C -------- Boundary Conditions ---
C
123 TN1( I) = 1.D0
C
GO TO 200
C
300 TN1( I) = 0.D0
C
GO TO 200
C
80 TN1( I) = TN( I) - Bx* ( 2.D0* TN( I+1) + 3.D0* TN( I)
1 - 6.D0* TN( I-1) + 1.D0 ) / 6.D0
2 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
GO TO 200
C
120 TN1( I) = TN( I) - Bx* ( 3.D0* TN( I+1) + 3.D0* TN( I)
1 - 7.D0* TN( I-1) + 1.D0 ) / 8.D0
2 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
GO TO 200
C
C ** B. C. *** for ' 6' ***
* 1 TN1( I) = TN( I)
* 2 - Bx * ( - TN ( I+2) + 12.D0* TN( I+1)
* 3 + 6.D0* TN( I) - 20.D0* TN( I-1)
* 4 + 3.D0* TN( I-2) ) / 24.D0 + RKX * DT/ DX / DX *
* 5 ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
140 TN1( I) = TN( I) - Bx* ( - TN( I+2) + 12.D0* TN( I+1)
1 + 6.D0* TN( I) - 20.D0* TN( I-1) + 3.D0 ) / 24.D0
2 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
GO TO 200
C
145 TN1( I) = TN( I) - Bx* ( 12.D0* TN( I+1) + 6.D0* TN( I)
1 - 20.D0* TN( I-1) + 3.D0* TN( I-2))/ 24.D0
2 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
GO TO 200
C
160 TN1( I) = TN( I) - Bx* ( - TN( I+2) + 8.D0* TN( I+1)
1 - 8.D0* TN( I-1) + 1.D0 )/ 12.D0
2 + Rx* ( TN( I+1) - 2.D0* TN( I) + TN( I-1) )
C
GO TO 200
C
165 TN1( I) = TN( I)
1 - Bx* ( 8.D0* TN( I+1) - 8.D0* TN( I-1) + TN( I-2) ) / 12.D0
2 + Rx* ( TN( I+1) - 2.D0*TN( I) + TN( I-1) )
C
GO TO 200
C
200 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE EXSL ( NNP, X, T, TX, DX, U, RKV, TMX )
C
C ----- To set initial 'T' sol. & Final Exact( TX) sol. ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION T( NNP), TX( NNP), X( NNP)
C
PI = 3.141592653589793D0
EL = 20.D0
C
XST = -1.D0
C
C ------------------------------------------------
DO 100 J = 1, NNP
C
X ( J) = XST + DX* DFLOAT( J-1)
T ( J) = 0.D0
TX( J) = 0.D0
C
100 CONTINUE
C ------------------------------------------------
C
C ----- Exact solution fior propagating temperature front ---
C
T ( 1) = 1.D0
TX( 1) = 1.D0
C
NEX = 50000
C
C ----------------------------------------------------------
DO 200 J = 2, NNP-1
C
IF ( X( J).LT. 0.D0) T( J) = 1.D0
IF ( DABS( X( J)).LT.0.0001D0) T( J) = 0.5D0
C
DEM = 0.D0
C
C ------------------------------------------------
DO 250 K = 1, NEX
C
VK = 2* K - 1
DUM = VK* PI/ EL
SNE = DSIN( DUM* ( X( J) - U* TMX) )
DIM = - RKV* DUM* DUM* TMX
C
IF ( DIM.LT.-20.D0) GO TO 220
C
DEM = DEM + ( SNE/ VK )* DEXP( DIM)
C
250 CONTINUE
C ------------------------------------------------
C
220 TX( J) = 0.5D0 - 2.D0* DEM / PI
C
200 CONTINUE
C
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE BNSL ( R, X, B, N, INT )
C
C ----- To use L & U factorisation to solve for X ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION R( 101), X( 101), B( 5,101)
C
IF ( INT.EQ.2) GO TO 3
C
C ----- INT = 1 : Tridiagonal system ---
C
NP = N - 1
C
C ------------------------------------------------
DO 100 J = 1, NP
C
JP = J + 1
R( JP) = R( JP) - B( 2,JP)* R( J)
C
100 CONTINUE
C ------------------------------------------------
X( N) = R( N)/ B( 3,N)

DO 200 J = 1, NP
C
JA = N - J
X( JA) = ( R( JA) - B( 4,JA)* X( JA+1) ) / B( 3,JA)
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
C
3 IBC = 1
NH = N / 2
C
IF ( 2* NH.EQ.N) R( N+1) = 0.D0
IF ( 2* NH.EQ.N) B( 2,N+1) = 0.D0
C
C ------------------------------------------------
DO 300 I = 1, 2
DO 300 J = 1, NH
C
JA = 2* J
C
C --------------------------------------
DO 340 K = 1, I
C
JB = JA - 1 + K
R( JB) = R( JB) - B( I,JB)* R( JB-1)
C
340 CONTINUE
C --------------------------------------
C
300 CONTINUE
C ------------------------------------------------
C
NEN = NH - IBC
C
X( N) = R( N)/ B( 3,N)
C
IF ( IBC.EQ.1) X( N-1)=( R( N-1) - B( 4,N-1)* X( N))/ B( 3,N-1)
C
C ------------------------------------------------
DO 400 J = 1, NEN
C
JA = N - 2*J + 1 - IBC
C
X( JA) = ( R( JA) - B( 4,JA)* X( JA+1) - B( 5,JA)* X( JA+2) )
1 / B( 3,JA)
C
X(JA-1) = ( R( JA) - B( 4,JA-1) * X( JA) )/ B( 3,JA-1)
C
400 CONTINUE
C ------------------------------------------------
C
RETURN
END
C *********************************************************************
C
SUBROUTINE BNFC ( B, N, INT )
C
C ----- To factorises Band Matrix arising Linear or
C quadratic into L, V ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION B( 5,101)
C
IF( INT.EQ.2) GO TO 150
C
C ----- INT = 1 : Linear elements : Tridiagonal system ---
C
NP = N - 1
C
C ------------------------------------------------
DO 100 J = 1, NP
C
JP = J + 1
C
B( 2,JP) = B( 2,JP)/ B( 3,J)
B( 3,JP) = B( 3,JP) - B( 2,JP)* B( 4,J)
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
C
C ----- INT = 2 : Quadratic elements: Tridaiagonal system ---
C
150 NH = N / 2
C
C ----------------------------------------------------------
DO 200 I = 1, 2
C
JS = 3 - I
C
C ------------------------------------------------
DO 250 J = JS, NH
C
JA = 2* ( J-1)
C
IF ( I.EQ.2) GO TO 270
C
JB = JA + 2
C
B( 1,JB) = B( 1,JB)/ B( 2,JB-1)
B( 2,JB) = B( 2,JB) - B( 1,JB)* B( 3,JB-1)
B( 3,JB) = B( 3,JB) - B( 1,JB)* B( 4,JB-1)
C
GO TO 250
C
270 JB = JA + 3
B( 2,JB-1) = B( 2,JB-1)/ B( 3,JB-2)
B( 3,JB-1) = B( 3,JB-1) - B( 2,JB-1)* B( 4,JB-2)
C
IF ( JB.GT.N) GO TO 250
C
B( 2,JB) = B( 2,JB) / B( 3,JB-1)
B( 3,JB) = B( 3,JB) - B( 2,JB)* B( 4,JB-1)
B( 4,JB) = B( 4,JB) - B( 2,JB)* B( 5,JB-1)
C
250 CONTINUE
C ------------------------------------------------
C
200 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************