*#######################################################################
*
      PROGRAM RRMAIM
*
*=======================================================================
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*-----------------------------------------------------------------------
      PARAMETER (NCT=50,MERR=5,MDEBUG=100)
      DIMENSION FZ(NCT),EQR(NCT),TC(NCT),PC(NCT),AF(NCT)
      DIMENSION AFZ(NCT),AFY(NCT),AFX(NCT),AEQR(NCT)
      DIMENSION IDEBUG(MDEBUG),SEVMX(MERR),SEV(MERR)
      REAL SECNDS
      LOGICAL TIMING,FOUND
      CHARACTER ANSWER,FSTAT*8
*=======================================================================

C --- DEFINE THE VERSION OF FSTAT MOST APPROPRIATE FOR YOUR SYSTEM

      FSTAT='UNKNOWN'
*     FSTAT='NEW'

C --- PROMPT FOR THE TYPE OF RUN

      WRITE(*,*) 'Is this a timing run (y/N)?'
      READ(*,'(A)') ANSWER
      IF (ANSWER.EQ.'Y' .OR. ANSWER.EQ.'y') THEN
         TIMING = .TRUE.
         WRITE(*,*)
     .      'Estimated MFLOPS on your machine (486DX/33 = 1 MFLOP):'
         READ(*,*) EMFLOP
         IF (EMFLOP.LT.0.0) EMFLOP=1.0
      ELSE
         TIMING = .FALSE.
      ENDIF
      OVRHD=0.0D0

C --- FILE OPERATIONS

      IO1=11
      IO2=12
      IO3=13
      IOUT=0
      IDMP=0
      IDBG=0
      NDEB=0
      DO 10 I=1,MDEBUG
         IDEBUG(I)=0
   10 CONTINUE

      OPEN(IO1,FILE='RRTEST.OUT',STATUS=FSTAT,IOSTAT=ICHK)
      IF(ICHK.NE.0)THEN
         IO1=0
         WRITE(*,*)'ERROR: Could not open RRTEST.OUT'
         STOP
      ENDIF

      IF(TIMING)THEN
         NREPS=MAX(1,NINT(2.0*EMFLOP))
         IPASS=0
         WRITE(*,*)'Now timing overhead...please be patient.'
      ELSE
         NREPS=1
         IPASS=1
         IOUT=IO1
         OPEN(IO2,FILE='RRTEST.ERR',STATUS=FSTAT,IOSTAT=ICHK)
         IF(ICHK.EQ.0)THEN
            IDMP=IO2
         ELSE
            WRITE(*,*)'WARNING: Could not open RRTEST.ERR'
         ENDIF
         INQUIRE(FILE='RRTEST.CHK',EXIST=FOUND)
         IF(FOUND)THEN
            OPEN(IO3,FILE='RRTEST.CHK',STATUS='OLD',IOSTAT=ICHK)
            IF(ICHK.EQ.0)THEN
               READ(IO3,*,IOSTAT=ICHK) (IDEBUG(I),I=1,MDEBUG)
               CLOSE(IO3)
               DO 20 I=1,MDEBUG
                  IF(IDEBUG(I).GT.0) NDEB=I
   20          CONTINUE
               IF(NDEB.GT.0)THEN
                  OPEN(IO3,FILE='RRTEST.DBG',STATUS=FSTAT,IOSTAT=ICHK)
                  IF(ICHK.EQ.0)THEN
                     IDBG=IO3
                     WRITE(IDBG,1300)
                  ELSE
                     WRITE(*,*)'WARNING: Could not open RRTEST.DBG'
                  ENDIF
               ENDIF
            ELSE
               WRITE(*,*)
     .            'WARNING: Could not open RRTEST.CHK (nor RRTEST.DBG)'
            ENDIF
         ENDIF
         WRITE(*,*)'Now executing...please be patient.'
      ENDIF

C --- SET THE PROBLEM CONSTANTS

      EPSF=1.0D-10
      EPSZ=1.0D-20
      EPSK=1.0D-50
      PK=10000.1D0
      TEM=700.0D0

C --- START THE MAIN LOOP

   50 CONTINUE

