*DECK HEAT02
      SUBROUTINE INITAL
C...
C...  METHOD OF LINES INTEGRATION OF BOUNDARY-VALUE ODES
C...
C...  THE FOLLOWING PROGRAMMING ILLUSTRATES HOW A SYSTEM OF FOUR
C...  BOUNDARY-VALUE ORDINARY DIFFERENTIAL EQUATIONS (BVODES) CAN
C...  BE INTEGRATED NUMERICALLY BY FIRST CONVERTING THE BVODES TO
C...  INITIAL-VALUE PDES, THEN INTEGRATING THE PDES TO STEADY STATE.
C...  THE STEADY STATE SOLUTION OF THE PDES IS THE SOLUTION TO THE
C...  BVODES.
C...
C...  THE FOLLOWING PARTIAL DIFFERENTIAL EQUATIONS (PDES) MODEL A
C...  THREE-PASS SHELL AND TUBE HEAT EXCHANGER
C...
C...   T1  + VT*T1  = CT*(TS - T1)
C...     T        X
C...
C...   T2  - VT*T2  = CT*(TS - T2)
C...     T        X
C...
C...   T3  + VT*T3  = CT*(TS - T3)
C...     T        X
C...
C...   TS  - VS*TS  = CS*((T1 - TS) + (T2 - TS) + (T3 - TS))
C...     T        X
C...
C...  WHERE
C...
C...   T1,T2,T3   TEMPERATURE OF THE FLUID IN TUBE PASSES 1, 2 AND 3,
C...              RESPECTIVELY
C...
C...   TS         TEMPERATURE OF THE FLUID IN THE SHELL
C...
C...   X          POSITION ALONG THE EXCHANGER
C...
C...   T          TIME
C...
C...   VT,VS      AVERAGE VELOCITY OF THE TUBE AND SHELL SIDE FLUIDS,
C...              RESPECTIVELY
C...
C...   CT,CS      CONSTANTS DEFINING THE HEAT TRANSFER RATES TO THE
C...              TUBE SIDE FLUID AND SHELL SIDE FLUID, RESPECTIVELY
C...
C...  THE EXCHANGER IS TO BE SIMULATED SUBJECT TO THE FOLLOWING CON-
C...  DITIONS AND PARAMETERS:
C...
C...   T1(X,0) = T2(X,0) = T3(X,0) = TS(X,0) = 25
C...
C...   T1(0,T) = 25, T2(XL,T) = T1(XL,T), T3(0,T) = T2(0,T)
C...
C...   TS(XL,T) = 220
C...
C...   VT = 1, VS = 0.5, CT = CS = 0.1, XL = 10, N = 11
C...
C...  THE EXITING TEMPERATURES, T3(XL,T) AND TS(0,T), ARE OF PARTICULAR
C...  INTEREST IN THIS SIMULATION.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=11)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    T1(NX),    T2(NX),    T3(NX),     TS(NX)
     2      /F/   T1T(NX),   T2T(NX),   T3T(NX),    TST(NX)
     3      /S/   T1X(NX),   T2X(NX),   T3X(NX),    TSX(NX)
     4      /C/        T0,        XL,        VT,         VS,
     5                 CT,        CS,       TS0,          N,
     6             XP(NX),        IP,       ISP
C...
C...  SET THE MODEL PARAMETERS
      T0=25.0D0
      XL=10.0D0
      VT=1.0D0
      VS=0.5D0
      CT=0.1D0
      CS=0.1D0
      TS0=220.0D0
      N=NX
C...
C...  INITIALIZE THE SPATIAL GRID USED IN SUBROUTINE PRINT
      DO 1 I=1,NX
         XP(I)=XL*DFLOAT(I-1)/DFLOAT(NX-1)
1     CONTINUE
C...
C...  MODEL INITIAL CONDITIONS
      DO 2 I=1,NX
         T1(I)=T0
         T2(I)=T0
         T3(I)=T0
         TS(I)=T0
2     CONTINUE
C...
C...  INITIALIZE THE DERIVATIVES
      CALL DERV
      IP=0
      ISP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=11)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    T1(NX),    T2(NX),    T3(NX),     TS(NX)
     2      /F/   T1T(NX),   T2T(NX),   T3T(NX),    TST(NX)
     3      /S/   T1X(NX),   T2X(NX),   T3X(NX),    TSX(NX)
     4      /C/        T0,        XL,        VT,         VS,
     5                 CT,        CS,       TS0,          N,
     6             XP(NX),        IP,       ISP
C...
C...  BOUNDARY CONDITIONS
      T1( 1)=T0
      TS(NX)=TS0
      T2(NX)=T1(NX)
      T3( 1)=T2(1)
C...
C...  SPATIAL DERIVATIVES
      CALL DSS020(0.0D0,XL,NX,T1,T1X, 1.0D0)
      CALL DSS020(0.0D0,XL,NX,T2,T2X,-1.0D0)
      CALL DSS020(0.0D0,XL,NX,T3,T3X, 1.0D0)
      CALL DSS020(0.0D0,XL,NX,TS,TSX,-1.0D0)
C...
C...  MODEL PDES
      DO 1 I=1,N
         T1T(I)=-VT*T1X(I)+CT*(TS(I)-T1(I))
         T2T(I)= VT*T2X(I)+CT*(TS(I)-T2(I))
         T3T(I)=-VT*T3X(I)+CT*(TS(I)-T3(I))
         TST(I)= VS*TSX(I)+CS*((T1(I)-TS(I))
     1         +(T2(I)-TS(I))+(T3(I)-TS(I)))
1       CONTINUE
C...
C...  ZERO THE TEMPORAL DERIVATIVES AT THE BOUNDARIES
      T1T(1)=0.0D0
      T2T(N)=0.0D0
      T3T(1)=0.0D0
      TST(N)=0.0D0
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=11)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    T1(NX),    T2(NX),    T3(NX),     TS(NX)
     2      /F/   T1T(NX),   T2T(NX),   T3T(NX),    TST(NX)
     3      /S/   T1X(NX),   T2X(NX),   T3X(NX),    TSX(NX)
     4      /C/        T0,        XL,        VT,         VS,
     5                 CT,        CS,       TS0,          N,
     6             XP(NX),        IP,       ISP
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION   T3P(101), TSP(101),   TP(101),
     1          T1SP(7,NX),          T2SP(7,NX),
     2          T3SP(7,NX),          TSSP(7,NX)
C...
C...  MONITOR THE CALCULATION
      WRITE(*,*)NORUN,T
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IF(IP.EQ.0)WRITE(NO,1)
1     FORMAT(9X,1HT,4X,8H T1(0,T),4X,8HTS(XL,T),
     1              4X,8HT3(XL,T),4X,8H TS(0,T))
C...
C...  PRINT THE SOLUTION
      IP=IP+1
      IF((IP-1)/10*10.EQ.(IP-1))THEN
         WRITE(NO,2)T,T1(1),TS(NX),T3(NX),TS(1)
2        FORMAT(F10.1,4F12.1)
      END IF
C...
C...  STORE THE SOLUTION FOR PLOTTING
      TP(IP)=T
      T3P(IP)=T3(NX)
      TSP(IP)=TS(1)
