      subroutine simul
*      FUNCTION SIMUL( N, A, X, EPS, INDIC, NRC)
C     THIS PROGRAM ORIGINALLY CAME FROM THE BOOK:
C     B. CARNAHAN, H. LUTHER, AND J. WILKES, APPLIED NUMERICAL
C     METHODS,  JOHN WILEY AND SONS, NEW YORK, 1969, P290
C        WHEN INDIC IS NEGATIVE, SIMUL COMPUTES THE INVERSE
C        OF THE N BY N MATRIX A IN PLACE.  WHEN INDIC IS ZERO,
C        SIMUL COMPUTES THE N SOLUTIONS X(1)...X(N) CORRESPONDING
C        TO THE SET OF LINEAR EQUATIONS WITH AUGMENTED MATRIX
C        OF COEFFICIENTS IN THE N BY N+1 ARRAY A AND IN
C        ADDITION COMPUTES TH INVERSE OF THE COEFFICIENT
C        MATRIX IN PLACE AS ABOVE.  IF INDIC IS POSITIVE,
C        THE SET OF LINEAR EQUATIONS IS SOLVED BUT THE
C        INVERSE IS NOT COMPUTED IN PLACE.  THE GAUSS-
C        JORDAN COMPLETE ELIMINATION METHOD IS EMPLOYED
C        WITH THE MAXIMUM PIVOT STRATEGY.  ROW AND
C        COLUMN SUBSCRIPTS OF SUCCESIVE PIVOT ELEMENTS ARE
C        SAVED IN ORDER IN THE IROW AND JCOL ARRAYS
C        RESPECTIVELY.  K IS THE PIVOT COUNTER, PIVOT IS
C        THE ALGEBRAIC VALUE OF THE PIVOT ELEMENT. MAX
C        THE NUMBER OF COLUMNS IN A AND DETER THE DETERMINANT
C        OF THE COEFFICIENT MATRIX.  THE SOLUTIONS ARE
C        DETERMINED IN THE N+1TH COLUMN OF A AND THEN
C        UNSCRAMBLED AND PUT IN ORDER IN X(1)...X(N)
C        USING THE PIVOT SUBSCRIPT INFORMATION AVAILABLE
C        IN THE IROW AND JCOL ARRAYS.  THE SIGN OF THE
C        DETERMINANT IS ADJUSTED, IF NECESSARY, BY DETERMINING
C        IF AN EVEN OR ODD NUMBER OF PAIRWISE INTERCHANGES IS
C        REQUIRED TO PUT THE ELEMENTS OF THE JORD ARRAY IN
C        ASCENDING SEQUENCE WHERE JORD(IROW(I))  =  JCOL(I).
C        IF THE INVERSE IS REQUIRED, IT IS UNSCRAMBLED IN PLACE
C        USING Y(1)...Y(N) AS TEMPORARY STORAGE.  THE VALUE OF
C        THE DETERMINANT IS GIVEN AS THE VALUE OF THE FUNCTION.
C        SHOULD THE POTENTIAL PIVOT OF LARGEST MAGNITUDE BE
C        SMALLER IN MAGNITUDE THAN EPS, THE MATRIX IS
C        CONSIDERED TO BE SINGULAR AND A TRUE ZERO IS RETURNED
C        AS THE VALUE OF THE FUNCTION.
C
      DIMENSION IROW(50),JCOL(50),JORD(50),Y(50)
      COMMON/ALL/N,NUMCYC,ntot
      COMMON/SCRT/E(9),ee(9,9),NH(4),qh(3),f(9),ff(9,10)
      COMMON/IO/NTOUT,NTIN,NEXTN
      common/sim/delta(10),eps,indic,nrc,deter
C
      MAX = ntot
      IO6 = 108
      IF(INDIC.GE.0) MAX = ntot+1
*
*      print*,'ntot is',ntot
*      do l=1,ntot
*        write(ntout,300)(ff(l,m),m=1,max)
*      enddo
*      print*,'eps is',eps
*      print*,'indic is',indic
*      print*,'the max array dimension is',nrc
C
C     .....  IS ntot LARGER THAN 50  .....
      IF(ntot.LE.50) GO TO 5
      WRITE(ntout,200)
      simult = 0.0
      RETURN
C
C     .....  BEGIN ELIMINATION PROCEDURE  .....
5     DETER = 1.
      DO 18 K = 1,ntot
*      do l=1,Ntot
*        write(ntout,300)(ff(l,m),m=1,max)
*      enddo
*      print*,'   '
      KM1 = K - 1
C     .....  SEARCH FOR THE PIVOT ELEMENT  .....
      PIVOT = 0.0
      DO 11 I = 1,ntot
      DO 11 J = 1,ntot