C --- INITIALIZE THE PROBLEM VARIABLES

      ITEST=0
      KONERR=0
      ITRTOT=0
      ITRMAX=0
      NPROBS=0
      IERTOT=0
      NEGMX=0
      NEGMY=0
      NEGTOT=0
      DO 70 I=1,MERR
         SEVMX(I)=-100.
   70 CONTINUE
      NSEV=0
      SEVTOT=0.D0

      FL=0.5D0
      FV=0.5D0

C --- START THE TIMER AND THEN START GENERATING TEST CASES

      TIME0=SECNDS(0.0)

      DO 500 L=5,NCT,5

         NC=L

         DO 100 I=1,NC
            QMW=14.0D0*FLOAT(I)+2.0D0
            TB=TWUTBP(QMW)
            SG=TWUSGP(TB)
            TC(I)=TWUTCP(TB)
            PC(I)=TWUPCP(TB)
            AF(I)=AFLK(TB,SG,TC(I),PC(I))
  100    CONTINUE

      DO 500 M=-2,12,2

         MM=MAX(0,M)
         MM=MIN(MM,10)
         DO 120 I=1,NC
            IF (I.EQ.1) THEN
               FZ(1)=MAX(EPSZ,ABS(FLOAT(MM)/10.0D0))
            ELSE
               FZ(I)=MAX(EPSZ,ABS((1.0D0-FZ(1))/FLOAT(NC-1)))
            ENDIF
  120    CONTINUE

            SZ=0.0
            DO 130 I=1,NC
               SZ=SZ+FZ(I)
  130       CONTINUE
            DO 140 I=1,NC
               FZ(I)=FZ(I)/SZ
               AFZ(I)=FZ(I)
  140       CONTINUE

         DO 300 K=-2,4
         DO 300 J=1,9
            PRS=FLOAT(J)*10.0D0**K
            IF (PRS.GE.PK) GOTO 300
            A=1.0D0-(PRS/PK)**0.7D0

            EQRMIN=1.0D0
            EQRMAX=1.0D0
            DO 200 I=1,NC
               TR=TEM/TC(I)
               PR=PRS/PC(I)
               EQR(I)=(PC(I)/PK)**(A-1.0D0)*EXP(5.37D0*A*(1.0D0+AF(I))
     .                        *(1.0D0-1.0D0/TR))/PR
               EQRMIN=MIN(EQRMIN,EQR(I))
               EQRMAX=MAX(EQRMAX,EQR(I))
  200       CONTINUE

            IF (EQRMIN.GE.1.0D0-1.0D-15 .OR.
     .          EQRMAX.LE.1.0D0+1.0D-15) THEN
               GOTO 300
            ENDIF

            IF (M.EQ.-2) THEN
               EQR(1)=1.0D0/EPSK
            ELSEIF (M.EQ.12) THEN
               EQR(NC)=EPSK
            ENDIF

            DO 205 I=1,NC
               AEQR(I)=EQR(I)
  205       CONTINUE