C...
C...  THE SPATIAL PROFILES ARE STORED EVERY 20TH CALL TO PRINT
      IF((IP-1)/20*20.EQ.(IP-1))THEN
      ISP=ISP+1
      DO 5 I=1,NX
         T1SP(ISP,I)=T1(I)
         T2SP(ISP,I)=T2(I)
         T3SP(ISP,I)=T3(I)
         TSSP(ISP,I)=TS(I)
5     CONTINUE
      END IF
C...
C...  PLOT THE SOLUTION AT THE END OF THE RUN.  THE ENTERING SHELL
C...  FLUID TEMPERATURE IS FIRST STORED TO SCALE THE TEMPERATURE PLOTS
4     IF(IP.LT.101)RETURN
      DO 6 I=1,NX
         T1SP(7,I)=TS0
         T2SP(7,I)=TS0
         T3SP(7,I)=TS0
         TSSP(7,I)=TS0
6     CONTINUE
      CALL SPLOTS(1,IP,TP,T3P)
      WRITE(NO,10)
10    FORMAT(1H ,//,31H T3(XL,T) VS T                 )
      CALL SPLOTS(1,IP,TP,TSP)
      WRITE(NO,11)
11    FORMAT(1H ,//,31H TS(0,T) VS T                  )
      CALL SPLOTS(7,NX,XP,T1SP)
      WRITE(NO,12)
12    FORMAT(1H ,//,31H T1(X,T) VS X, T = 0, 20,... 100)
      CALL SPLOTS(7,NX,XP,T2SP)
      WRITE(NO,13)
13    FORMAT(1H ,//,31H T2(X,T) VS X, T = 0, 20,... 100)
      CALL SPLOTS(7,NX,XP,T3SP)
      WRITE(NO,14)
14    FORMAT(1H ,//,31H T3(X,T) VS X, T = 0, 20,... 100)
      CALL SPLOTS(7,NX,XP,TSSP)
      WRITE(NO,15)
15    FORMAT(1H ,//,31H TS(X,T) VS X, T = 0, 20,... 100)
C...
C...  THE FOLLOWING CALLS TO SUBROUTINES BVODE1 AND BVODE2 (BOUNDARY-
C...  VALUE ORDINARY DIFFERENTIAL EQUATIONS, VERSIONS 1 AND 2) SOLVE
C...  THE STEADY STATE PROBLEM (T1  = T2  = T3  = TS  = 0) DIRECTLY, FOR
C...                              T     T     T     T
C...  COMPARISON WITH THE DYNAMIC SOLUTION AT ESSENTIALLY INFINITE TIME
C...  (STEADY STATE).  SUBROUTINES BVODE1 AND BVODE2 DIFFER IN THE WAY
C...  IN WHICH THE BOUNDARY CONDITIONS OF THE STEADY STATE BVODE PROBLEM
C...  ARE IMPLEMENTED.  THE COMMENTS IN BVODE1 AND BVODE2 GIVE ALL OF THE
C...  DETAILS
      CALL BVODE1
      CALL BVODE2
      RETURN
      END
      SUBROUTINE BVODE1
C...
C...  SUBROUTINE BVODE1 NUMERICALLY INTEGRATES THE BOUNDARY-VALUE
C...  ORDINARY DIFFERENTIAL EQUATIONS (BVODES) OBTAINED BY SETTING
C...  THE TIME DERIVATIVES TO ZERO IN THE THREE-PASS SHELL AND TUBE
C...  HEAT EXCHANGER PDES.  THE BVODES ARE EQUATIONS (1) TO (4)
C...  BELOW.  THE METHOD OF SOLUTION IS TO REPLACE THE SPATIAL
C...  DERIVATIVES (IN X) WITH THREE-POINT DIFFERENCE APPROXIMATIONS.
C...  THE RESULTING SET OF LINEAR ALGEBRAIC EQUATIONS IS THEN SOLVED
C...  BY CALLS TO SUBROUTINES DECOMP AND SOLVE.  THE DETAILS OF THE
C...  NUMERICAL SOLUTION ARE GIVEN IN THE FOLLOWING COMMENTS.
C...
C...  SUBROUTINE BVODE1 IS BASED ON DIRECT USE OF THE BOUNDARY CONDI-
C...  TIONS AS ALGEBRAIC EQUATIONS AT THE BOUNDARIES.  THIS CONTRASTS
C...  WITH SUBROUTINE BVODE2 IN WHICH EQUATIONS (1) TO (4) ARE APPROXI-
C...  MATED AT THE BOUNDARIES BY FINITE DIFFERENCES, THEN THE BOUNDARY
C...  CONDITIONS ARE SUBSTITUTED INTO THESE APPROXIMATIONS TO PRODUCE
C..   THE ALGEBRAIC EQUATIONS AT THE BOUNDARIES.
C...
C...     VT*DT1/DX = CT*(TS - T1)                                  (1)
C...
C...    -VT*DT2/DX = CT*(TS - T2)                                  (2)
C...
C...     VT*DT3/DX = CT*(TS - T3)                                  (3)
C...
C...    -VS*DTS/DX = CS*((T1 - TS) + (T2 - TS) + (T3 - TS))        (4)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  ACCESS THE INPUT/OUTPUT DEVICE NUMBERS
      COMMON/IO/NI,NO
C...
C...  DIMENSION THE ARRAYS FOR THE LINEAR ALGEBRAIC EQUATIONS
      DIMENSION A(44,44), B(44), IPVT(44), WORK(44)
C...
C...  ZERO THE COEFFICIENT MATRIX SO THAT THE SUBSEQUENT PROGRAMMING
C...  WILL DEFINE ONLY THE NONZERO ELEMENTS
      DO 13 I=1,44
      DO 13 J=1,44
      A(I,J)=0.0D0
13    CONTINUE
C...
C...  ******************************************************************
C...
C...  EQUATION (1)
C...
C...  GRID POINT 1 (ALGEBRAIC EQUATION (AE) 1)
C...
C...     T1(1) = 25
C...
C...  GRID POINT 2 (AE 2)
C...
C...     VT*(T1(3) - T1(1))/(2*DX) = CT*(TS(2) - T1(2))
C...
C...     -VT/(2*DX)*T1(1) + CT*T1(2) + VT/(2*DX)*T1(3) - CT*TS(2) = 0
C...
C...  GRID POINT 3 (AE 3)
C...
C...     VT*(T1(4) - T1(2))/(2*DX) = CT*(TS(3) - T1(3))
C...
C...     -VT/(2*DX)*T1(2) + CT*T1(3) + VT/(2*DX)*T1(4) - CT*TS(3) = 0
C...
C...  GRID POINT 4 (AE 4)
C...
C...     VT*(T1(5) - T1(3))/(2*DX) = CT*(TS(4) - T1(4))
C...
C...     -VT/(2*DX)*T1(3) + CT*T1(4) + VT/(2*DX)*T1(5) - CT*TS(4) = 0
C...        .                                                 .
C...        .                                                 .
C...        .                                                 .
C...
C...  GRID POINT 9 (AE 9)
C...
C...     VT*(T1(10) - T1(8))/(2*DX) = CT*(TS(9) - T1(9))
C...
C...     -VT/(2*DX)*T1(8) + CT*T1(9) + VT/(2*DX)*T1(10) - CT*TS(9) = 0
C...
C...  GRID POINT 10 (AE 10)
C...
C...     VT*(T1(11) - T1(9))/(2*DX) = CT*(TS(10) - T1(10))
C...
C...     -VT/(2*DX)*T1(9) + CT*T1(10) + VT/(2*DX)*T1(11) - CT*TS(10) = 0
C...
C...  GRID POINT 11 (AE 11)
C...
C...     VT*(3*T1(11) - 4*T1(10) + T1(9))/(2*DX) = CT*(TS(11) - T1(11))
C...
C...    VT/(2*DX)*T1(9) - 4*VT/(2*DX)*T1(10) + (3*VT/(2*DX) + CT)*T1(11)
C...
C...                     - CT*TS(11) = 0
C...
C...  SET THE MODEL PARAMETERS
         DX=10.0D0/DFLOAT(11-1)
         VT=1.0D0
         VS=0.5D0
         C1=VT/(2.0D0*DX)
         C2=VS/(2.0D0*DX)
         CT=0.1D0
         CS=0.1D0
C...
C...  ******************************************************************
C...
C...  DEFINE THE COEFFICIENT MATRIX FOR THE FINITE DIFFERENCE APPROXI-
C...  MATION OF EQUATION (1)
C...
C...  AE 1
         A(1,1)=1.0D0
C...
C...  AES 2 TO 10
         DO 1 I=2,10
         A(I,I-1)=-C1
         A(I,I  )= CT
         A(I,I+1)= C1
         IS=3*11+I
         A(I,IS )=-CT
1        CONTINUE
C...
C...  AE 11
         A(11, 9)=C1
         A(11,10)=-4.0D0*C1
         A(11,11)= 3.0D0*C1+CT
         A(11,44)=-CT
C...
C...  ******************************************************************
C...
C...  EQUATION (2)
C...
C...  GRID POINT 11 (AE 22)
C...
C...     T2(11) = T1(11)
C...
C...  GRID POINT 10 (AE 21)
C...
C...     -VT*(T2(11) - T2(9))/(2*DX) = CT*(TS(10) - T2(10))
C...
C...     VT/(2*DX)*T2(9) + CT*T2(10) - VT/(2*DX)*T2(11) - CT*TS(10) = 0
C...        .                                                 .
C...        .                                                 .
C...        .                                                 .
C...
C...  GRID POINT 2 (AE 13)
C...
C...     -VT*(T2(3) - T2(1))/(2*DX) = CT*(TS(2) - T2(2))
C...
C...     VT/(2*DX)*T2(1) + CT*T2(2) - VT/(2*DX)*T2(3) - CT*TS(2) = 0
C...
C...  GRID POINT 1 (AE 12)
C...
C...     -VT*(-T2(3) + 4*T2(2) - 3*T2(1))/(2*DX) = CT*(TS(1) - T2(1))
C...
C...     (3*VT/(2*DX)+CT)*T2(1) - 4*VT/(2*DX)*T2(2) + VT/(2*DX)*T2(3)
C...
C...                         - CT*TS(1) = 0
C...
C...  AE 12
         A(12,12)= 3.0D0*C1+CT
         A(12,13)=-4.0D0*C1
         A(12,14)= C1
         A(12,34)=-CT
C...
C...  AES 13 TO 21
         DO 2 I=13,21
         A(I,I-1)= C1
         A(I,I  )= CT
         A(I,I+1)=-C1
         IS=2*11+I
         A(I,IS )=-CT
2        CONTINUE
C...
C...  AE 22
         A(22,11)=-1.0D0
         A(22,22)= 1.0D0
C...
C...  ******************************************************************
C...
C...  EQUATION (3)
C...
C...  GRID POINT 1 (AE 23)
C...
C...     T3(1) = T2(1)
C...
C...  GRID POINT 2 (AE 24)
C...
C...     VT*(T3(3) - T3(1))/(2*DX) = CT*(TS(2) - T3(2))
C...
C...     -VT/(2*DX)*T3(1) + CT*T3(2) + VT/(2*DX)*T3(3) - CT*TS(2) = 0
C...        .                                                 .
C...        .                                                 .
C...        .                                                 .
C...
C...  GRID POINT 10 (AE 32)
C...
C...     VT*(T3(11) - T3(9))/(2*DX) = CT*(TS(10) - T3(10))
C...
C...    -VT/(2*DX)*(T3(9) + CT*T3(10) + VT*T3(11)/(2*DX) - CT*TS(10) = 0
C...
C...  GRID POINT 11 (AE 33)
C...
C...     VT*(3*T1(11) - 4*T3(10) + T3(9))/(2*DX) = CT*(TS(11) - T3(11))
C...
C...     VT/(2*DX)*T3(9) - 4*VT/(2*DX)*T3(10) + (3*VT/(2*DX)
C...
C...                + CT*T3(11) - CT*TS(11) = 0
C...
C...  AE 23
         A(23,12)=-1.0D0
         A(23,23)= 1.0D0
C...
C...  AES 24 TO 32
         DO 3 I=24,32
         A(I,I-1)=-C1
         A(I,I  )= CT
         A(I,I+1)= C1
         IS=11+I
         A(I,IS )=-CT
3        CONTINUE
C...
C...  AE 33
         A(33,31)=C1
         A(33,32)=-4.0D0*C1
         A(33,33)= 3.0D0*C1+CT
         A(33,44)=-CT
C...
C...  ******************************************************************
C...
C...  EQUATION (4)
C...
C...  GRID POINT 11 (AE 44)
C...
C...     TS(11) = 220
C...
C...  GRID POINT 10 (AE 43)
C...
C...     -VS*(TS(11) - TS(9))/(2*DX) = CS*((T1(10) - TS(10))
C...
C...          + (T2(10) - TS(10)) + (T3(10) - TS(10)))
C...
C...     VS/(2*DX)*TS(9) + 3*CS*TS(10) - VS/(2*DX)*TS(11) - CS*T1(10)
C...
C...                     - CS*T2(10) - CS*T3(10) = 0
C...        .                                                 .
C...        .                                                 .
C...        .                                                 .
C...
C...  GRID POINT 2 (AE 35)
C...
C...     -VS*(TS(3) - TS(1))/(2*DX) = CS*((T1(2) - TS(2))
C...
C...          + (T2(2) - TS(2)) + (T3(2) - TS(2)))
C...
C...     VS/(2*DX)*TS(1) + 3*CS*TS(2) - VS/(2*DX)*TS(3) - CS*T1(2)
C...
C...                   - CS*T2(2) - CS*T3(2) = 0
C...
C...  GRID POINT 1 (AE 34)
C...
C...     -VS*(-TS(3) + 4*TS(2) - 3*TS(1))/(2*DX) = CS*(T1(1) - TS(1))
C...
C...               + (T2(1) - TS(1)) + (T3(1) - TS(1)))
C...
C...    (3*VS/(2*DX) + 3*CS)*TS(1) - 4*VS/(2*DX)*TS(2) + VS/(2*DX)*TS(3)
C...
C...                 - CS*T1(1) - CS*T2(1) - CS*T3(1) = 0
C...
C...  AE 34
         A(34,34)= 3.0D0*C2+3.0D0*CS
         A(34,35)=-4.0D0*C2
         A(34,36)= C2
         A(34, 1)=-CS
         A(34,12)=-CS
         A(34,23)=-CS
C...
C...  AES 35 TO 43
         DO 4 I=35,43
         A(I,I-1)= C2
         A(I,I  )=3.0D0*CS
         A(I,I+1)=-C2
         IS=I-3*11
         A(I,IS )=-CS
         IS=I-2*11
         A(I,IS )=-CS
         IS=I-1*11
         A(I,IS )=-CS
4        CONTINUE
C...
C...  AE 44
         A(44,44)=1.0D0
C...
C...  ******************************************************************
C...
C...  RIGHT HAND SIDE (RHS) VECTOR
C...
C...  AE 1
         B(1)=25.0D0
C...
C...  AE 2 TO 43
         DO 9 I=2,43
         B(I)=0.0D0
9        CONTINUE
C...
C...  AE 44
         B(44)=220.0D0
C...
C...  ******************************************************************
C...
C...  MAP THE COEFFICIENT MATRIX, A, TO DISPLAY ITS STRUCTURE
      N=44
      CALL MAP(N,A)
C...
C...  PRINT THE COEFFICIENT MATRIX TO CONFIRM THE VALUES OF THE
C...  INDIVIDUAL ELEMENTS
      CALL CHECK(N,A,IPVT)
C...
C...  SOLVE THE LINEAR ALGEBRAIC SYSTEM
      CALL DECOMP(N,N,A,COND,IPVT,WORK)
C...
C...  CHECK FOR A SINGULAR OR ILL-CONDITIONED SYSTEM
      IF(COND.LT.1.0D+20)GO TO 10
C...
C...  THE CONDITION NUMBER OF THE COEFFICIENT MATRIX, A, EXCEEDS 1.0D+20
C...  SO A WARNING MESSAGE IS PRINTED AND EXECUTION IS TERMINATED
      WRITE(NO,11)COND
11    FORMAT(26H CONDITION NUMBER OF A IS ,E12.4)
      STOP
C...
C...  COMPLETE AND PRINT THE SOLUTION
10    CALL SOLVE(N,N,A,B,IPVT)
      WRITE(NO,5)(B(I),I=1,11)
5     FORMAT(1H1,20H SOLUTION TO BVODE 1,/,10X,5F10.3,/,6F10.3,///)
      WRITE(NO,6)(B(I),I=12,22)
6     FORMAT(1H ,20H SOLUTION TO BVODE 2,/,10X,5F10.3,/,6F10.3,///)
      WRITE(NO,7)(B(I),I=23,33)
7     FORMAT(1H ,20H SOLUTION TO BVODE 3,/,10X,5F10.3,/,6F10.3,///)
      WRITE(NO,8)(B(I),I=34,44)
8     FORMAT(1H ,20H SOLUTION TO BVODE 4,/,10X,5F10.3,/,6F10.3,///)
      RETURN
      END
      SUBROUTINE BVODE2
C...
C...  SUBROUTINE BVODE2 NUMERICALLY INTEGRATES THE BOUNDARY-VALUE
C...  ORDINARY DIFFERENTIAL EQUATIONS (BVODES) OBTAINED BY SETTING
C...  THE TIME DERIVATIVES TO ZERO IN THE THREE-PASS SHELL AND TUBE
C...  HEAT EXCHANGER PDES.  THE BVODES ARE EQUATIONS (1) TO (4)
C...  BELOW.  THE METHOD OF SOLUTION IS TO REPLACE THE SPATIAL
C...  DERIVATIVES (IN X) WITH THREE-POINT DIFFERENCE APPROXIMATIONS.
C...  THE RESULTING SET OF LINEAR ALGEBRAIC EQUATIONS IS THEN SOLVED
C...  BY CALLS TO SUBROUTINES DECOMP AND SOLVE.  THE DETAILS OF THE
C...  NUMERICAL SOLUTION ARE GIVEN IN THE FOLLOWING COMMENTS.
C...
C...  SUBROUTINE BVODE2 IS BASED ON FINITE DIFFERENCE APPROXIMATIONS
C...  OF EQUATIONS (1) TO (4) AT THE BOUNDARIES WITH THE FOUR BOUNDARY
C...  CONDITIONS SUBSTITUTED INTO THESE APPROXIMATIONS TO PRODUCE THE
C...  ALGEBRAIC EQUATIONS AT THE BOUNDARIES.  THIS CONTRASTS WITH SUB-
C...  ROUTINE BVODE1 IN WHICH THE BOUNDARY CONDITIONS ARE USED DIRECTLY
C...  AS ALGEBRAIC EQUATIONS.
C...
C...     VT*DT1/DX = CT*(TS - T1)                                  (1)
C...
C...    -VT*DT2/DX = CT*(TS - T2)                                  (2)
C...
C...     VT*DT3/DX = CT*(TS - T3)                                  (3)
C...
C...    -VS*DTS/DX = CS*((T1 - TS) + (T2 - TS) + (T3 - TS))        (4)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  ACCESS THE INPUT/OUTPUT DEVICE NUMBERS
      COMMON/IO/NI,NO
C...
C...  DIMENSION THE ARRAYS FOR THE LINEAR ALGEBRAIC EQUATIONS
      DIMENSION A(44,44), B(44), IPVT(44), WORK(44)
C...
C...  ZERO THE COEFFICIENT MATRIX SO THAT THE SUBSEQUENT PROGRAMMING
C...  WILL DEFINE ONLY THE NONZERO ELEMENTS
      DO 13 I=1,44
      DO 13 J=1,44
      A(I,J)=0.0D0
13    CONTINUE
C...
C...  ******************************************************************
C...
C...  EQUATION (1)
C...
C...  GRID POINT 1 (ALGEBRAIC EQUATION (AE) 1)
C...
C...     VT*(-T1(3) + 4*T1(2) - 3*T1(1))/(2*DX) = CT*(TS(1) - T1(1))
C...
C...     4*VT/(2*DX)*T1(2) - VT/(2*DX)*T1(3) - CT*TS(1) =
C...
C...                  (3*VT/(2*DX) - CT)*25
C...
C...  GRID POINT 2 (AE 2)
C...
C...     VT*(T1(3) - T1(1))/(2*DX) = CT*(TS(2) - T1(2))
C...
C...     -VT/(2*DX)*T1(1) + CT*T1(2) + VT/(2*DX)*T1(3) - CT*TS(2) = 0
C...
C...  GRID POINT 3 (AE 3)
C...
C...     VT*(T1(4) - T1(2))/(2*DX) = CT*(TS(3) - T1(3))
C...
C...     -VT/(2*DX)*T1(2) + CT*T1(3) + VT/(2*DX)*T1(4) - CT*TS(3) = 0
C...
C...  GRID POINT 4 (AE 4)
C...
C...     VT*(T1(5) - T1(3))/(2*DX) = CT*(TS(4) - T1(4))
C...
C...     -VT/(2*DX)*T1(3) + CT*T1(4) + VT/(2*DX)*T1(5) - CT*TS(4) = 0
C...        .                                                 .
C...        .                                                 .
C...        .                                                 .
C...
C...  GRID POINT 9 (AE 9)
C...
C...     VT*(T1(10) - T1(8))/(2*DX) = CT*(TS(9) - T1(9))
C...
C...     -VT/(2*DX)*T1(8) + CT*T1(9) + VT/(2*DX)*T1(10) - CT*TS(9) = 0
C...
C...  GRID POINT 10 (AE 10)
C...
C...     VT*(T1(11) - T1(9))/(2*DX) = CT*(TS(10) - T1(10))
C...
C...     -VT/(2*DX)*T1(9) + CT*T1(10) + VT/(2*DX)*T1(11) - CT*TS(10) = 0
C...
C...  GRID POINT 11 (AE 11)
C...
C...     VT*(3*T1(11) - 4*T1(10) + T1(9))/(2*DX) = CT*(TS(11) - T1(11))
C...
C...    VT/(2*DX)*T1(9) - 4*VT/(2*DX)*T1(10) + (3*VT/(2*DX) + CT)*T1(11)
C...
C...                     - CT*TS(11) = 0
C...
C...  SET THE MODEL PARAMETERS
         DX=10.0D0/DFLOAT(11-1)
         VT=1.0D0
         VS=0.5D0
         C1=VT/(2.0D0*DX)
         C2=VS/(2.0D0*DX)
         CT=0.10D0
         CS=0.10D0
C...
C...  ******************************************************************
C...
C...  DEFINE THE COEFFICIENT MATRIX FOR THE FINITE DIFFERENCE APPROXI-
C...  MATION OF EQUATION (1)
C...
C...  AE 1
         A(1,1)= 0.0D0
         A(1,2)= 4.0D0*C1
         A(1,3)= -C1
         A(1,34)=-CT
C...
C...  AES 2 TO 10
         DO 1 I=2,10
         A(I,I-1)=-C1
         A(I,I  )= CT
         A(I,I+1)= C1
         IS=3*11+I
         A(I,IS )=-CT
1        CONTINUE
C...
C...  AE 11
         A(11, 9)= C1
         A(11,10)=-4.0D0*C1
         A(11,11)= 3.0D0*C1+CT
         A(11,44)=-CT
C...
C...  ******************************************************************
C...
C...  EQUATION (2)
C...
C...  GRID POINT 11 (AE 22)
C...
C...     -VT*(3*T2(11) - 4*T2(10) + T2(9))/(2*DX) = CT*(TS(11) - T2(11))
C...
C...    -VT/(2*DX)*T2(9) + 4*VT/(2*DX)*T2(10) + (-3*VT/(2*DX)+CT)*T1(11)
C...
C...                            -CT*TS(11) = 0
C...
C...  GRID POINT 10 (AE 21)
C...
C...     -VT*(T2(11) - T2(9))/(2*DX) = CT*(TS(10) - T2(10))
C...
C...     VT/(2*DX)*T2(9) + CT*T2(10) - VT/(2*DX)*T2(11) - CT*TS(10) = 0
C...        .                                                 .
C...        .                                                 .
C...        .                                                 .
C...
C...  GRID POINT 2 (AE 13)
C...
C...     -VT*(T2(3) - T2(1))/(2*DX) = CT*(TS(2) - T2(2))
C...
C...     VT/(2*DX)*T2(1) + CT*T2(2) - VT/(2*DX)*T2(3) - CT*TS(2) = 0
C...
C...  GRID POINT 1 (AE 12)
C...
C...     -VT*(-T2(3) + 4*T2(2) - 3*T2(1))/(2*DX) = CT*(TS(1) - T2(1))
C...
C...     (3*VT/(2*DX)+CT)*T2(1) - 4*VT/(2*DX)*T2(2) + VT/(2*DX)*T2(3)
C...
C...                         - CT*TS(1) = 0
C...
C...  AE 12
         A(12,12)= 3.0D0*C1+CT
         A(12,13)=-4.0D0*C1
         A(12,14)= C1
         A(12,34)=-CT
C...
C...  AES 13 TO 21
         DO 2 I=13,21
         A(I,I-1)= C1
         A(I,I  )= CT
         A(I,I+1)=-C1
         IS=2*11+I
         A(I,IS )=-CT
2        CONTINUE
C...
C...  AE 22
         A(22,20)=-C1
         A(22,21)= 4.0D0*C1
         A(22,11)=-3.0D0*C1+CT
         A(22,44)=-CT
C...
C...  ******************************************************************
C...
C...  EQUATION (3)
C...
C...  GRID POINT 1 (AE 23)
C...
C...     VT*(-T3(3) + 4*T3(2) - 3*T3(1))/(2*DX) = CT*(TS(1) - T3(1))
C...
C...     4*VT/(2*DX)*T3(2) - VT/(2*DX)*T3(3) + (-3*VT/(2*DX)+CT)*T2(1)
C...
C...                         - CT*TS(1) = 0
C...
C...  GRID POINT 2 (AE 24)
C...
C...     VT*(T3(3) - T3(1))/(2*DX) = CT*(TS(2) - T3(2))
C...
C...     -VT/(2*DX)*T3(1) + CT*T3(2) + VT/(2*DX)*T3(3) - CT*TS(2) = 0
C...        .                                                 .
C...        .                                                 .
C...        .                                                 .
C...
C...  GRID POINT 10 (AE 32)
C...
C...     VT*(T3(11) - T3(9))/(2*DX) = CT*(TS(10) - T3(10))
C...
C...    -VT/(2*DX)*(T3(9) + CT*T3(10) + VT*T3(11)/(2*DX) - CT*TS(10) = 0
C...
C...  GRID POINT 11 (AE 33)
C...
C...     VT*(3*T1(11) - 4*T3(10) + T3(9))/(2*DX) = CT*(TS(11) - T3(11))
C...
C...     VT/(2*DX)*T3(9) - 4*VT/(2*DX)*T3(10) + (3*VT/(2*DX)
C...
C...                + CT*T3(11) - CT*TS(11) = 0
C...
C...  AE 23
         A(23,23)= 0.0D0
         A(23,24)= 4.0D0*C1
         A(23,25)=-C1
         A(23,12)=-3.0D0*C1+CT
         A(23,34)=-CT
C...
C...  AES 24 TO 32
         DO 3 I=24,32
         A(I,I-1)=-C1
         A(I,I  )= CT
         A(I,I+1)= C1
         IS=11+I
         A(I,IS )=-CT
3        CONTINUE
C...
C...  AE 33
         A(33,31)= C1
         A(33,32)=-4.0D0*C1
         A(33,33)= 3.0D0*C1+CT
         A(33,44)=-CT
C...
C...  ******************************************************************
C...
C...  EQUATION (4)
C...
C...  GRID POINT 11 (AE 44)
C...
C...     -VS*(3*TS(11) - 4*TS(10) + TS(9))/(2*DX) = CS*(T1(11) - TS(11))
C...
C...              + (T2(11) - TS(11)) + (T3(11) - TS(11)))
C...
C...     -VS/(2*DX)*TS(9) + 4*VS/(2*DX)*TS(10) - CS*(T1(11) + T2(11)
C...
C...              + T3(11)) = (3*VS/(2*DX) - 3*CS)*220
C...
C...  GRID POINT 10 (AE 43)
C...
C...     -VS*(TS(11) - TS(9))/(2*DX) = CS*((T1(10) - TS(10))
C...
C...          + (T2(10) - TS(10)) + (T3(10) - TS(10)))
C...
C...     VS/(2*DX)*TS(9) + 3*CS*TS(10) - VS/(2*DX)*TS(11) - CS*T1(10)
C...
C...                     - CS*T2(10) - CS*T3(10) = 0
C...        .                                                 .
C...        .                                                 .
C...        .                                                 .
C...
C...  GRID POINT 2 (AE 35)
C...
C...     -VS*(TS(3) - TS(1))/(2*DX) = CS*((T1(2) - TS(2))
C...
C...          + (T2(2) - TS(2)) + (T3(2) - TS(2)))
C...
C...     VS/(2*DX)*TS(1) + 3*CS*TS(2) - VS/(2*DX)*TS(3) - CS*T1(2)
C...
C...                   - CS*T2(2) - CS*T3(2) = 0
C...
C...  GRID POINT 1 (AE 34)
C...
C...     -VS*(-TS(3) + 4*TS(2) - 3*TS(1))/(2*DX) = CS*(T1(1) - TS(1))
C...
C...               + (T2(1) - TS(1)) + (T3(1) - TS(1)))
C...
C...    (3*VS/(2*DX) + 3*CS)*TS(1) - 4*VS/(2*DX)*TS(2) + VS/(2*DX)*TS(3)
C...
C...                 - CS*T1(1) - CS*T2(1) - CS*T3(1) = 0
C...
C...  AE 34
         A(34,34)= 3.0D0*C2+3.0D0*CS
         A(34,35)=-4.0D0*C2
         A(34,36)= C2
         A(34, 1)=-CS
         A(34,12)=-CS
         A(34,23)=-CS
C...
C...  AES 35 TO 43
         DO 4 I=35,43
         A(I,I-1)= C2
         A(I,I  )= 3.0D0*CS
         A(I,I+1)=-C2
         IS=I-3*11
         A(I,IS )=-CS
         IS=I-2*11
         A(I,IS )=-CS
         IS=I-1*11
         A(I,IS )=-CS
4        CONTINUE
C...
C...  AE 44
         A(44,44)=0.0D0
         A(44,43)=4.0D0*C2
         A(44,42)=-C2
         A(44,11)=-CS
         A(44,22)=-CS
         A(44,33)=-CS
C...
C...  ******************************************************************
C...
C...  RIGHT HAND SIDE (RHS) VECTOR
C...
C...  AE 1
         B(1)=(3.0D0*C1-CT)*25.0D0
C...
C...  AE 2 TO 43
         DO 9 I=2,43
         B(I)=0.0D0
9        CONTINUE
C...
C...  AE 44
         B(44)=(3.0D0*C2-3.0D0*CS)*220.0D0
C...
C...  ******************************************************************
C...
C...  MAP THE COEFFICIENT MATRIX, A, TO DISPLAY ITS STRUCTURE
      N=44
      CALL MAP(N,A)
C...
C...  PRINT THE COEFFICIENT MATRIX TO CONFIRM THE VALUES OF THE
C...  INDIVIDUAL ELEMENTS
      CALL CHECK(N,A,IPVT)
C...
C...  SOLVE THE LINEAR ALGEBRAIC SYSTEM
      CALL DECOMP(N,N,A,COND,IPVT,WORK)
C...
C...  CHECK FOR A SINGULAR OR ILL-CONDITIONED SYSTEM
      IF(COND.LT.1.0D+20)GO TO 10
C...
C...  THE CONDITION NUMBER OF THE COEFFICIENT MATRIX, A, EXCEEDS 1.0D+20
C...  SO A WARNING MESSAGE IS PRINTED AND EXECUTION IS TERMINATED
      WRITE(NO,11)COND
11    FORMAT(26H CONDITION NUMBER OF A IS ,E12.4)
      STOP
C...
C...  COMPLETE AND PRINT THE SOLUTION
10    CALL SOLVE(N,N,A,B,IPVT)
      WRITE(NO,5)(B(I),I=1,11)
5     FORMAT(1H1,20H SOLUTION TO BVODE 1,/,10X,5F10.3,/,6F10.3,///)
      WRITE(NO,6)(B(I),I=12,22)
6     FORMAT(1H ,20H SOLUTION TO BVODE 2,/,10X,5F10.3,/,6F10.3,///)
      WRITE(NO,7)(B(I),I=23,33)
7     FORMAT(1H ,20H SOLUTION TO BVODE 3,/,10X,5F10.3,/,6F10.3,///)
      WRITE(NO,8)(B(I),I=34,44)
8     FORMAT(1H ,20H SOLUTION TO BVODE 4,/,10X,5F10.3,/,6F10.3,///)
      RETURN
      END
      SUBROUTINE MAP(N,A)
C...
C...  SUBROUTINE MAP IS A BASIC MAPPING ROUTINE FOR LINEAR ALGEBRAIC
C...  EQUATIONS.  IT CAN ALSO MAP THE JACOBIAN MATRIX OF A SYSTEM OF
C...  NONLINEAR ALGEBRAIC OR ORDINARY DIFFERENTIAL EQUATIONS.  THE
C...  SIZE OF THE MATRIX TO BE MAPPED IS LIMITED TO 49 X 49 (THIS IS
C...  AN ARBITRARY CHOICE, AND THE ROUTINE CAN EASILY BE MODIFIED TO
C...  ACCOMMODATE LARGER MATRICES).  FOR A MORE SOPHISTICATED MAPPING
C...  ROUTINE, THE READER IS REFERRED TO SUBROUTINE JMAP AVAILABLE IN
C...  DSS/2.
C...
C...     N      ORDER (ROW AND COLUMN DIMENSIONS) OF THE COEFFICIENT
C...            MATRIX TO BE MAPPED (INPUT)
C...
C...     A      TWO DIMENSIONAL ARRAY TO STORE THE COEFFICIENT MATRIX
C...            TO BE MAPPED (INPUT)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  VARIABLE DIMENSION THE ARRAY FOR THE COEFFICIENT MATRIX
      DIMENSION A(N,N)
C...
C...  ABSOLUTE DIMENSION THE ARRAYS FOR THE SYMBOLS, FIRST AND SECOND
C...  LINES OF THE MAP
      DIMENSION SYMBOL(11), HLINE1(51), HLINE2(50)
C...
C...  ACCESS THE INPUT/OUTPUT DEVICE NUMBERS
      COMMON/IO/NI,NO
C...
C...  DEFINE THE SYMBOLS PRINTED IN THE MAP
      DATA SYMBOL/2H 1,2H 2,2H 3,2H 4,2H 5,2H 6,2H 7,2H 8,2H 9,
     1            2H 0,2H  /,XCHAR/2H X/
C...
C...  DEFINE THE SYMBOLS PRINTED IN THE FIRST (TENS) ROW OF THE MAP
      DO 1 I=1,5
      DO 1 J=1,10
      JS=10*(I-1)+J
      IF(I.EQ.1)HLINE1(JS)=SYMBOL(11)
      IF(I.EQ.2)HLINE1(JS)=SYMBOL(1)
      IF(I.EQ.3)HLINE1(JS)=SYMBOL(2)
      IF(I.EQ.4)HLINE1(JS)=SYMBOL(3)
      IF(I.EQ.5)HLINE1(JS)=SYMBOL(4)
1     CONTINUE
C...
C...  DEFINE THE SYMBOLS PRINTED IN THE SECOND (UNITS) ROW OF THE MAP
      DO 2 I=1,5
      DO 2 J=1,10
      JS=10*(I-1)+J
      HLINE2(JS)=SYMBOL(J)
2     CONTINUE
C...
C...  PRINT A LABEL FOR THE HORIZONTAL AND VERTICAL LINES OF THE
C...  MAP
      WRITE(NO,8)
8     FORMAT(1H1,15X,40H COLUMN INDEX IS FOR THE UNKNOWNS -----+ ,//,
     1 1H ,4X,1H.,/,1H ,4X,1H.,/,
     2 1H ,4X,1H.,   40H ROW INDEX IS FOR THE EQUATIONS          ,/,
     3 1H ,4X,1H.,/,1H ,4X,1H.,/,
     4 1H ,4X,1H+,/)
C...
C...  PRINT THE FIRST (TENS) AND SECOND (UNITS) ROWS OF THE MAP
      NP1=N+1
      WRITE(NO,3)(HLINE1(I),I=1,NP1)
3     FORMAT(1H , 8X,50A2)
      WRITE(NO,4)(HLINE2(I),I=1,N)
4     FORMAT(1H ,10X,50A2)
C...
C...  PRINT THE ITH ROW OF THE MAP, 1 LE I LE N
      DO 6 I=1,N
C...
C...  INITIALLY, BLANK ALL J CHARACTERS OF THE ITH ROW OF THE MAP
      DO 5 J=1,N
      HLINE1(J)=SYMBOL(11)
C...
C...  IF THE JTH COLUMN ELEMENT OF THE COEFFICIENT MATRIX, A, IN THE
C...  ITH ROW EXCEEDS 1.0D-10 (ARBITRARILY SELECTED), SET THE CHARACTER
C...  TO X TO BE PRINTED IN THE I,J POSITION OF THE MAP
      IF(ABS(A(I,J)).GT.1.0D-10)HLINE1(J)=XCHAR
5     CONTINUE
C...
C...  PRINT THE ITH ROW OF THE MAP
      WRITE(NO,7)I,(HLINE1(J),J=1,N)
7     FORMAT(1H ,/,4X,I2,6X,50A2)
C...
C...  PRINT ALL N ROWS OF THE MAP BY GOING THROUGH DO LOOP 6 N TIMES
6     CONTINUE
      RETURN
      END
      SUBROUTINE CHECK(N,A,NEQ)
C...
C...  SUBROUTINE CHECK PRINTS OUT THE NUMERICAL VALUES OF THE ELEMENTS
C...  OF THE COEFFICIENT MATRIX, A.  IT IS HIGHLY SPECIFIC TO THIS
C...  PARTICULAR PROBLEM, E.G., DO LOOPS HAVE NUMERICAL UPPER LIMITS
C...  CORRESPONDING TO N = 44.
C...
C...  ARGUMENT LIST
C...
C...     N      ORDER (ROW AND COLUMN DIMENSIONS) OF THE COEFFICIENT
C...            MATRIX TO BE PRINTED (INPUT)
C...
C...     A      TWO DIMENSIONAL ARRAY CONTAINING THE COEFFICIENT MATRIX
C...            TO BE PRINTED (INPUT)
C...
C...     NEQ  CONE DIMENSIONAL, INTEGER WORK ARRAY (INPUT)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  VARIABLE DIMENSION THE ARRAYS USED BY SUBROUTINE CHECK
      DIMENSION A(N,N), NEQ(N)
C...
C...  ACCESS THE INPUT/OUTPUT DEVICE NUMBERS
      COMMON/IO/NI,NO
C...
C...  STORE THE COLUMN SUBSCRIPTS OF THE COEFFICIENT MATRIX FOR
C...  SUBSEQUENT PRINTING
      DO 1 I=1,44
      NEQ(I)=I
1     CONTINUE
C...
C...  PRINT THE ALGEBRAIC APPROXIMATIONS OF EACH OF THE BOUNDARY-VALUE
C...  ODES (BVODES)
      DO 6 NEQN=1,4
      WRITE(NO,2)NEQN
2     FORMAT(1H1,14H ODE EQUATION ,I2,/)
C...
C...  STEP THROUGH THE 11 ALBEGRAIC APPROXIMATIONS FOR EACH BVODE
      DO 3 I=1,11
      IAE=11*(NEQN-1)+I
C...
C...  PRINT THE ITH ROW OF THE COEFFICIENT MATRIX FOR BVODE NEQN
      WRITE(NO,4)IAE,(NEQ(J),J=1,22),(A(I,J),J=1,22)
4     FORMAT(4H AE ,I2,2X,22I5,/,8X,22F5.1  )
      WRITE(NO,5)(NEQ(J),J=23,44),(A(I,J),J=23,44)
5     FORMAT(          8X,22I5,/,8X,22F5.1,/)
C...
C...  GO TO THE NEXT ROW OF THE COEFFICIENT MATRIX VIA DO LOOP 3
3     CONTINUE
C...
C...  GO TO THE NEXT BVODE VIA DO LOOP 6
6     CONTINUE
      RETURN
      END
*DECK DECOMP
      SUBROUTINE DECOMP(NDIM,N,A,COND,IPVT,WORK)
C
      INTEGER NDIM,N
      DOUBLE PRECISION A(NDIM,N),COND,WORK(N)
      INTEGER IPVT(N)
C
C     DECOMPOSES A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION
C     AND ESTIMATES THE CONDITION OF THE MATRIX.
C
C     USE SOLVE TO COMPUTE SOLUTIONS TO LINEAR SYSTEMS.
C
C     INPUT..
C
C        NDIM = DECLARED ROW DIMENSION OF THE ARRAY CONTAINING  A.
C
C        N = ORDER OF THE MATRIX.
C
C        A = MATRIX TO BE TRIANGULARIZED.
C
C     OUTPUT..
C
C        A  CONTAINS AN UPPER TRIANGULAR MATRIX  U  AND A PERMUTED
C          VERSION OF A LOWER TRIANGULAR MATRIX  I-L  SO THAT
C          (PERMUTATION MATRIX)*A = L*U .
C
C        COND = AN ESTIMATE OF THE CONDITION OF  A .
C           FOR THE LINEAR SYSTEM  A*X = B, CHANGES IN  A  AND  B
C           MAY CAUSE CHANGES  COND  TIMES AS LARGE IN  X .
C           IF  COND+1.0 .EQ. COND , A IS SINGULAR TO WORKING
C           PRECISION.  COND  IS SET TO  1.0D+32  IF EXACT
C           SINGULARITY IS DETECTED.
C
C        IPVT = THE PIVOT VECTOR.
C           IPVT(K) = THE INDEX OF THE K-TH PIVOT ROW
C           IPVT(N) = (-1)**(NUMBER OF INTERCHANGES)
C
C     WORK SPACE..  THE VECTOR  WORK  MUST BE DECLARED AND INCLUDED
C                   IN THE CALL.  ITS INPUT CONTENTS ARE IGNORED.
C                   ITS OUTPUT CONTENTS ARE USUALLY UNIMPORTANT.
C
C     THE DETERMINANT OF A CAN BE OBTAINED ON OUTPUT BY
C        DET(A) = IPVT(N) * A(1,1) * A(2,2) * ... * A(N,N).
C
      DOUBLE PRECISION EK, T, ANORM, YNORM, ZNORM
      INTEGER NM1, I, J, K, KP1, KB, KM1, M
      DOUBLE PRECISION DABS, DSIGN
C
      IPVT(N) = 1
      IF (N .EQ. 1) GO TO 80
      NM1 = N - 1
C
C     COMPUTE 1-NORM OF A
C
      ANORM = 0.0D0
      DO 10 J = 1, N
         T = 0.0D0
         DO 5 I = 1, N
            T = T + DABS(A(I,J))
    5    CONTINUE
         IF (T .GT. ANORM) ANORM = T
   10 CONTINUE
C
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
      DO 35 K = 1,NM1
         KP1= K+1
C
C        FIND PIVOT
C
         M = K
         DO 15 I = KP1,N
            IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
   15    CONTINUE
         IPVT(K) = M
         IF (M .NE. K) IPVT(N) = -IPVT(N)
         T = A(M,K)
         A(M,K) = A(K,K)
         A(K,K) = T
C
C        SKIP STEP IF PIVOT IS ZERO
C
         IF (T .EQ. 0.0D0) GO TO 35
C
C        COMPUTE MULTIPLIERS
C
         DO 20 I = KP1,N
             A(I,K) = -A(I,K)/T
   20    CONTINUE
C
C        INTERCHANGE AND ELIMINATE BY COLUMNS
C
         DO 30 J = KP1,N
             T = A(M,J)
             A(M,J) = A(K,J)
             A(K,J) = T
             IF (T .EQ. 0.0D0) GO TO 30
             DO 25 I = KP1,N
                A(I,J) = A(I,J) + A(I,K)*T
   25        CONTINUE
   30    CONTINUE
   35 CONTINUE
C
C     COND = (1-NORM OF A)*(AN ESTIMATE OF 1-NORM OF A-INVERSE)
C     ESTIMATE OBTAINED BY ONE STEP OF INVERSE ITERATION FOR THE
C     SMALL SINGULAR VECTOR.  THIS INVOLVES SOLVING TWO SYSTEMS
C     OF EQUATIONS, (A-TRANSPOSE)*Y = E  AND  A*Z = Y  WHERE  E
C     IS A VECTOR OF +1 OR -1 CHOSEN TO CAUSE GROWTH IN Y.
C     ESTIMATE = (1-NORM OF Z)/(1-NORM OF Y)
C
C     SOLVE (A-TRANSPOSE)*Y = E
C
      DO 50 K = 1, N
         T = 0.0D0
         IF (K .EQ. 1) GO TO 45
         KM1 = K-1
         DO 40 I = 1, KM1
            T = T + A(I,K)*WORK(I)
   40    CONTINUE
   45    EK = 1.0D0
         IF (T .LT. 0.0D0) EK = -1.0D0
         IF (A(K,K) .EQ. 0.0D0) GO TO 90
         WORK(K) = -(EK + T)/A(K,K)
   50 CONTINUE
      DO 60 KB = 1, NM1
         K = N - KB
         T = 0.0D0
         KP1 = K+1
         DO 55 I = KP1, N
            T = T + A(I,K)*WORK(K)
   55    CONTINUE
         WORK(K) = T
         M = IPVT(K)
         IF (M .EQ. K) GO TO 60
         T = WORK(M)
         WORK(M) = WORK(K)
         WORK(K) = T
   60 CONTINUE
C
      YNORM = 0.0D0
      DO 65 I = 1, N
         YNORM = YNORM + DABS(WORK(I))
   65 CONTINUE
C
C     SOLVE A*Z = Y
C
      CALL SOLVE(NDIM, N, A, WORK, IPVT)
C
      ZNORM = 0.0D0
      DO 70 I = 1, N
         ZNORM = ZNORM + DABS(WORK(I))
   70 CONTINUE
C
C     ESTIMATE CONDITION
C
      COND = ANORM*ZNORM/YNORM
      IF (COND .LT. 1.0D0) COND = 1.0D0
      RETURN
C
C     1-BY-1
C
   80 COND = 1.0D0
      IF (A(1,1) .NE. 0.0D0) RETURN
C
C     EXACT SINGULARITY
C
   90 COND = 1.0D+32
      RETURN
      END
      SUBROUTINE SOLVE(NDIM, N, A, B, IPVT)
C
      INTEGER NDIM, N, IPVT(N)
      DOUBLE PRECISION A(NDIM,N),B(N)
C
C   SOLUTION OF LINEAR SYSTEM, A*X = B .
C   DO NOT USE IF DECOMP HAS DETECTED SINGULARITY.
C
C   INPUT..
C
C     NDIM = DECLARED ROW DIMENSION OF ARRAY CONTAINING A .
C
C     N = ORDER OF MATRIX.
C
C     A = TRIANGULARIZED MATRIX OBTAINED FROM DECOMP .
C
C     B = RIGHT HAND SIDE VECTOR.
C
C     IPVT = PIVOT VECTOR OBTAINED FROM DECOMP .
C
C   OUTPUT..
C
C     B = SOLUTION VECTOR, X .
C
      INTEGER KB, KM1, NM1, KP1, I, K, M
      DOUBLE PRECISION T
C
C     FORWARD ELIMINATION
C
      IF (N .EQ. 1) GO TO 50
      NM1 = N-1
      DO 20 K = 1, NM1
         KP1 = K+1
         M = IPVT(K)
         T = B(M)
         B(M) = B(K)
         B(K) = T
         DO 10 I = KP1, N
             B(I) = B(I) + A(I,K)*T
   10    CONTINUE
   20 CONTINUE
C
C     BACK SUBSTITUTION
C
      DO 40 KB = 1,NM1
         KM1 = N-KB
         K = KM1+1
         B(K) = B(K)/A(K,K)
         T = -B(K)
         DO 30 I = 1, KM1
             B(I) = B(I) + A(I,K)*T
   30    CONTINUE
   40 CONTINUE
   50 B(1) = B(1)/A(1,1)
      RETURN
      END
