C----------------------------------------------------------------------------
      FUNCTION C(T)
C     TRUNCATE INPUT STRING TO 81 CHARACTERS
C     ENTRY CD RETURNS STRING WITH '$' APPENDED AFTER LAST NON-BLANK
C     CHARACTER TO CONFORM TO DISSPLA PACKAGE REQUIREMENTS
C
      PARAMETER (NCHRS=81)
      CHARACTER C*(NCHRS),CD*(NCHRS),T*(*)
      INTEGER STRIM
      C=T
      RETURN
C
      ENTRY CD(T)
C     BE SURE STRING IS TERMINATED BY A $
      C=T
C     LIMIT THE NUMBER OF CHARACTERS TO PARAMETER 'NCHRS'
      NC=MIN0(NCHRS,STRIM(C))
      NCP=MIN0(NCHRS,NC+1)
      IF(C(NC:NC).NE.'$') C(NCP:NCP)='$'
      END
C----------------------------------------------------------------------------
      SUBROUTINE CLEAR(A,N)
C     SET N ELEMENTS OF ARRAY A TO ZERO
      DOUBLE PRECISION A(N)
      IF(N.LE.0) RETURN
      DO    10 I=1,N
10    A(I)=0.0
      RETURN
      END
      SUBROUTINE COMPRS