C --- START EXTRA FEED CALCULATIONS FOR EQR CLOSE TO 1

            ICHK=0
  210       CONTINUE
            IF (PRS.GE.9000.0) THEN
               IF (M.EQ.-2 .OR. M.EQ.12) GOTO 300
               ICHK=ICHK+1
               IF (ICHK.GT.98) GOTO 300
               IF (ICHK.GT.1) GOTO 235

               EQRMIN=1.0D0
               EQRMAX=1.0D0
               SAFY=0.0D0
               SAFX=0.0D0
               DO 220 I=1,NC
                  AFY(I)=FZ(I)*EQR(I)
                  AFX(I)=FZ(I)/EQR(I)
                  SAFY=SAFY+AFY(I)
                  SAFX=SAFX+AFX(I)
  220          CONTINUE
               DO 230 I=1,NC
                  AFY(I)=AFY(I)/SAFY
                  AFX(I)=AFX(I)/SAFX
                  AEQR(I)=AFY(I)/AFX(I)
                  EQRMIN=MIN(EQRMIN,AEQR(I))
                  EQRMAX=MAX(EQRMAX,AEQR(I))
  230          CONTINUE
               IF (EQRMIN.EQ.1.0D0 .OR. EQRMAX.EQ.1.0D0) GOTO 300
               FVMIN=1.00000001D0/(1.0-EQRMAX)
               FVMAX=0.99999999D0/(1.0-EQRMIN)
               DFV=1.0D0/90.0D0

  235          CONTINUE
               IF (ICHK.EQ.1) THEN
                  AFV=FVMIN
                  AFL=1.0D0-AFV
               ELSEIF (ICHK.EQ.2) THEN
                  AFV=-EPSF
                  AFL=1.0D0-AFV
               ELSEIF (ICHK.EQ.3) THEN
                  AFV=0.0D0
                  AFL=1.0D0-AFV
               ELSEIF (ICHK.EQ.4) THEN
                  AFV=EPSF
                  AFL=1.0D0-AFV
               ELSEIF (ICHK.EQ.5) THEN
                  AFL=EPSF
                  AFV=1.0D0-AFL
               ELSEIF (ICHK.EQ.6) THEN
                  AFL=0.0D0
                  AFV=1.0D0
               ELSEIF (ICHK.EQ.7) THEN
                  AFL=-EPSF
                  AFV=1.0D0-AFL
               ELSEIF (ICHK.EQ.8) THEN
                  AFV=FVMAX
                  AFL=1.0D0-AFV
               ELSE
                  AFV=0.0D0+DFV*FLOAT(ICHK-8)
                  AFL=1.0D0-AFV
               ENDIF

               SFZ=0.0D0
               DO 240 I=1,NC
                  AFZ(I)=AFV*AFY(I) + AFL*AFX(I)
                  IF (AFZ(I).LT.0.0D0) AFZ(I)=EPSZ
                  SFZ=SFZ+AFZ(I)
  240          CONTINUE
               DO 250 I=1,NC
                  AFZ(I)=AFZ(I)/SFZ
  250          CONTINUE

            ENDIF

C --- SEND THIS TEST CASE TO THE SOLVER DRIVER

            ITEST=ITEST+1
            CALL RRTEST(IPASS,ITEST,NREPS,NC,AFZ,AEQR,FL,FV,IOUT,IDMP,
     .                  IDBG,NDEB,IDEBUG,ITER,NEGX,NEGY,SEV,IERR)

C --- TALLY THE RESULTS

            IF (ITER.LT.0) KONERR=KONERR+1
            ITER=ABS(ITER)
            ITRTOT=ITRTOT+ITER
            IF(ITER.GT.ITRMAX) ITRMAX=ITER
            IF (IERR.NE.0) THEN
               NPROBS = NPROBS + 1
               IERTOT = IERTOT + IERR
               NEGMX = MAX(NEGMX,NEGX)
               NEGMY = MAX(NEGMY,NEGY)
               NEGTOT = NEGTOT + NEGX + NEGY
            ENDIF
            DO 280 I = 1, MERR
               IF(SEV(I).GT.SEVMX(I)) SEVMX(I) = SEV(I)
               IF(SEV(I).GT.0.D0)THEN
                  NSEV = NSEV + 1
                  SEVTOT = SEVTOT + SEV(I)
               ENDIF
  280       CONTINUE

C --- GO BACK FOR ANOTHER TEST CASE

            IF (ICHK.NE.0) GOTO 210

  300    CONTINUE

  500 CONTINUE

C --- TESTS COMPLETED; STOP THE TIMER

      TIME = SECNDS(0.0) - TIME0

C --- IF THIS PASS JUST TIMED THE OVERHEAD, THEN REPEAT IT FOR RESULTS

      IF(IPASS.EQ.0)THEN
         OVRHD=TIME
         IPASS=1
         WRITE(*,*)'Now timing solution...please be patient.'
         NSECS=OVRHD*(1.D0+0.520D0*NREPS)/(1.D0+0.024D0*NREPS)
         WRITE(*,510)NSECS,10*NSECS
  510    FORMAT(' Estimated time remaining:',I4,' to',I5,' seconds.')
         GO TO 50
      ELSE
         TIME=TIME-OVRHD
         OVRHD=OVRHD+OVRHD
      ENDIF