C     .....  SCAN IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBSCRIPTS  ...
      IF(K.EQ.1) GO TO 9
      DO 8 ISCAN = 1,KM1
      DO 8 JSCAN = 1,KM1
      IF(I.EQ.IROW(ISCAN)) GO TO 11
      IF(J.EQ.JCOL(JSCAN)) GO TO 11
8     CONTINUE
9     IF(ABS(ff(i,j)).LE.ABS(PIVOT)) GO TO 11
      PIVOT = ff(I,J)
      IROW(K) = I
      JCOL(K) = J
11    CONTINUE
C
C     .....  INSURE THAT SELECTED PIVOT IS LARGER THAN EPS  .....
      IF(ABS(PIVOT).GT.EPS) GO TO 13
      simult = 0.0
      RETURN
C
C     .....  UPDATE THE DETERMINANT VALUE  .....
13    IROWK = IROW(K)
      JCOLK = JCOL(K)
      DETER = DETER*PIVOT
C     .....  NORMALIZE PIVOT ROW ELEMENTS  .....
      DO 14 J = 1,MAX
      ff(IROWK,J) = ff(IROWK,J)/PIVOT
14    CONTINUE
*      do l=1,Ntot
*        write(ntout,300)(ff(l,m),m=1,max)
*      enddo
*      print*,'   '
C
C     .....  CARRY OUT ELIMINATION AND DEVELOP INVERSE .....
      ff(IROWK,JCOLK) = 1./PIVOT
      DO 18 I = 1,ntot
      AIJCK = ff(I,JCOLK)
      IF(I.EQ.IROWK) GO TO 18
      ff(I,JCOLK) = - AIJCK/PIVOT
      DO 17 J = 1,MAX
      IF(J.NE.JCOLK) ff(I,J) = ff(I,J) - AIJCK*ff(IROWK,J)
17    CONTINUE
18    CONTINUE
*      do l=1,Ntot
*        write(ntout,300)(ff(l,m),m=1,max)
*      enddo
*      print*,'   '
C
C     .....  ORDER SOLUTION VALUES (IF ANY) AND CREATE JORD ARRAY  .....
      DO 20 I = 1,ntot
      IROWI = IROW(I)
      JCOLI = JCOL(I)
      JORD(IROWI) = JCOLI
      IF(INDIC.GE.0) delta(JCOLI) = ff(IROWI,MAX)
20    CONTINUE
*      print*,'the delta solutions determined by simul are: '
*      print*,(delta(i),i=1,ntot)
C
C     .....  ADJUST SIGN OF DETERMINANT  .....
      INTCH = 0
      NM1 = ntot - 1
      DO 22 I = 1,NM1
      IP1 = I + 1
      DO 22 J = IP1,ntot
      IF(JORD(J).GE.JORD(I)) GO TO 22
      JTEMP = JORD(J)
      JORD(J) = JORD(I)
      JORD(I) = JTEMP
      INTCH = INTCH + 1
22    CONTINUE
      IF(INTCH.NE.INTCH/2*2) DETER = - DETER
C
C     .....  IF INDIC IS POSITIVE RETURN WITH RESULTS  .....
24    IF(INDIC.LE.0) GO TO 26
      simult = DETER
      RETURN
C
C     .....  IF INDIC IS NEGATIVE OR ZERO, UNSCRAMBLE THE INVERSE
C            FIRST BY ROWS  .....
26    DO 28 J = 1,ntot
      DO 27 I = 1,ntot
      IROWI = IROW(I)
      JCOLI = JCOL(I)
      Y(JCOLI) = ff(IROWI,J)
27    CONTINUE
      DO 28 I = 1,ntot
      ff(I,J) = Y(I)
28    CONTINUE
*      do l=1,Ntot
*        write(ntout,300)(ff(l,m),m=1,max)
*      enddo
*      print*,'   '
C     .....  THEN BY COLUMNS  .....
      DO 30 I = 1,ntot
      DO 29 J = 1,ntot
      IROWJ = IROW(J)
      JCOLJ = JCOL(J)
      Y(IROWJ) = ff(I,JCOLJ)
29    CONTINUE
      DO 30 J = 1,ntot
      ff(I,J) = Y(J)
30    CONTINUE
*      do l=1,ntot
*        write(ntout,300)(ff(l,m),m=1,max)
*      enddo
*      print*,'   '
*      print*,'the solutions are'
*      write(ntout,300)(delta(i),i=1,ntot)
C
C     .....  RETURN FOR INDIC NEGATIVE OR ZERO  .....
      simult = DETER
      RETURN
C
C     .....  FORMAT FOR OUTPUT STATEMENT  .....
200   FORMAT(10H0N TOO BIG )
 300  format(1h 9f10.6)
      END