C     INITIALIZE THE METACODE FILE(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
C     NDEVTXT is rewound here because this is the first graphics code
C     routine.
C
      REWIND NDEVTXT
      WRITE(NDEVTXT,100) 20
100   FORMAT(I2)
      END
C----------------------------------------------------------------------------
      SUBROUTINE SETTRC
C     INITIALIZE /TRAC/ LABELED COMMON
C     COMMON DECK "TRAC" INSERTED HERE
      LOGICAL GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY
      COMMON /TRAC/ GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY,SMT,OSMT
      COMMON/TRAC/ROLD(20),DROLD(20),TOLD,ZDOT,D2Z,RAD,RAD1
      GROUND=.FALSE.
      SURF=.FALSE.
      END
C----------------------------------------------------------------------------
      SUBROUTINE CURVE(X,Y,N,IC)
C     DRAWS A CURVE THROUGH A SET OF DATA POINTS WITH OR WITHOUT
C     MARKER SYMBOLS(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
      DOUBLE PRECISION X(N),Y(N)

      if (n.eq.0) return

      WRITE(NDEVTXT,100) 29,N,IC
      if (n.ge.2) then
           WRITE(NDEVTXT,200) (X(I),Y(I),X(I+1),Y(I+1),I=1,N-1,2)
C..........Handle the case of an odd number of data points.
           if (MOD(n,2).ne.0) WRITE(NDEVTXT,300) X(N),Y(N)
      else if (n.eq.1) then
           WRITE(NDEVTXT,300) X(1),Y(1)
      endif
100   FORMAT(I2,2I5)
200   FORMAT(4E14.6)
300   FORMAT(2E14.6)
      END
C----------------------------------------------------------------------------
      SUBROUTINE DDEND
C     EMPTIES PLOT BUFFER AND RELEASES PLOTTING COMMAND FILE(DDPLOT)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
      COMMON/KNKN/KNBP,KNVC,KNDT
      COMMON/DDLIM/MXIX,MXIY,MNIX,MNIY
C
      WRITE(NDEVTXT,100) -1
100   FORMAT(I2)
      END
C----------------------------------------------------------------------------
      SUBROUTINE DDTEXT(N,TEXT)
C**********************************************************************
C     WRITES A CHARACTER ARRAY PACKED A10(DDPLOT)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
 
      implicit double precision (a-h,o-z)
 
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
      CHARACTER TEXT*(*)
      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY
      COMMON/KNKN/KNBP,KNVC,KNDT
100   FORMAT(A10,2I5,I2,1X,8A10)
      KNDT=KNDT+1
c%      WRITE(ndevtxt) 3,IX,IY

      IX = 0
      IY = 0
      IOR = 0

      write (ndevtxt,33)3,ix,iy
33    format (i2,2i4)
      LN=MIN0(N*10,LEN(TEXT))
      IF(LN.EQ.N*10) THEN
c%         WRITE(NDEVGRP) IOR,N,N*10,(TEXT(I:I),I=1,LN)
      write (ndevtxt,34)ior,n*10,(text(i:i),i=1,ln)
34    format (2i3/80a1)
      ELSE
C%         WRITE(NDEVGRP)IOR,N,N*10,(TEXT(I:I),I=1,LN)
c%     1   ,(' ',I=LN+1,N*10)
      write (ndevtxt,35)ior,n*10,(text(i:i),i=1,ln)
     1     ,(' ',i=ln+1,n*10)
35    format (2i3/80a1/80a1)
      ENDIF
      END

C----------------------------------------------------------------------------
      SUBROUTINE DFCNST(DRMACH,MAX)
      DOUBLE PRECISION DRMACH(MAX)
C     RETURN MAX MACHINE CONSTANTS(OR THE NUMBER AVAILABLE)
C     IN ARRAY DRMACH.

C     THIS MODULE PROVIDES CONSTANTS WHICH ARE EXPLICITLY DEPENDENT
C     ON MACHINE ARCHITECTURE AND SHOULD BE EXAMINED WHEN MAKING ANY
C     MIGRATION OF THE RAY TRACING ROUTINES.
C
C
      DOUBLE PRECISION RMACH(5)
      common /crmach/rmach
C
C     SINGLE PRECISION MACHINE CONSTANTS
C      (SEE DIGITAL SIGNAL PROCESSING, IEEE PRESS 1979 P S-7)
C
C     THE CONSTANTS ARE DERIVED FROM THE FOLLOWING
C     FORM FOR FLOATING POINT NUMBERS: 
C
C        SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C     FOR MOST MACHINES THE BASE B=2.
C     FOR THE CYBER 60-BIT MACHINES
C        T=48, EMIN=-974 AND EMAX=1070.
C     FOR THE PDP-11 AND VAX-11 MACHINES
C        T=24, EMIN=-127 AND EMAX=127
C
C     RMACH(1)  SMALLEST POSITIVE MAGNITUDE = B**(EMIN-1)
C     RMACH(2)  LARGEST MAGNITUDE           = B**EMAX*(1-B**(-T))
C     RMACH(3)  SMALLEST RELATIVE SPACING   = B**(-T)
C     RMACH(4)  LARGEST RELATIVE SPACING    = B**(1-T)
C     RMACH(5)                              = LOG10(B=2)
C
C     CONSTANTS HERE ARE FOR CDC 6000/7000/CYBER SERIES COMPUTERS
C     PLEASE SUBSTITUTE FOR OTHER MACHINES
C
C
C      INTEGER SMALL(2)
C      INTEGER LARGE(2)
C      INTEGER RIGHT(2)
C      INTEGER DIVER(2)
C      INTEGER LOG10(2)
CC
C      EQUIVALENCE (RMACH(1),SMALL(1))
C      EQUIVALENCE (RMACH(2),LARGE(1))
C      EQUIVALENCE (RMACH(3),RIGHT(1))
C      EQUIVALENCE (RMACH(4),DIVER(1))
C      EQUIVALENCE (RMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR CYBER 170/180 COMPUTERS
C
C*****      DATA RMACH(1) / O"00014000000000000000" /
C*****      DATA RMACH(2) / O"37767777777777777777" /
C*****      DATA RMACH(3) / O"16404000000000000000" /
C*****      DATA RMACH(4) / O"16414000000000000000" /
C*****      DATA RMACH(5) / O"17164642023241175720" /
C
C***REFERENCES  FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR
C                 A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-
C                 MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,
C                 PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  R1MACH
C     MACHINE CONSTANTS FOR VAX 11/780
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1), SMALL(2) /        128,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /
C     DATA DIVER(1), DIVER(2) /       9472,           0 /
C     DATA LOG10(1), LOG10(2) /  546979738,  -805796613 /
C
C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB /
C
C   MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING)
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1), SMALL(2) /         16,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /      15552,           0 /
C     DATA DIVER(1), DIVER(2) /      15568,           0 /
C     DATA LOG10(1), LOG10(2) /  1142112243, 2046775455 /
C
C     DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F /
C
C     MACHINE CONSTANTS FOR THE CRAY 1
C
C     DATA SMALL(1) / 201354000000000000000B /
C     DATA SMALL(2) / 000000000000000000000B /
C
C     DATA LARGE(1) / 577767777777777777777B /
C     DATA LARGE(2) / 000007777777777777774B /
C
C     DATA RIGHT(1) / 376434000000000000000B /
C     DATA RIGHT(2) / 000000000000000000000B /
C
C     DATA DIVER(1) / 376444000000000000000B /
C     DATA DIVER(2) / 000000000000000000000B /
C
C     DATA LOG10(1) / 377774642023241175717B /
C     DATA LOG10(2) / 000007571421742254654B /
C
C     The only time this is used is by GET function which is not even
C     called by the PROFIL program.  This RMACH value is included
C     just for "completeness".