C --- OUTPUT THE RESULTS

      SEVMAX = SEVMX(1)
      DO 550 I = 2, MERR
         IF(SEVMX(I).GT.SEVMAX) SEVMAX = SEVMX(I)
  550 CONTINUE
      WRITE(IO1,600)
      WRITE(IO1,610)ITEST,KONERR,NPROBS,IERTOT,NEGTOT,
     .              SEVMAX,SEVTOT/MAX(1,NSEV),
     .              ITRMAX,FLOAT(ITRTOT)/ITEST
      IF(TIMING)WRITE(IO1,620)NREPS,TIME,OVRHD
      WRITE(IO1,630)TIME+OVRHD
      IF (IDMP.GT.0)
     .   WRITE(IDMP,1400) NEGMY,NEGMX,(SEVMX(I),I=1,MERR)
  600 FORMAT(/' THE RACHFORD-RICE CONTEST RESULTS'/1X,52('-'))
  610 FORMAT( ' Total Number of Test Cases               =',I10,
     .       /' Reported Number of Convergence Problems  =',I10,
     .       /' Cases with Material Balance Errors       =',I10,
     .       /' Total Number of Material Balance Errors  =',I10,
     .       /' Number of Negative Mole Fractions        =',I10,
     .       /' Maximum "Severity" of Errors             =',F10.1,
     .       /' Average of All Positive "Severities"     =',F10.1,
     .       /' Maximum Number of Iterations Reported    =',I10,
     .       /' Average Number of Iterations Reported    =',F10.3)
  620 FORMAT( ' Number of Timed Repetitions              =',I10,
     .       /' Solution  Time (seconds)                 =',F10.3,
     .       /' Overhead  Time (seconds)                 =',F10.3)
  630 FORMAT( ' Total Run Time (seconds)                 =',F10.3)
 1300 FORMAT(
     .'* Format:'/
     .'* ITEST  NC  ITER  '/
     .'* FV(initial)  FL(initial)  FV(final)  FL(final) '/
     .'* K(i) '/
     .'* Z(i) '/
     .'* Y(i) '/
     .'* X(i) ')
 1400 FORMAT(/'  Max:',2I6,5F8.1)

      STOP
      END

*#######################################################################
*
      SUBROUTINE RRTEST(IPASS,ITEST,NREPS,NC,FZI,EQRI,FLI,FVI,IOUT,IDMP,
     +                  IDBG,NDEB,IDEBUG,ITER,NEGX,NEGY,SEV,IERR)
*
*=======================================================================
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FZI(*),EQRI(*),IDEBUG(*),SEV(*)
*-----------------------------------------------------------------------
      PARAMETER (NCT=50,MERR=5)
      DIMENSION FZ(NCT),EQR(NCT),FX(NCT),FY(NCT),TOL(MERR),ERR(MERR)
*-----------------------------------------------------------------------
      PARAMETER (TOLR=1.D-15,EPSM=2.22D-16,SMALL=1.D-74,SEV0=-10.0D0)
      PARAMETER (TOLXY=TOLR,TOLK=TOLR,TOLF=TOLR,TOLZ=TOLR)
*=======================================================================
*
*  INITIALIZE FX AND FY, IN CASE IPASS = 0
*
      DO I = 1, NC
         FX(I) = FZI(I)
         FY(I) = FZI(I)
      ENDDO
*
*  REPEAT THE CALCULATION NREPS TIMES, FOR TIMING PURPOSES
*
      DO IREP = 1, NREPS
*
*  COPY FZI AND EQRI INTO SCRATCH ARRAYS, IN CASE RRSOLV CHANGES THEM
*
         DO I = 1, NC
            FZ(I) = FZI(I)
            EQR(I) = EQRI(I)
         ENDDO
*
*  COPY INITIAL GUESSES AND CALL THE SOLVER
*
         FL = FLI
         FV = FVI
         ITER = 0
         IF(IPASS.GT.0) CALL RRSOLV(NC,FZ,EQR,FX,FY,FL,FV,ITER)
      ENDDO
