C
C
C     LIST OF SUBROUTINES IN THIS FILE 
C
C     CB_SRK : CALCULATE DENSITY SOLUTION OF CUBIC EOS 
C			 (SRK) USING NEWTON-RAPHSON METHOD
C
C*******************************************************************
C
C KDB THEMOPHYSICAL PROPERTIES CALCULATION FORTRAN ROUTINE LIBRARY 
C
C [NAME   ]  CB_SRK 
C
C [TYPE   ]  FORTRAN SUBROUTINE                              
C
C [PURPOSE]  CALCULATE DENSITY SOLUTION OF CUBIC EOS(SRK)
C			USING NEWTON-RAPHSON METHOD
C
C [USAGE  ]  CALL CB_SRK(NC,T,P,RHO,Y,IPHASE,IST) 
C
C [ARGUMENTS]                                                                  
C     NC      : COMPONENT NUMBER (1-50) TO CALCULATE DENSITY(INTEGER, INPUT) 
C     T       : TEMPERATURE IN KELVIN (REAL*8, INPUT) 
C	P		: PRESSURE IN kPa (REAL*8, INPUT)
C	RHO		: DENSITY IN MOL/CM3 (REAL*8, OUTPUT) 
C	Y(50)	: MOLE FRACTION (REAL*8, INPUT)
C	IPHASE	: PHASE IDENTIFIER ( 0-GAS, 1-LIQUID)
C     IST     : STATUS OF CALCULATION (INTEGER, OUTPUT) 
C                = 0	    : NORMAL TERMINATION 
C                = 1001  	: NO CONVERGENCE 
C			   = 1002   : FAILURE OF CALCULATION	
C
C [COMMENTS]
C     
C [REQUIRED COMMON BLOCKS]
C	NONE
C                                                                  
C [REQUIRED SUBROUTINES OR FUNCTIONS] 
C	NONE 
C
C [REFERENCE]
C     RAGUH RAMAN, CHEMICAL PROCESS COMPUTATIONS, 1985,  
C		ELSEVIER APPLEID SCIENCE PUBLISHERS, LONDON AND NEWYORK
C	
C [REVISION INFORMATION]
C     1.REVISED BY Y.S.KIM, KOREA UNIVERSITY, 2002
C*******************************************************************
	
      SUBROUTINE CB_SRK(NC,T,P,RHO,Y,IPHASE,IST)
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	REAL*8 M,Y(50),ALFA(50) 
      COMMON/HC_PROP/WT(50),TB(50),TF(50),TC(50),PC(50),VC(50),ZC(50)
	1	,ACCF(50),WSRK(50),VEST(50),ZRA(50),SOLP(50),VOLP(50)
     2	,QI(50),RI(50),DM(50) 

C
C		INITIALIZATION OF PARAMETERS
C

	IST = 0
	ITMAX = 100
	RF = 0.1 
	AM = 0
	BM = 0
	DO 20 I = 1,NC
	TR = T/TC(I)
	M = 0.48508 + ACCF(I)*(1.55171-0.15613*ACCF(I))
	IF (TR .LE. 1.8) ALFA(I) = (1.+M*(1.-DSQRT(TR)))**2
	IF (TR .GT. 1.8) ALFA(I) = DEXP(2.*M*(1.-TR**0.5))
  20  CONTINUE
      DO 30 I = 1,NC
	BM = BM + Y(I)*TC(I)/PC(I)
	DO 30 J = 1,NC
	AIJ = TC(I)*TC(J)*DSQRT(ALFA(I)*ALFA(J)/(PC(I)*PC(J)))
  30  AM = AM + Y(I)*Y(J)*AIJ
      AM = 0.42748*AM*P/T**2
	BM = 0.08664*BM*P/T

C
C		INITIALIZATION OF DENSITIY
C

	RHOHS = 1.D0/BM
	RHOMC = 1.D0/3.8473221D0/BM
	IF (IPHASE.EQ.0) THEN
		RHO = RHOMC/4.D0
		Z = P/8314.47D0/T/RHO
	ELSE
		IF(T .LT. 200.D0) AAY = 5.D-2
		IF(T .GE. 200.D0) AAY = 3.D-1
		RHO = RHOHS-AAY*(RHOHS-RHOMC)
		Z = P/8314.47D0/T/RHO
	ENDIF

C
C		SOLVE THE SRK EOS USING NEWTON-RAPHSON METHOD
C

	ITER = 0   
  40  F =  -AM*BM - Z*((BM*(BM+1.)-AM)+Z*(1. - Z))
      IF (DABS(F) .LE. 5.D-8) GO TO 50
	FP = -(BM*(BM+1)-AM)-Z*(2.-3.*Z)
	DZ = -F/FP
	IF (DABS(DZ) .GT. Z*RF) DZ = DSIGN(Z*RF,DZ)
	Z = Z + DZ
	ITER = ITER + 1
  
	IF (ITER .LT. ITMAX) THEN
		GO TO 40
	ELSE
		IST = 1001
		RETURN
	ENDIF

  50	CONTINUE 

C		
C		SOLVE QUADRATIC FOR OTHER ROOTS
C		USE MAX. OF ROOTS FOR GAS PHASE AND MIN. FOR LIQUID PHASE
C

      DISCR = (1. +3.*Z)*(1.-Z) + 4.*(BM*(BM+1.)-AM)
	IF ( DISCR.LT.0) THEN
		IST = 1002
		RETURN
	END IF
	Z1 = 0.5*(1.-Z +DSQRT(DISCR))
	Z2 = 1. -Z -Z1
	IF (IPHASE .EQ. 0)    Z = DMAX1(Z, Z1, Z2)
	IF (IPHASE .EQ. 1)    Z = DMIN1(Z, Z1, Z2)
	RHO = P/Z/8314.47D0/T
	RETURN 
      END