c.....smallest relative spacing
      rmach(3) = 2.d0**(-52)

C     The following is equivalenced to BIGVAL in READW function, although
C     it is also probably never going to be used by the PROFIL program.

C.....largest 4-byte integer
      rmach(5) = 2147483647.d0

      CALL RMOVE(DRMACH,RMACH,MIN0(MAX,5))
C
      END
C-----------------------------------------------------------------------
      SUBROUTINE DFSYS
C*****CHANGED THE CHARACTER DECLARATION SO IT WOULD COMPILE. JAH.
C*****
C     This entire routine has been commented out except for the SYSINI
C     entry point (which actually does nothing but is there in case
C     any special system calls are required.	Most of these
C     functions have no importance for PROFIL calculations or plotting.
C
c*****	CHARACTER*(*) DAT, TIM*(*)
c*****	CHARACTER DATX*10
c*****	CHARACTER TIMX*10
C
C     ALL ROUTINES WHICH RELY ON SYSTEM BASED FUNCTIONS, SUCH AS TIME
C     OF DAY OR MAY BENEFIT FROM SYSTEM AVAILABLE ROUTINES FOR BETTER
C     PERFORMANCE, SUCH AS 'RMOVE' HAVE BEEN COLLECTED IN THIS ROUTINE.
C
C*****	    ENTRY RMOVE(X,Y,N)
C*****	RMOVE was changed to a separate subroutine.
C
C     MOVE 'N' COMPONENTS OF ARRAY 'Y' TO ARRAY 'X'.
C     THIS LOGIC WILL HANDLE OVERLAP CASES ONLY FOR X(M)=X(M+N1)
C     WITH N1>0, M=1,N.
C*****      REAL X(N),Y(N)
C
C      DO 10 I=1,N
C10    X(I)=Y(I)
C     RETURN
C
C     ENTRY IMOVE(IX,IY,N)
C*****	    ENTRY IMOVE(X,Y,N)
C***** IMOVE was changed to a separate subroutine.
C
C     SPECIAL ENTRY FOR INTEGER ARRAYS FOR MACHINES HAVING DIFFERENT
C     WORD SIZES FOR INTEGERS THAN REALS.
C     INTEGER IX(N),IY(N)
C
C      IF(N.LE.0) RETURN
C     USE CYBER BUILT-IN ROUTINE
C*****      CALL MOVLEV(Y,X,N)
C
C     THE EQUIVALENT FORTRAN CODING FOR THIS ROUTINE ARE THE FOLLOWING
C     TWO LINES.
C      DO 20 I=1,N
C20    IX(I)=IY(I)
C
C*****	RETURN
C*****	ENTRY SYSSEC(SECS)
C
C     OBTAIN CURENT TIME OF DAY IN SECONDS FROM SYSTEM ROUTINE
C      SECS=0
C     SYSTEM CALL NEEDED FOR THE CYBER NOS OP SYS
C*****      CALL SECOND(SECS)
C     SYSTEM CALL NEEDED FOR CRAY COMPUTER UNDER CTSS
C      CALL TIMEUSED(ICPU,IO,ISYS,MEM)
C      SECS=ICPU/1.E6
C*****	RETURN
C
C*****	ENTRY SYSTIM(TIM)
C
C     OBTAIN CURENT TIME OF DAY AS A CHARACTER STRING FROM SYSTEM ROUTI
C      TIM='TIME'
C     SYSTEM CALL NEEDED FOR THE CYBER NOS OP SYS
C*****      CALL TIME(TIM)
C     SYSTEM CALL NEEDED FOR THE CRAY CTSS OP SYS
C     CALL TIMEDATE(TIM,DATX,MACH)
C*****	RETURN
C
C*****	ENTRY SYSDAT(DAT)
C
C     OBTAIN CURENT DATE AS A CHARACTER STRING FROM SYSTEM ROUTINE
C      DAT='DATE'
C     SYSTEM CALL NEEDED FOR THE CYBER NOS OP SYS
C*****      CALL DATE(DAT)
C     SYSTEM CALL NEEDED FOR THE CRAY CTSS OP SYS
C     CALL TIMEDATE(TIMX,DAT,MACH)
C*****	RETURN
C     STANDARD SYSTEM INITIALIZATION TO BE CALLED FROM MAIN PROGRAM
      ENTRY SYSINI
C     SYSTEM INITIALIZATION CALL FOR CTSS ON CRAY COMPUTER
C     CALL LINK("READ5,UNIT5=TTY,PRINT6,UNIT6=TTY//")
      RETURN
C
C*****	ENTRY MORTEM
C*****	STOP 'POST MORTEM'
      END
C----------------------------------------------------------------------------
      SUBROUTINE DONEPL
C     TERMINATES PLOTTING SESSION(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
C
      WRITE(NDEVTXT,100) 31
100   FORMAT(I2)
      END
C----------------------------------------------------------------------------
      SUBROUTINE ENDPL(N)
C     TERMINATES PLOTTING ON CURRENT PAGE(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
C
      WRITE(NDEVTXT,100) 30,N
100   FORMAT(I2,I3)
      END
C-----------------------------------------------------------------------------
      SUBROUTINE MX1ALF(T1,T2)
C     SPECIFY USE OF ALTERNATE CHARACTER SET NUMBER 1(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
      CHARACTER*10 T1,T2

      WRITE(NDEVTXT,100) 11, T1, T2
100   FORMAT(I2,TR1,A,1X,A1)
      END
C----------------------------------------------------------------------------
      SUBROUTINE MX2ALF(T1,T2)
C     SPECIFY USE OF ALTERNATE CHARACTER SET NUMBER 2(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
      CHARACTER*10 T1,T2

      WRITE(NDEVTXT,200) 12,T1,T2
200   FORMAT(I2,TR1,A,1X,A1)
      END
C----------------------------------------------------------------------------
      FUNCTION NCLOSE(AM,AX,AI,NTICKS,AXIS)
C
C     ADDITIONAL LOGIC TO INSURE THAT DIGITS ARE NOT PACKED
C     TOO CLOSE TOGETHER. RULE OF THUMB: KEEP NUMBER OF DIGITS PER
C     INCH LESS THAN 4.
C
      implicit double precision (a-h,o-z)

      NCLOSE=0
      IF(AI.LE.0.0) RETURN
C
      IF(AI .GT. 0.1*(AX-AM)) THEN
         NCLOSE=1
         RETURN
      ENDIF
C
      NDS=dlog10(dmax1(dabs(AM),dabs(AX)))
C
      NDS=MOD(NDS,5)+3
      IF(AM.LT.0.0) NDS=NDS+1
      NCLOSE=(NDS*(AX-AM)/AI/AXIS-1.0)/4.0+1
      IF(NCLOSE.LE.1) RETURN
      AI=NCLOSE*AI
      NTICKS=4
      RETURN
      END
C-------------------------------------------------------------------------
      FUNCTION ND2B(INDEC)
C
C     CONVERT A NUMERIC DECIMAL DIGIT STRING TO A BIT STRING
C     WITH EACH BIT SET BY A CORRESPONDING DECIMAL DIGIT.
C
      ND2B=0
      IF(INDEC.LE.0) RETURN
      M=INDEC
      MB=1
C
10    IF(MOD(M,10).NE.0) ND2B=ND2B+MB
      MB=MB*2
      M=M/10
      IF(M.GT.0) GO TO 10
      END
C----------------------------------------------------------------------------
      SUBROUTINE NOBRDR
C     SUPPRESSES DRAWING THE PAGE BORDER(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT

      WRITE(NDEVTXT,100) 36
100   FORMAT(I2)
      END
C------------------------------------------------------------------------------
      FUNCTION NUMSTG(V,N,FRM)
C     CONVERTS A NUMERIC VALUE TO A STRING
C     USING FORMAT PROVIDED IN CALL
      CHARACTER NUMSTG*80,FRM*(*)
      INTEGER V(N)
C
      NUMSTG=' '
      WRITE(NUMSTG,FRM) V
      RETURN
      END
C----------------------------------------------------------------------------
      SUBROUTINE OPNREP(IUN,FNAME)
C     THIS OPEN OPERATION ALLOWS FOR AN EXISTING VERSION OF
C     A FILE. IF THE FILE EXISTS IT IS OVERWRITTEN.
C     THE LOGIC USED IS HOPEFULLY MACHINE INDEPENDENT.
      CHARACTER*(*) FNAME
      CALL ZAPFIL(IUN,FNAME)
      OPEN(IUN,FILE=FNAME,STATUS='NEW')
      REWIND IUN
      RETURN
      ENTRY OPNURP(IUN,FNAME)
      CALL ZAPFIL(IUN,FNAME)
      OPEN(IUN,FILE=FNAME,STATUS='NEW',FORM='UNFORMATTED')
      REWIND IUN
      END
C----------------------------------------------------------------------------
      SUBROUTINE ZAPFIL(IUN,FNAME)
C     TEST FOR EXISTENCE OF FILE AND DELETE
      CHARACTER FNAME*(*)
      LOGICAL EX
      INQUIRE(FILE=FNAME,EXIST=EX)
      IF(EX) THEN
         OPEN(UNIT=IUN,FILE=FNAME,STATUS='OLD')
         CLOSE(UNIT=IUN,STATUS='DELETE')
      ENDIF
      END
C----------------------------------------------------------------------------
      SUBROUTINE OVERRD(VAR,TST,DEFT,NFLG,NFVLEQ,NFVLNE)
C
C     OVERRIDE SUPPORT ROUTINE
C     TEST 'VAR' AGAINST 'TST', IF EQUAL SET TO DEFAULT 'DEFT'
C     AND ALSO SET INTEGER FLAG 'NFLG' TO VALUE 'NFVLEQ' ELSE 'NFVLNE'
C
      implicit double precision (a-h,o-z)

      IF(VAR.EQ.TST) THEN
         VAR=DEFT
         NFLG=NFVLEQ
      ELSE
         NFLG=NFVLNE
      ENDIF
      RETURN
      END
C----------------------------------------------------------------------------
      SUBROUTINE PAGE(X,Y)
C     DEFINE LIMITS OF PHYSICAL PAGE(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
      WRITE(NDEVTXT,100) 23,X,Y
100   FORMAT(I2,2F6.2)
      END
C----------------------------------------------------------------------------
      SUBROUTINE PHYSOR(X,Y)
C     DEFINE LOCATION OF PHYSICAL ORIGIN OF PLOTTING AREA(DISSPLA)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
      WRITE(NDEVTXT,100) 22,X,Y
100   FORMAT(I2,2F5.2)
      END
C****************************************************************************
      SUBROUTINE GAUSEL (C,NRD,NRR,NCC,NSF)
C**************************************************************************** 
C     MATRIX INVERSION BY METHOD OF GAUSSIAN PIVOT
C     SOLVES NCC-NRR SETS OF NRR EQUATIONS
C********* SAME AS SUBROUTINE GAUSSEL WRITTEN BY L. DAVID LEWIS ********
C.....This routine is included in the main body of PROFILE even though its
C     only use is when either the CTABLE sound speed model or the TTABLE
C     temperature model is used.

      implicit double precision (a-h,o-z)

      DIMENSION C(NRD,NCC),L(128,2)
C          BITS = 2.**-18
      DATA BITS/3.8146972656d-6/
      NR=NRR
      NC=NCC
      IF(NC.LT.NR.OR.NR.GT.128.OR.NR.LE.0) STOP 5000
C
C          INITIALIZE.
      NSF=0
      NRM=NR-1
      NRP=NR+1
      D=1.
      LSD=1
      DO 1 KR=1,NR
      L(KR,1)=KR
1     L(KR,2)=0
      IF(NR.EQ.1) GO TO 42
C
C          ELIMINATION PHASE.
      DO 41 KP=1,NRM
      KPP=KP+1
      PM=0.d0

      MPN=0
C
C          SEARCH COLUMN KP FROM DIAGONAL DOWN FOR MAX PIVOT.
      DO 2 KR=KP,NR
      LKR=L(KR,1)
      PT=dABS(C(LKR,KP))
      IF(PT.LE.PM) GO TO 2
      PM=PT
      MPN=KR
      LMP=LKR
2     CONTINUE
C
C          IF MAX PIVOT IS ZERO, MATRIX IS SINGULAR.
      IF(MPN.EQ.0) GO TO 9
      NSF=NSF+1
      IF(MPN.EQ.KP) GO TO 3
C
C          NEW ROW NUMBER KP HAS MAX PIVOT.
      LSD=-LSD
      L(KP,2)=L(KP,1)
      L(MPN,1)=L(KP,1)
      L(KP,1)=LMP
C
C          ROW OPERATIONS TO ZERO COLUMN KP BELOW DIAGONAL.
3     MKP=L(KP,1)
      P=C(MKP,KP)
      D=D*P
      DO 41 KR=KPP,NR
      MKR=L(KR,1)
      Q=C(MKR,KP)/P
      IF(Q.EQ.0.d0) GO TO 41
C
C          SUBTRACT Q * PIVOT ROW FROM ROW KR.
      DO 4 LC=KPP,NC
      R=Q*C(MKP,LC)
      C(MKR,LC)=C(MKR,LC)-R
4     IF(dABS(C(MKR,LC)).LT.dABS(R)*BITS) C(MKR,LC)=0.d0
41    CONTINUE
C
C          LOWER RIGHT HAND CORNER.
42    LNR=L(NR,1)
      P=C(LNR,NR)
      IF(P.EQ.0.d0) GO TO 9
      NSF=NSF+1
      D=D*P*LSD
      IF(NR.EQ.NC) GO TO 8
C
C          BACK SOLUTION PHASE.
      DO 61 MC=NRP,NC
      C(LNR,MC)=C(LNR,MC)/P
      IF(NR.EQ.1) GO TO 61
      DO 6 LL=1,NRM
      KR=NR-LL
      MR=L(KR,1)
      KRP=KR+1
      DO 5 MS=KRP,NR
      LMS=L(MS,1)
      R=C(MR,MS)*C(LMS,MC)
      C(MR,MC)=C(MR,MC)-R
5     IF(dABS(C(MR,MC)).LT.dABS(R)*BITS) C(MR,MC)=0.d0
6     C(MR,MC)=C(MR,MC)/C(MR,KR)
61    CONTINUE
C
C          SHUFFLE SOLUTION ROWS BACK TO NATURAL ORDER.
      DO 71 LL=1,NRM
      KR=NR-LL
      MKR=L(KR,2)
      IF(MKR.EQ.0) GO TO 71
      MKP=L(KR,1)
      DO 7 LC=NRP,NC
      Q=C(MKR,LC)
      C(MKR,LC)=C(MKP,LC)
7     C(MKP,LC)=Q
71    CONTINUE
C
C          NORMAL AND SINGULAR RETURNS.  GOOD SOLUTION COULD HAVE D=0.
8     C(1,1)=D
      GO TO 91
9     C(1,1)=0.d0
   91 RETURN
         END
  