*
*  CALCULATE THE MATERIAL BALANCE ERRORS
*
      NEGX = 0
      NEGY = 0
      DO I = 1, MERR
         TOL(I) = TOLR
         ERR(I) = 0.D0
      ENDDO
      TOL(1) = TOL(1) + NC*EPSM
      TOL(2) = TOL(1)
      DO I = 1, NC
         IF(FX(I).LT.0.D0) NEGX = NEGX + 1
         IF(FY(I).LT.0.D0) NEGY = NEGY + 1
         ERR(1) = ERR(1) + FY(I)
         ERR(2) = ERR(2) + FX(I)
         ERR(3) = MAX(ERR(3),ABS((FL*FX(I)+FV*FY(I))-FZI(I))/
     +               (ABS(FL*FX(I))+ABS(FV*FY(I))+FZI(I)))
         ERR(4) = MAX(ERR(4),ABS(FY(I)-EQRI(I)*FX(I))/
     +               (ABS(FY(I))+ABS(EQRI(I)*FX(I))+SMALL))
      ENDDO
      ERR(5) = ABS(FL + FV - 1)/(ABS(FL)+ABS(FV)+1)
      ERR(1) = ABS(ERR(1)-1.D0)
      ERR(2) = ABS(ERR(2)-1.D0)
*
*  TEST FOR UNMET TOLERANCES
*
      IERR = 0
      IF(NEGX.GT.0) IERR = IERR + 1
      IF(NEGY.GT.0) IERR = IERR + 1
      SEVM = 0.D0
      DO I = 1, MERR
         IF(ERR(I).GT.TOL(I)) IERR = IERR + 1
         SEV(I)=MAX(LOG10(MAX(SMALL,ERR(I))/MAX(SMALL,TOL(I))),SEV0)
         IF(SEV(I).GT.SEVM) SEVM = SEV(I)
      ENDDO
*
*  SUMMARY OUTPUT
*
      IF (IOUT.GT.0) THEN
         IF (ITEST.EQ.1) WRITE(IOUT,1000)
         IF (IERR.EQ.0) THEN
            WRITE(IOUT,1010) ITEST,ITER,FV,FL,IERR
         ELSE
            WRITE(IOUT,1010) ITEST,ITER,FV,FL,IERR,SEVM
         ENDIF
      ENDIF
 1000 FORMAT(/
     .' Summary Report:'/
     .' -------------------------------------------------------------'/
     .'                                                     Maximum  '/
     .'  Test      #                                 #      Positive '/
     .'    #     Iters.       V            L       Errors   Severity '/
     .' -------  ------  -----------  -----------  ------  ----------')
 1010 FORMAT(I7,I8,1X,1P,2E13.4,I6,0P,F11.1)
*
*  DETAILED ERROR OUTPUT
*
      IF (IDMP.GT.0) THEN
         IF (ITEST.EQ.1) WRITE(IDMP,1100)
         WRITE(IDMP,1110) ITEST,NEGY,NEGX,(SEV(I),I=1,MERR)
      ENDIF
 1100 FORMAT(/
     .' Error Report:'/
     .' ----------------------------------------------------------'/
     .'                      Severity = log10(residual/tolerance) '/
     .'        Occurrences        (should never be positive)      '/
     .'        =========== ======================================='/
     .' Test #  y<0   x<0     y       x       z       K      V+L  '/
     .' ------ ----- ----- ------- ------- ------- ------- -------')
 1110 FORMAT(3I6,5F8.1)
*
*  DEBUGGING OUTPUT
*
      IF (IDBG.GT.0) THEN
         DO 1200 I=1,NDEB
            IF (ITEST.EQ.IDEBUG(I)) THEN
               WRITE(IDBG,1400) ITEST,NC,ITER
               WRITE(IDBG,1410) FVI,FLI,FV,FL
               WRITE(IDBG,1410) (EQR(J),J=1,NC)
               WRITE(IDBG,1410) (FZ(J),J=1,NC)
               WRITE(IDBG,1410) (FY(J),J=1,NC)
               WRITE(IDBG,1410) (FX(J),J=1,NC)
               GOTO 1300
            ENDIF
 1200    CONTINUE
 1300    CONTINUE
 1400    FORMAT(3I8)
 1410    FORMAT(1P,4E24.15)
      ENDIF
*
*  MAKE SURE THE RETURNED PHASE FRACTIONS SUM TO 1
*
      IF(FL.LT.FV)THEN
         FLI = FL
         FVI = 1 - FL
      ELSE
         FVI = FV
         FLI = 1 - FV
      ENDIF
      RETURN
      END

C
C---- Twu correlation for molecular weight of normal paraffins based
C---- on normal boiling point.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.
C
      FUNCTION TWUMWP(TB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      T=LOG(TB/(10.44D0-0.0052D0*TB))
      ICTR=0
  100 CONTINUE
      ICTR=ICTR+1
      IF (ICTR.GT.50) RETURN
      TBC=EXP(5.71419D0+2.71579D0*T-0.286590D0*T*T-
     .    39.8544D0/T-0.122488D0/T/T)-24.7522D0*T+35.3115D0*T*T
      F=TBC-TB
      DFDT=(2.71579D0-0.573180D0*T+39.8544D0/T/T+0.244976D0/T**3)*
     .      EXP(5.71419D0+2.71579D0*T-0.286590D0*T*T-
     .      39.8544D0/T-0.122488D0/T/T)-24.7522D0+70.623D0*T
      IF (ABS(F).LT.0.000001D0) GOTO 200
      T=T-F/DFDT
      GOTO 100
  200 CONTINUE
      TWUMWP=EXP(T)
      RETURN
      END
C
C---- Twu correlation for boiling point of normal paraffins based
C---- on molecular weight.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.D0
C
      FUNCTION TWUTBP(QMW)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      T=LOG(QMW)
      TWUTBP=EXP(5.71419D0+2.71579D0*T-0.286590D0*T*T-
     .       39.8544D0/T-0.122488D0/T/T)-24.7522D0*T+35.3115D0*T*T
      RETURN
      END
C
C---- Twu correlation for critical temperature of normal paraffins based
C---- on normal boiling point.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.D0
C
      FUNCTION TWUTCP(TB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      TBS=TB/100.0D0
      TBTC=0.533272D0+0.191017D-3*TB+0.779681D-7*TB*TB-
     .     0.284376D-10*TB**3+0.959468D2/TBS**13
      TWUTCP=TB/TBTC
      RETURN
      END
C
C---- Twu correlation for critical pressure of normal paraffins based
C---- on normal boiling point.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.D0
C
      FUNCTION TWUPCP(TB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      TCP=TWUTCP(TB)
      A=1.0D0-TB/TCP
      TWUPCP=(3.83354D0+1.19629D0*SQRT(A)+34.8888D0*A+
     .       36.1952D0*A*A+104.193D0*A**4)**2
      RETURN
      END
C
C---- Twu correlation for specific gravity of normal paraffins based
C---- on normal boiling point.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.D0
C
      FUNCTION TWUSGP(TB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      TCP=TWUTCP(TB)
      A=1.0D0-TB/TCP
      TWUSGP=0.843593D0-0.128624D0*A-3.36159D0*A**3-13749.5D0*A**12
      RETURN
      END
C
C---- Twu correlation for critical volume of normal paraffins based
C---- on normal boiling point.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.D0
C
      FUNCTION TWUVCP(TB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      TCP=TWUTCP(TB)
      A=1.0D0-TB/TCP
      TWUVCP=1.0D0/(1.0D0-(0.419869D0-0.505839D0*A-1.56436D0*A**3-
     .       9481.70D0*A**14))**8
      RETURN
      END
C
C---- Twu correlation for critical temperature of petroleum fractions
C---- based on normal boiling point and specific gravity.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.D0
C
      FUNCTION TWUTCF(TB,SG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      TCP=TWUTCP(TB)
      SGP=TWUSGP(TB)
      DSG=EXP(5.0D0*(SGP-SG))-1.0D0
      SQTB=SQRT(TB)
      FT=DSG*(-0.362456D0/SQTB+(0.0398285D0-0.948125D0/SQTB)*DSG)
      TWUTCF=TCP*((1.0D0+2.0D0*FT)/(1.0D0-2.0D0*FT))**2
      RETURN
      END
C
C---- Twu correlation for critical volume of petroleum fractions
C---- based on normal boiling point and specific gravity.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.D0
C
      FUNCTION TWUVCF(TB,SG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      VCP=TWUVCP(TB)
      SGP=TWUSGP(TB)
      DSG=EXP(4.0D0*(SGP*SGP-SG*SG))-1.0D0
      SQTB=SQRT(TB)
      FV=DSG*(0.466590D0/SQTB+(-0.182421D0+3.01721D0/SQTB)*DSG)
      TWUVCF=VCP*((1.0D0+2.0D0*FV)/(1.0D0-2.0D0*FV))**2
      RETURN
      END
C
C---- Twu correlation for critical pressure of petroleum fractions
C---- based on normal boiling point and specific gravity.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.D0
C
      FUNCTION TWUPCF(TB,SG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      TCP=TWUTCP(TB)
      PCP=TWUPCP(TB)
      SGP=TWUSGP(TB)
      VCP=TWUVCP(TB)
      VCF=TWUVCF(TB,SG)
      TCF=TWUTCF(TB,SG)
      DSG=EXP(0.5D0*(SGP-SG))-1.0D0
      SQTB=SQRT(TB)
      FP=DSG*(2.53262D0-46.1955D0/SQTB-0.00127885D0*TB+
     .   (-11.4277D0+252.140D0/SQTB+0.00230535D0*TB)*DSG)
      TWUPCF=PCP*(TCF/TCP)*(VCP/VCF)*
     .       ((1.0D0+2.0D0*FP)/(1.0D0-2.0D0*FP))**2
      RETURN
      END
C
C---- Twu correlation for molecular weight of petroleum fractions
C---- based on normal boiling point and specific gravity.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.D0
C
      FUNCTION TWUMWF(TB,SG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      QMWP=TWUMWP(TB)
      SGP=TWUSGP(TB)
      DSG=EXP(5.0D0*(SGP-SG))-1.0D0
      SQTB=SQRT(TB)
      X=ABS(0.0123420D0-0.328086D0/SQTB)
      FM=DSG*(X+(-0.0175691D0+0.193168D0/SQTB)*DSG)
      DMW=LOG(QMWP)*((1.0D0+2.0D0*FM)/(1.0D0-2.0D0*FM))**2
      TWUMWF=EXP(DMW)
      RETURN
      END
C
C---- Twu correlation for boiling point when molecular weight and
C---- specific gravity are known.
C---- Twu, C.H.: Fluid Phase Equilibria, 16 (1984) 137-150.D0
C
      FUNCTION TWUTBF(QMW,SG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      TBEST=1.01D0*TWUTBP(QMW)
      ICTR=0
  100 CONTINUE
      ICTR=ICTR+1
      IF (ICTR.GT.20) THEN
         TWUTBF=TWUTBP(QMW)
         RETURN
      ENDIF
      QMWF=TWUMWF(TBEST,SG)
      F=(QMWF-QMW)/QMW
      IF (ABS(F).LT.0.0001D0) THEN
         TWUTBF=TBEST
         RETURN
      ENDIF
      TBKEEP=TBEST
      IF (ICTR.EQ.1) TBEST=TBEST*(1.0D0-F)
      IF (ICTR.GT.1) TBEST=TBEST-F*(TBOLD-TBEST)/(FOLD-F)
      TBOLD=TBKEEP
      FOLD=F
      GOTO 100
      END
C
C---- Ingolf Soreide Boiling Point Correlation
C
      FUNCTION TBFIS(QMW,SG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      TBFIS=1928.3D0
     .     -1.695D5*EXP(-4.922D-3*QMW-4.7685D0*SG+3.462D-3*QMW*SG)/
     .      QMW**0.03522D0*SG**3.266D0
      RETURN
      END
C
C---- Lee-Kesler correlation for acentric factor based on normal boiling
C---- point and specific gravity.
C---- Lee, B.I. and Kesler, M.G.: AIChE J. 21 (1975) pg. 510.D0
C
      FUNCTION AFLK(TB,SG,TCF,PCF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      TBR=TB/TCF
      TBR6=TBR**6
      DTBR=LOG(TBR)
      PBR=14.696D0/PCF
      AFLK=(LOG(PBR)-5.92714D0+6.09648D0/TBR+
     .     1.28862D0*DTBR-0.169347D0*TBR6)/
     .     (15.2518D0-15.6875D0/TBR-13.4721D0*DTBR+0.43577D0*TBR6)
      RETURN
      END
*
*     Lahey timer routine
*
      real function secnds (start)
      implicit real (a-h,o-z)
      integer*4 iticks
      call timer(iticks)
      secnds=0.01*float(iticks)-start
      return
      end

*
*     Salford timer routine (accurate to 1 microsecond)
*
*      real function secnds (start)
*      implicit real (a-h,o-z)
*      double precision sec
*      logical align
*      align=.false.
*      sec=high_res_clock@(align)
*      secnds=sec-start
*      return
*      end
