CDate: Thu Nov 26, 1992 5:05 am EST CFrom: michael C EMS: INTERNET / MCI ID: 376-5414 C MBX: [email protected] C CTO: * Curtis Whitson / MCI ID: 522-6476 C **==FSTEST.FOR C C C Multiphase PT-flash program C C Copyright: Michael L. Michelsen C Department of Chemical Engineering C Technical University of Denmark C C Bygning 229, DTH, DK-2800 Lyngby, DENMARK C C C ****************************************************** C * THIS POGRAM MUST NOT BE DISTRIBUTED WITHOUT * C * THE WRITTEN PERMISSION OF MICHAEL L MICHELSEN * C ****************************************************** C C MAIN DIMENSIONING VARIABLES ARE HANDLED BY PARAMETER SPECS: C C MDIM: MAX. NO. OF COMPONENTS, CURRENT VALUE 20 C NFMAX: MAX. NO. OF PHASES, CURRENT NO: 5 (EXCL. PURE SOLIDS) C PROGRAM FSTEST C C SAMPLE MAIN PROGRAM FOR MULTIPHASE FLASH XFLASH C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,NFMAX=5) C C INFIL IS UNITNO. FOR INPUT FILE C PARAMETER (INFIL=21,IOUT=22) DIMENSION FEED(MDIM) , PHASES(NFMAX) , ZFACT(NFMAX) , & ZFRAC(MDIM,5) , ZFEED(MDIM) DIMENSION IFTYP(MDIM) , SOLIDS(MDIM) DIMENSION WKAREA((MDIM*NFMAX)*(MDIM*NFMAX)) DIMENSION LIST(MDIM) LOGICAL SELECT(MDIM) CHARACTER*64 FILNAM CHARACTER*4 FLTP(2) , FASTYP DATA FLTP/' LIQ' , ' VAP'/ CHW WRITE (*,*) ' ENTER 0 FOR CONSOLE INPUT, 1 FOR FILE INPUT ' CHW READ (*,*) INTYP INTYP=1 CHW IF (INTYP .NE. 0) THEN WRITE(*,*) ' ENTER NAME FOR INPUT FILE ' READ(*,'(A)') FILNAM OPEN (INFIL, FILE = FILNAM) OPEN (22, FILE = 'MLM3P.RES',STATUS='UNKNOWN') ENDIF IF (INTYP .EQ. 0) THEN WRITE(*,9010) READ(*,*)MODEL WRITE(*,*)' ENTER NO. OF COMPONENTS (MAX.15) ' READ(*,*)NCOMP ELSE READ(INFIL,*) MODEL READ(INFIL,*) NCOMP ENDIF IF( MODEL .EQ.1 )THEN C C HYBRID MODEL C CALL HYBINP(NCOMP) DO 50 I = 1 , NCOMP SELECT(I) = .TRUE. 50 CONTINUE ELSE C C EQUATION OF STATE C IF (INTYP .EQ. 0) THEN WRITE(*,*)' ENTER MODEL: 0=SRK, 1=PR' READ(*,*)IEQ WRITE(*,*)' ENTER LIST OF COMPONENTS ' READ(*,*)(LIST(I),I=1,NCOMP) WRITE(*,*)' DEBUG UNIT NO? (NONPOSITIVE VALUE: NO DEBUG) ' READ(*,*)ND ELSE READ(INFIL,*) IEQ READ(INFIL,*) ND C C DEBUG PARAMETER VALUE: C C <= 0: NO DETAILED PRINTOUT C 6: CONSOLE PRINTOUT C OTHER: PRINTOUT ON CORRESPONDING UNIT NO., FILENAME DEBUG C ENDIF CALL EOSINP(NCOMP,SELECT,IEQ,LIST,INTYP,INFIL) ENDIF IF( ND.GT.0 .AND. ND.NE.6 )OPEN(ND,FILE='DEBUG') IF (INTYP .EQ. 0) THEN C C OVERALL FEED COMPOSITION C WRITE(*,*)' ENTER FEED ' READ(*,*)(ZFEED(I),I=1,NCOMP) ELSE READ(INFIL,*) (ZFEED(I),I=1,NCOMP) ENDIF IF (INTYP .EQ. 0) THEN WRITE(*,*) ' ITYP=0 USES DEFAULT INITIAL ESTIMATES ' WRITE(*,*) ' ITYP=1 USES USER PROVIDED ' WRITE(*,*) ' ENTER ITYP, T AND P ' READ(*,*) ITYP,T,P C C ITYP = 0: NO INITIAL ESTIMATES AVAILABLE C ITYP > 0: INITIAL ESTIMATES USED C IF (ITYP .GT. 0) THEN WRITE(*,*) ' ENTER NO. OF PHASES ' READ(*,*) NOFAS WRITE(*,*) ' ENTER LIST OF PHASE TYPES (1 = LIQ, -1=VAP )' READ(*,*) (IFTYP(K), K=1,NOFAS) WRITE(*,*) ' ENTER PHASE AMOUNTS ' READ(*,*) (PHASES(K),K=1,NOFAS) WRITE(*,*) ' ENTER SOLID AMT. AND MOLE FRACTIONS : ' DO 305 I=1,NCOMP WRITE(*,306) I 306 FORMAT(' COMP. NO .',I3) READ(*,*) SOLIDS(I),(ZFRAC(I,K),K=1,NOFAS) 305 CONTINUE ENDIF ELSE READ(INFIL,*) ITYP,T,P IF (ITYP .GT. 0) THEN READ(INFIL,*) NOFAS READ(INFIL,*) (IFTYP(K), K=1,NOFAS) READ(INFIL,*) (PHASES(K),K=1,NOFAS) DO 315 I=1,NCOMP READ(INFIL,*) SOLIDS(I),(ZFRAC(I,K),K=1,NOFAS) 315 CONTINUE ENDIF ENDIF 350 CONTINUE MAXLOC = 6 CALL XFLASH(MODEL,ITYP,NCOMP,SELECT,T,P,ZFEED,NOFAS,PHASES,JVAP, & MAXLOC,BETVAL,ZFRAC,ZFACT,IFTYP,SOLIDS,ND,IERCOD,WKAREA) WRITE(22,*)' ERROR CODE = ' , IERCOD WRITE(22,9030)T , P , JVAP , BETVAL DO 400 I = 1 , NCOMP WRITE(22,9020)I , SOLIDS(I) , (ZFRAC(I,J),J=1,NOFAS) 400 CONTINUE WRITE(22,9040) DO 500 J = 1 , NOFAS IF( IFTYP(J).EQ.1 )THEN FASTYP = FLTP(1) ELSE FASTYP = FLTP(2) ENDIF WRITE(22,9050)FASTYP , PHASES(J) , ZFACT(J) 500 CONTINUE C C PROVIDE POSSIBILITY FOR RERUN WITH OLD RESULTS AS INITIAL ESTIMATES C IF (INTYP .EQ. 0) THEN WRITE(*,*) ' ENTER ITYP,T,P (T=0 TO STOP ' READ(*,*) ITYP,T,P ELSE READ(INFIL,*) ITYP,T,P ENDIF IF (T .GT. 0.D0) GOTO 350 C--- TERMINATE RUN IF( ND.GT.0 .AND. ND.NE.6 ) CLOSE(ND) IF (INTYP .NE. 0) CLOSE(INFIL) 9010 FORMAT(/,' ENTER MODEL TYPE ',/,' MODEL =0: SRK/PR ',/, & ' MODEL =1: SPLIT MODEL ',/) 9020 FORMAT(2X,I3,2X,6F12.7) 9030 FORMAT(/,' TEMP = ',F9.3,' PRES = ',F9.3,/,' VAPF = ',I3, & ' FRACTION = ',F9.3,/, & ' CMP SOLID FLUID MOLE FRACTIONS ',/) 9040 FORMAT(/,' TYPE PHASE AMOUNT COMP.FACTOR Z ',/) 9050 FORMAT(1X,A4,2X,F15.8,8X,F10.5) END **==XFLASH.FOR C C C START FLASH C C C SUBROUTINE XFLASH(MODEL,ITYP,NCOMP,SELECT,TEMP,PRES,FEED,NOFAS, & PHASES,JVAP,MAXLOC,BETVAL,ZFRAC,CMPFAC,IFTYP,SOLIDS, & NDEBUG,IRETC,XMAT) C C MULTIPHASE P-T FLASH ROUTINE C C INPUT PARAMETERS: C C MODEL: 0 FOR EOS, 1 FOR SPLIT MODEL C ITYP: TYPE OF CALCULATION; C ITYP=0: NO INITIAL ESTIMATES C ITYP=1: INITIAL ESTIMATES PROVIDED C ITYP=2: INITIAL ESTIMATES PROVIDED, NO STABILITY C ANALYSIS PERFORMED C NCOMP: NO. OF COMPONENTS IN MIXTURE C SELECT: LOGICAL ARRAY THAT DETERMINES WHETHER COMPONENTS ARE USED C AS BASE COMPOSITIONS FOR STABILITY SEARCH C C TEMP, FLASH TEMPERATURE AND PRESSURE C PRES: C FEED: OVERALL MOLAR COMPOSITION OF MIXTURE C NDEBUG: OUTPUT UNIT FOR DEBUG PRINT;NEGATIVE NDEBUG: NO PRINT C MAXLOC: STOPS SEARCH FOR NEW PHASES WHEN MAXLOC IS REACHED C XMAT: SUBPROGRAM WORK AREA OF SIZE AT LEAST NTT*NTT, WHERE C NTT IS MAXFAS*MAXCOMP C C OUTPUT: NOFAS: NO. OF EQUILIBRIUM PHASES C PHASES: PHASE AMOUNTS (MOLES) C JVAP: INDEX FOR VAPOUR PHASE C BETVAL: VAPOUR FRACTION C ZFRAC: PHASE COMPOSITION ARRAY (MOLEFRACS) C CMPFAC: PHASE COMPRESSIBILITY FACTORS C IFTYP: PHASE TYPE, 1=LIQ, -1=VAP C SOLIDS: SOLID COMPONENT AMOUNTS C IRETC: ERROR CODE; 0 = NO ERROR, 1 = NO CONVERGENCE, C 2 = MAX. PHASES EXCEEDED C C FOR ITYP=1,2 THE INITIAL ESTIMATE OF THE NUMBER OF PHASES, C THEIR AMOUNT AND THEIR COMPOSITION IS TAKEN C FROM NOFAS,PHASES AND ZFRAC, RESPECTIVELY. C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,NFMAX=5,NEQM=(NFMAX-1)*MDIM) PARAMETER(NTTAL=NFMAX*MDIM) PARAMETER(ITMAX=10,FASTOL=1.D-8,STABTL=1.D-12,GLIM=1.D-17) DIMENSION FEED(*) , PHASES(*) , ZFRAC(MDIM,*) , CMPFAC(*) , & XMAT(*) DIMENSION IFTYP(*) , SOLIDS(*) C C ALL CALCULATIONS MUST BE PERFORMED IN DOUBLE PRECISION C SECOND PARAMETER BLOCK IS TOLERANCE SPECS; C C IF FASTOL IS CHANGED BY A FACTOR OF 10, STABTL AND GLIM MUST C BE CHANGED BY A FACTOR OF 100 C LOGICAL STRTUP , SELECT(*) , ACTIVE , OUTCLL , EOSTYP DIMENSION Y(MDIM) , YVAL(NTTAL) , GRAD(NEQM) , DLX(MDIM) DIMENSION WORK(5,NEQM) C C THE COMMON-AREAS IN XFLASH ARE USED AS FOLLOWS: C C COMMON/KEEP/: USED TO TRANSFER ESSENTIAL INFORMATION BETWEEN SUBROUTINES. C C NFAS: CURRENT NO. OF PHASES C T,P : TEMPERATURE AND PRESSURE C GNUL: REDUCED GIBBS ENERGY OF INITIAL ESTIMATE C ZFEED: FEED MOLE FRACTIONS C FR: LOG (FUGACITY/P) FOR CURRENT SOLUTION ESTIMATE C XVL: DISTRIBUTION ARRAY; XVL(I,J) IS FRACTION OF FEED AMOUNT C OF COMPONENT I PRESENT IN PHASE J C SFAS: FRACTION OF TOTAL PRESENT IN PHASE C ZVAL: CURRENT COMPRESSIBILITY OF PHASE C IFLU: STATE OF PHASE: 1=LIQ,0=UNKNOWN,-1=VAP C EOSTYP: MODEL TYPE; .TRUE.=EOS; .FALSE.=HYBRID C ISTAB: STATE FOR TRIAL PHASE FOR STABILITY ANALYSIS (-1,0,1) C C COMMON/FPROP/: USED TO RECEIVE THERMODYNAMIC PROPS, I.E. LOG FU-\ C GACITY COEFFICINETS AND DERIVATIVES C C FUG: LOG COMPONENT FUGACITI COEFFICIENTS C FUGX: SCALED COMPOSITION DERIVATIVES OF FUG C C COMMON/AKTIV/: USED TO TRANSFER INFORMATION ABOUT 'TRACE' COMPONENTS C C ACTIVE: MARKS COMPONENT AS TRACE IN ALL PHASES C IVAR: INDEX FOR CURRENT LIST OF NON-TRACE VARIABLES C NACT: TOTAL NO. OF NON-TRACE VARIABLES C EPSV: TOLERANCE FOR TRACE COMPONENT SELECTION C C COMMON/ERRDEB/: SELECTS UNIT FOR EXTENSIVE (DEBUG) PRINTOUT C C C COMMON/OPTIM/: SAVES CURRENT ESTIMATE WITH MIN. GIBBS ENERGY C C VBEST: CURRENT VALUES OF INDEPENDENT VARIABLES (XVL-ARRAY) C GMIN: SMALLEST VALUE OF G/(RT) FOUND SO FAR C COMMON /KEEP / NFAS , T , P , GNUL , ZFEED(MDIM) , FR(MDIM) , & XVL(MDIM,NFMAX) , SFAS(NFMAX) , ZVAL(NFMAX) , & IFLU(MDIM) , EOSTYP , ISTFAS , ISTAB COMMON /FPROP / FUG(MDIM) , FUGX(MDIM,MDIM) COMMON /AKTIV / ACTIVE(MDIM) , IVAR(NEQM+MDIM) , NACT , EPSV COMMON /OPTIM / VBEST(MDIM*(NFMAX+1)) , GMIN COMMON /ERRDEB/ NDPRT COMMON /SOLCL / FUGAR(MDIM) , ISOLC(MDIM) , ISOLPC COMMON /TYPEOS/ MODTYP EXTERNAL FSTAB , PHASN C-----NORMALIZE COMPOSITION IRETC = 0 MODTYP = MODEL NDPRT = NDEBUG STRTUP = .TRUE. EOSTYP = (MODEL.EQ.0) OUTCLL = .FALSE. T = TEMP P = PRES EPSV = 1.D-4 TOTFED = 0.D0 DO 100 I = 1 , NCOMP TOTFED = TOTFED + FEED(I) 100 CONTINUE DO 200 I = 1 , NCOMP ZFEED(I) = FEED(I)/TOTFED 200 CONTINUE C C GET PROPERTIES OF SOLIDS C CALL SOLFUG(NCOMP,T,P,ISOLC,FUGAR) IF( ITYP.GT.0 )THEN C C USER-PROVIDED INITIAL ESTIMATE IS AVAILABLE: C ISOLPC = 0 TOTEST = 0.D0 DO 250 J = 1 , NOFAS TOTEST = TOTEST + PHASES(J) 250 CONTINUE DO 300 I = 1 , NCOMP TOTEST = TOTEST + SOLIDS(I) IF (SOLIDS(I) .GT. 0) ISOLPC = 1 300 CONTINUE DO 350 J = 1 , NOFAS SFAS(J) = PHASES(J)/TOTEST 350 CONTINUE NFAS = NOFAS C--- DIRECT PASS FOR PERTURBATION CALCULATION IF( ITYP.GT.1 )GOTO 1200 GNUL = 0.D0 CALL REFINE(NCOMP,IFTYP,ZFRAC,GVALUE,SOLIDS) IF( NFAS.GT.1 )THEN IF( NDPRT.GT.0 )THEN ZX = 0. WRITE(NDPRT,9010)T , P , ZX ENDIF GNUL = GVALUE GOTO 1200 ENDIF ENDIF ISOLPC = 0 CALL FUGMOD(NCOMP,3,0,IC,T,P,ZX,ZFEED,FUG,FUGX) IF( EOSTYP )THEN IFLU(1) = 0 ELSEIF( IC.GT.0 )THEN IFLU(1) = 1 ELSE IFLU(1) = -1 ENDIF PHASES(1) = 1.D0 CMPFAC(1) = ZX EPSV = 1.D-4 C----EPSV IS USED TO SELECT 'TRACE' COMPONENTS C IF WE HAVE MARGINAL INSTABILITY,EPSV IS DECREASED C IF( NDPRT.GT.0 )WRITE(NDPRT,9010)T , P , ZX GNUL = 0. C-----GIBBS ENERGY OF SINGLE PHASE SFAS(1) = 1. DO 400 I = 1 , NCOMP C------XVL(I,J) IS FRACTION OF COMPONENT I IN PHASE J SOLIDS(I) = 0.D0 XVL(I,1) = 1. C------ZERO MOLE FRACTION NOT PERMITTED (LOG-TERMS) ZFEED(I) = ZFEED(I) + 1.D-20 ZFRAC(I,1) = ZFEED(I) DLX(I) = DLOG(ZFEED(I)) FR(I) = DLX(I) + FUG(I) IF( FUGAR(I).LT.FR(I) )THEN ISOLPC = 1 SOLIDS(I) = 1.D-16 ENDIF GNUL = GNUL + ZFEED(I)*FR(I) 400 CONTINUE IF( NDPRT.GT.0 )WRITE(NDPRT,9020)(I,ZFEED(I),I=1,NCOMP) NFAS = 1 IF( ISOLPC.GT.0 )GOTO 1200 C C PURE IS PURE COMPONENT BASED STABILITY TEST 500 CALL PURE(STRTUP,NCOMP,SELECT,ZFRAC,Y,YV1) C C NEGATIVE YV1: INSTABILITY VERIFIED C POSITIVE YV1: INSTABILITY POSSIBLE, FURTHER TESTING REQUIRED C IF( YV1.LE.-.01 )THEN IF( NDPRT.GT.0 )WRITE(NDPRT,9060) GOTO 900 ENDIF DO 600 I = 1 , NCOMP C C MODIFIED TO SQUARE ROOT FORM C C Y(I) = 2.D0*DSQRT(Y(I)) 600 CONTINUE XLAM = 0.25 IF( NDPRT.GT.0 )WRITE(NDPRT,9030) DO 700 I = 1 , NCOMP IVAR(I) = I ACTIVE(I) = .TRUE. 700 CONTINUE NACT = NCOMP CALL MURML(NCOMP,1,ITMAX,XLAM,STABTL,FUN,Y,GRAD,XMAT,WORK,FSTAB) DO 800 I = 1 , NCOMP Y(I) = Y(I)**2/4 800 CONTINUE C C SET BACK AGAIN C IF( DABS(FUN).LT..01 )EPSV = 1.D-3*DSQRT(DABS(FUN)) IF( NDPRT.GT.0 )WRITE(NDPRT,9040)NACT IF( FUN.GT.-FASTOL )THEN C C TRY ONE-DIMENSIONAL SEARCH FROM PREVIOUS SOLUTION C ISTAB = 1 CALL ONEDIR(NCOMP,ZFRAC,YV1,Y) IF( YV1.GE.-FASTOL )THEN C----STILL NO INDICATION OF INSTABILITY; STOP HERE IF( NDPRT.GT.0 )WRITE(NDPRT,9050) IF( .NOT.OUTCLL )CALL OUT(DLX,NCOMP,IFTYP,ZFRAC,CMPFAC, & PHASES,SOLIDS,JVAP,BETVAL) GOTO 1400 ENDIF ENDIF IF( NDPRT.GT.0 )WRITE(NDPRT,9060) 900 NFAS = NFAS + 1 IF( NFAS.GT.NFMAX )THEN IRETC = 2 NFAS = NFMAX GOTO 1400 ENDIF DO 1000 J = 2 , NFAS JN = NFAS + 2 - J IFLU(JN) = IFLU(JN-1) SFAS(JN) = SFAS(JN-1) DO 950 I = 1 , NCOMP ZFRAC(I,JN) = ZFRAC(I,JN-1) 950 CONTINUE 1000 CONTINUE DO 1100 I = 1 , NCOMP ZFRAC(I,1) = Y(I) 1100 CONTINUE SFAS(1) = 1.D-16 IFLU(1) = ISTAB C-----LIMIT STEPSIZE 1200 XLAM = .2 IF( NFAS.EQ.2 )XLAM = .5 CALL DSUB(NCOMP,ZFRAC,GNORM,GLIM,NFAS,SOLIDS) IF( GNORM.GE.GLIM )THEN SOLSUM = 0. DO 1250 I = 1 , NCOMP SOLSUM = SOLSUM + SOLIDS(I) 1250 CONTINUE ISL = 0 IF( SOLSUM.GT.0. )THEN ISOLPC = 1 ISL = 1 NFAS = NFAS + 1 SFAS(NFAS) = SOLSUM ELSE ISOLPC = 0 ENDIF M = (NFAS-1)*NCOMP DO 1300 I = 1 , M C YVAL(I)=VBEST(I) IF( I.LE.NCOMP )ACTIVE(I) = .TRUE. IVAR(I) = I 1300 CONTINUE NACT = M C-----CHECK FOR REDUCTION OF NUMBER OF PHASES CALL MURML(M,1,ITMAX,XLAM,GLIM,FUN,YVAL,GRAD,XMAT,WORK,PHASN) IF( XLAM.GT.1.D0 )THEN IRETC = 1 ELSE IRETC = 0 ENDIF IF( ISL.GT.0 )THEN IF( SFAS(NFAS).EQ.0.D0 )THEN NFAS = NFAS - 1 ISOLPC = 0 DO 1310 I = 1 , NCOMP SOLIDS(I) = 0.D0 1310 CONTINUE ELSE DO 1320 I = 1 , NCOMP SOLIDS(I) = VBEST(NCOMP*(NFAS-1)+I)*ZFEED(I) 1320 CONTINUE NFAS = NFAS - 1 ENDIF ENDIF ENDIF EPSV = 1.D-4 IF( NDPRT.GT.0. )WRITE(NDPRT,9040)NACT C-----PRINT RESULTS OUTCLL = .TRUE. CALL OUT(DLX,NCOMP,IFTYP,ZFRAC,CMPFAC,PHASES,SOLIDS,JVAP,BETVAL) C-----AND RETURN FOR ADDITIONAL TESTING IF( ITYP.LE.1 .AND. NFAS .LT. MAXLOC)GOTO 500 1400 DO 1500 J = 1 , NFAS PHASES(J) = PHASES(J)*TOTFED 1500 CONTINUE DO 1600 I = 1 , NCOMP SOLIDS(I) = SOLIDS(I)*TOTFED 1600 CONTINUE NOFAS = NFAS RETURN 9010 FORMAT(/,' PT-FLASH, T = ',F8.2,'KEL, P = ',F8.2,' ATM, Z =', & F10.5) 9020 FORMAT(/,' COMPOSITION ',/,(I3,2X,F12.8)) 9030 FORMAT(' ') 9040 FORMAT(/,' VARIABLES ACTIVE ',I4/) 9050 FORMAT(/,' SOLUTION IS NOW STABLE * * * ') 9060 FORMAT(/,' SYSTEM IS STILL UNSTABLE, PHASE SPLIT PERFORMED '/) END **==OUT.FOR C C INTERMEDIATE OUTPUT ROUTINE C SUBROUTINE OUT(XMOL,NCOMP,IFTYP,ZFRAC,ZCOMP,PHASAM,SOLIDS,JVAP, & BETVAL) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C SUBROUTINE OUT OUTPUTS CURRENT RESULTS AND RESETS THE C MIXTURE FUGACITIES FOR A NEW STABILITY ANALYSIS C C XMOL: WORKSPACE FOR MOLE FRACTIONS C NCOMP: NO. OF COMPONENTS IN MIXTURE (I) C IFTYP: PHASE TYPE (LIQ=1,VAP=-1) (O) C ZFRAC: ARRAY OF MOLE FRACTIONS (O) C ZCOMP: COMPRESSIBILITY FACTORS (O) C PHASAM: PHASE AMOUNT (RELATIVE) (O) C SOLIDS: VECTOR OF SOLIDS (I) C JVAP: INDEX FOR VAPOUR PHASE (O) C BETVAL: VAPOUR PHASE FRACTION (O) C PARAMETER(MDIM=20,NFMAX=5) DIMENSION IFTYP(*) , XMOL(*) , ZFRAC(MDIM,*) , ZCOMP(*) , & PHASAM(*) DIMENSION SOLIDS(*) LOGICAL SOLFAS , EOSTYP COMMON /ERRDEB/ NDPRT COMMON /KEEP / NFAS , T , P , GNUL , ZFEED(MDIM) , FR(MDIM) , & XVL(MDIM,NFMAX) , SFAS(NFMAX) , ZVAL(NFMAX) , & IFLU(MDIM) , EOSTYP , ISTFAS , ISTAB COMMON /FPROP / FUG(MDIM) , FUGX(MDIM,MDIM) SOLFAS = .FALSE. BETVAL = 0.D0 JVAP = 0 DO 200 J = 1 , NFAS SUM = 0. DO 50 I = 1 , NCOMP XMOL(I) = XVL(I,J)*ZFEED(I)/SFAS(J) SUM = SUM + XMOL(I) 50 CONTINUE SUM = DLOG(SUM) CALL FUGMOD(NCOMP,1,IFLU(J),IC,T,P,ZX,XMOL,FUG,FUGX) ZVAL(J) = ZX IF( .NOT.EOSTYP )THEN IFTYP(J) = IFLU(J) ELSEIF( IC.GT.0 )THEN IFTYP(J) = 1 ELSE IFTYP(J) = -1 ENDIF DO 100 I = 1 , NCOMP IF( SOLIDS(I).NE.0.D0 )SOLFAS = .TRUE. IF( XMOL(I).GT.0.D0 )THEN XMOL(I) = DLOG(XMOL(I)) FRI = XMOL(I) - SUM + FUG(I) ELSE FRI = 1000.D0 ENDIF IF( J.EQ.1 )THEN FR(I) = FRI ELSE IF( FRI.LT.FR(I) )FR(I) = FRI ENDIF 100 CONTINUE 200 CONTINUE INCN = 0 IF( NDPRT.GT.0 )WRITE(NDPRT,9010)NFAS IF( NDPRT.GT.0 )WRITE(NDPRT,9020)(J,SFAS(J),J=1,NFAS) IF( NDPRT.GT.0 )WRITE(NDPRT,9030)(J,ZVAL(J),J=1,NFAS) IF( NDPRT.GT.0 )WRITE(NDPRT,9040)(J,J=1,NFAS) C C IDENTIFY VAPOUR PHASE,IF ANY C DO 300 J = 1 , NFAS IF( IFTYP(J).LT.0 )THEN INCN = INCN + 1 JVAP = J BETVAL = SFAS(J) ENDIF 300 CONTINUE IF( INCN.GT.1 )WRITE(NDPRT,*)' MORE THAN ONE VAP ' DO 400 I = 1 , NCOMP DO 350 J = 1 , NFAS ZFRAC(I,J) = XVL(I,J)*ZFEED(I)/SFAS(J) 350 CONTINUE IF( NDPRT.GT.0 )WRITE(NDPRT,9050)I , (ZFRAC(I,J),J=1,NFAS) 400 CONTINUE DO 500 J = 1 , NFAS ZCOMP(J) = ZVAL(J) PHASAM(J) = SFAS(J) 500 CONTINUE IF( SOLFAS )THEN IF( NDPRT.GT.0 )WRITE(NDPRT,9060) DO 550 I = 1 , NCOMP IF( SOLIDS(I).GT.0 )THEN IF( NDPRT.GT.0 )WRITE(NDPRT,9070)I , SOLIDS(I) ENDIF 550 CONTINUE ENDIF IF( NDPRT.GT.0 )WRITE(22,9080)JVAP , BETVAL RETURN 9010 FORMAT(/,' FLASH RESULT ',I3,' PHASES * * * * ') 9020 FORMAT(/,' TOTAL AMOUNT ',4(1X,I2,F12.8)) 9030 FORMAT(' COMPRESS ',4(1X,I2,F12.8)) 9040 FORMAT(/,' COMPOSITION ',4(4X,I2,9X)) 9050 FORMAT(' X(',I2,') ',4(1X,G14.7)) 9060 FORMAT(' ') 9070 FORMAT(' COMP. NO. ',I3,' SOLID AMOUNT ',F10.7) 9080 FORMAT(' VAPOUR PHASE = ',I3,' AMOUNT: ',F10.7,/) END **==REFINE.FOR SUBROUTINE REFINE(NCOMP,IFTYP,ZFRAC,GVALUE,SOLIDS) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C REFINE PROVIDES A RESTART CAPABILITY WHERE THE RESULT OF A C PREVIOUS CALCULATION (OR SOME OTHER EXTERNALLY GENERATED C COMPOSITION ESTIMATE) IS USED TO INITIATE CALCULATIONS WITH C MORE THAN ONE PHASE C C NCOMP: NO. OF COMPONENTS IN MIXTURE C IFTYP: ARRAY OF PHASE TYPE INDICATORS; C USE 0 IF AN EOS-MODEL IS USED, AND C 1 FOR LIQUID, -1 FOR VAPOUR WITH A MIXED MODEL C ZFRAC: ARRAY OF PHASE MOLE FRACTIONS C C OUTPUT: GVALUE IS THE GIBBS ENERGY OF THE REVISED ESTIMATE C SOLIDS: SOLID PHASE AMOUNTS C PARAMETER(MDIM=20,NFMAX=5) DIMENSION IFTYP(*) , ZFRAC(MDIM,*) , SOLIDS(*) DIMENSION CR(MDIM,NFMAX+1) DIMENSION IPOS(MDIM) LOGICAL SOLIN , EOSTYP COMMON /KEEP / NFAS , T , P , GNUL , ZFEED(MDIM) , FR(MDIM) , & XVL(MDIM,NFMAX) , SFAS(NFMAX) , ZVAL(NFMAX) , & IFLU(MDIM) , EOSTYP , ISTFAS , ISTAB COMMON /SOLCL / FUGAR(MDIM) , ISOLC(MDIM) , ISOLPC SOLIN = (ISOLPC.GT.0) C C SET STATES TEMPORARILY C DO 100 J = 1 , NFAS IFLU(J) = IFTYP(J) 100 CONTINUE CALL MAXFRC(NCOMP,NFAS,ZFRAC,IPOS) CALL NTOPI(NCOMP,IPOS,ZFRAC,GTOT,CR,GNORM,SOLIDS) CALL FSPLIT(NCOMP,NFAS,ZFEED,SFAS,SOLIDS,CR,FUGAR,ZFRAC,20,ICODE, & ISOLC,SOLIN) C C C LOOP TO REMOVE PHASES PRESENT IN ZERO AMOUNT C 200 DO 300 K = 1 , NFAS IF( SFAS(K).EQ.0.D0 )GOTO 500 300 CONTINUE C C ALL PHASES PRESENT IN NONZERO AMOUNTS C ISOLPC = 0 DO 400 I = 1 , NCOMP IF( SOLIDS(I).GT.0.D0 )ISOLPC = 1 400 CONTINUE IF( NFAS.GT.1 .OR. ISOLPC.GT.0 )THEN CALL MAXFRC(NCOMP,NFAS,ZFRAC,IPOS) CALL NTOPI(NCOMP,IPOS,ZFRAC,GTOT,CR,GNORM,SOLIDS) GVALUE = GTOT ENDIF C C RESET FOR EOS C IF( EOSTYP )THEN DO 450 J = 1 , NFAS IFLU(J) = 0 450 CONTINUE ENDIF RETURN 500 NFAS = NFAS - 1 IF( K.LE.NFAS )THEN C C NOT LAST PHASE REMOVED; OTHERS HAVE TO MOVE DOWN C DO 550 J = K , NFAS IFLU(J) = IFLU(J+1) SFAS(J) = SFAS(J+1) DO 520 I = 1 , NCOMP ZFRAC(I,J) = ZFRAC(I,J+1) 520 CONTINUE 550 CONTINUE ENDIF GOTO 200 END **==DSUB.FOR SUBROUTINE DSUB(NCOMP,ZFRAC,GNORM,GLIM,NFA,SOLIDS) C THE DSUB-ROUTINE ATTEMPTS TO CONVERGE THE MULTIPHASE FLASH BY C DIRECT SUBSTITUTIN, ACCELRATED BY THE GENERAL DOMINANT EIGEN- C VALUE METHOD. WE USE 3 EXTRAPOLATION STEPS, EACH AFTER 5 STEPS C OF DIRECT SUBSTITUTION C C DSUB ALSO TESTS FOR ELIMINATION OF PHASES DURING THE ITERATIVE C CALCULATION C C INPUT: C NCOMP: NO. OF COMPONENTS IN MIXTURE C ZFRAC: ARRAY OF PHASE MOLE FRACTIONS C GLIM: GRADIENT NORM TOLERANCE C C OUTPUT C GNORM: GRADIENT NORM ON EXIT C NFA: NO. OF PHASES ON EXIT (ROUTINE MAY ELIMINATE PHASES) C SOLIDS: CALCULATED AMOUNTS OF SOLIDS C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,NFMAX=5) DIMENSION ZFRAC(MDIM,*) , SOLIDS(*) DIMENSION CROLD(MDIM,NFMAX+1) DIMENSION CR(MDIM,NFMAX+1) , SSAVE(MDIM) , IPOS(MDIM) DIMENSION WORK(MDIM*(NFMAX+1),3) C C ARRAY DIMENSIONED BELOW NEED NO CHANGE IN SIZE C THEY ARE WORK AREAS FOR THE GDOM-CALCULATION C DIMENSION B(3,3) , BAU(3) , IP(3) LOGICAL SOLIN , EOSTYP COMMON /ERRDEB/ NDPRT COMMON /KEEP / NFAS , T , P , GNUL , ZFEED(MDIM) , FR(MDIM) , & XVL(MDIM,NFMAX) , SFAS(NFMAX) , ZVAL(NFMAX) , & IFLU(MDIM) , EOSTYP , ISTFAS , ISTAB COMMON /OPTIM / VBEST(MDIM*(NFMAX+1)) , GMIN COMMON /SOLCL / FUGAR(MDIM) , ISOLC(MDIM) , ISOLPC C---M IS SPECIFIER FOR GDOM 100 SOLIN = (ISOLPC.GT.0) LMAX = 4 200 NTOT = NCOMP*(NFAS+1) NFA = NFAS C C----GET FUGACITY COEFFICIENTS FROM INITIAL ESTIMATE CALL MAXFRC(NCOMP,NFAS,ZFRAC,IPOS) CALL NTOPI(NCOMP,IPOS,ZFRAC,GTOT,CR,GNORM,SOLIDS) IF( .NOT.SOLIN .AND. ISOLPC.GT.0 )GOTO 100 C---SAVE FUGACITY COEFFICIENTS CALL COPAR(NCOMP,NFAS+1,MDIM,CR,CROLD) C---FUGACITY COEFFICIENTS SAVED IN CROLD C----AND SAVE VECTOR IF( NDPRT.GT.0 )WRITE(NDPRT,9010) DO 600 L = 1 , LMAX ISO = 0 ICOR = 0 250 ISO = ISO + 1 DO 350 KK = 1 , 5 C C SAVE PHASE AMOUNTS AND CALCULATE NEW DISTRIBUTION C C---CR MODIFIED C---CROLD HOLDS OLD FUGACITY COEFFICIENTS C C CONTINUE WITH EVALUATION C CALL FSPLIT(NCOMP,NFAS,ZFEED,SFAS,SOLIDS,CR,CR(1,NFAS+1), & ZFRAC,20,ICODE,ISOLC,SOLIN) C C GET NEW GIBBS ENERGY AND NEW FUGACITY COEFFICIENTS C C-- CR RECALCULATED CALL NTOPI(NCOMP,IPOS,ZFRAC,GTOT,CR,GNORM,SOLIDS) C C CHECK IF EXTRAPOLATE IMPROVES G; OTHERWISE REDUCE STEP C SAVE G=VALUE FIRST TIME C IF( L.EQ.1 .AND. KK.EQ.1 )GMIN = GTOT IF( KK.EQ.1. .AND. GTOT.GT.GMIN+1.D-10 )THEN DO 260 I = 1 , NCOMP DO 255 J = 1 , NFAS + 1 C-------RESTART FROM LAST VALUE, IF EXTRAPOLATE IS WORSE C IN FIRST STEPS, TRY TO REDUCE, BUT IN THE END, JUST EXIT C IF( ISO.GT.2 .OR. L.EQ.LMAX )THEN CR(I,J) = CROLD(I,J) ELSE CR(I,J) = (CR(I,J)+CROLD(I,J))/2.D0 ENDIF 255 CONTINUE 260 CONTINUE ICOR = ICOR + 1 IF( ICOR.LE.4 )GOTO 250 GOTO 600 ELSEIF( GTOT.LT.GMIN+1.D-10 )THEN C C NEW STEP IS 'BETTER' WITHIN TOLERANCE; SAVE IT C GMIN = GTOT DO 270 J = 1 , NFAS DO 265 I = 1 , NCOMP IF( J.EQ.1 )SSAVE(I) = SOLIDS(I) VBEST(I+(J-1)*NCOMP) = XVL(I,J) 265 CONTINUE 270 CONTINUE IF( .NOT.SOLIN .AND. ISOLPC.GT.0 )GOTO 100 C TEST CONVERGENCE C C C TEST ELIMINATION OF PHASE C JOUT = 0 DO 280 J = 1 , NFAS IF( SFAS(J).EQ.0.D0 )JOUT = J 280 CONTINUE IF( JOUT.NE.0 )THEN IF( NDPRT.GT.0 )WRITE(NDPRT,9020) NFAS = NFAS - 1 C----LAST PHASE ELIMINATED ? IF( JOUT.LE.NFAS )THEN C----NO, A MIDDLE PHASE; REARRANGE C----- RESET YIELD ARRAYS TO SMALLER NUMBER OF PHASES SFAS(JOUT) = SFAS(NFAS+1) IFKP = IFLU(JOUT) IFLU(JOUT) = IFLU(NFAS+1) IFLU(NFAS+1) = IFKP DO 282 I = 1 , NCOMP EE = ZFRAC(I,JOUT) ZFRAC(I,JOUT) = ZFRAC(I,NFAS+1) ZFRAC(I,NFAS+1) = EE EE = VBEST(I+(JOUT-1)*NCOMP) VBEST(I+(JOUT-1)*NCOMP) = VBEST(I+NFAS*NCOMP) VBEST(I+NFAS*NCOMP) = EE 282 CONTINUE ENDIF GOTO 200 ENDIF IF( GNORM.LT.GLIM )GOTO 700 IF( L.EQ.LMAX )GOTO 600 ENDIF DO 300 I = 1 , NCOMP DO 290 J = 1 , NFAS + 1 GRADI = CR(I,J) IV = I + (J-1)*NCOMP IF( KK.GT.2 )WORK(IV,KK-2) = GRADI - CROLD(I,J) CROLD(I,J) = GRADI 290 CONTINUE 300 CONTINUE 350 CONTINUE C C 140 IS THE END OF THE DS-STEPS FOR THE EXTRAPOLATION CYCLE C C NEXT COMES THE EXTRAPOLATION CALCULATION C DO 400 I = 1 , 3 DO 380 J = I , 3 BIJ = 0. DO 360 K = 1 , NTOT BIJ = BIJ + WORK(K,I)*WORK(K,J) 360 CONTINUE B(I,J) = BIJ B(J,I) = BIJ 380 CONTINUE 400 CONTINUE C------CONVERGENCE; DANGER OF UNDERFLOW IN EXTRAPOLATION IF( B(3,3).LT.1.D-16 )GOTO 700 DO 450 J = 1 , 2 B(3,J) = 0. BAU(J) = 0. 450 CONTINUE BAU(3) = 1. B(3,3) = 1. C C SOLVE LINEAR EQNS FOR EXTRAPOLATE C CALL LU(3,3,IP,B) CALL BACK(3,3,IP,B,BAU) SNY = 0. DO 500 J = 1 , 3 SNY = SNY + BAU(J) 500 CONTINUE AFAC = -BAU(1) - BAU(2) BFAC = -BAU(1) SLIM = .5 CPR = 0. SMAX = SLIM DO 550 J = 1 , NFAS + 1 DO 520 I = 1 , NCOMP IN = (J-1)*NCOMP + I YNEWI = (AFAC*WORK(IN,3)+BFAC*WORK(IN,2))/SNY CR(I,J) = CR(I,J) + YNEWI DS = YNEWI IF( DABS(DS).GT.SMAX )SMAX = DABS(DS) CPR = CPR + DS*WORK(IN,3) 520 CONTINUE 550 CONTINUE IF( NDPRT.GT.0 )WRITE(NDPRT,9030)GMIN , (BAU(J),J=1,3) IF( CPR.LT.0. .OR. SMAX.GT.SLIM )THEN C C CHECK LENGTH OF STEP AND THAT DIRECTION IS OK C FAC = SLIM/SMAX IF( CPR.LT.0. )FAC = -FAC DO 580 I = 1 , NCOMP DO 560 J = 1 , NFAS + 1 CR(I,J) = CROLD(I,J) + FAC*(CR(I,J)-CROLD(I,J)) 560 CONTINUE 580 CONTINUE ENDIF C C 230 MARKS THE END OF THE OUTER LOOP C 600 CONTINUE 700 DO 800 J = 1 , NFAS ESUM = 0.D0 DO 750 I = 1 , NCOMP IF( J.EQ.1 )THEN SOLIDS(I) = SSAVE(I) VBEST(NFAS*NCOMP+I) = SSAVE(I)/ZFEED(I) ENDIF INDX = I + (J-1)*NCOMP XVL(I,J) = VBEST(INDX) ESUM = ESUM + ZFEED(I)*XVL(I,J) 750 CONTINUE SFAS(J) = ESUM 800 CONTINUE IF( NDPRT.GT.0 )WRITE(NDPRT,9040)GMIN , GNORM RETURN 9010 FORMAT(' GDEM INITIATED : ',/, & ' OBJECT.FUN CF1 CF2 CF3 ') 9020 FORMAT(/,' PHASE REMOVED IN DSUB ',/) 9030 FORMAT(2X,D14.6,3X,4F9.4) 9040 FORMAT(/,' END DSUB GMIN =',D25.15,' GNORM =',D15.3/) END **==NTOPI.FOR SUBROUTINE NTOPI(NCOMP,IPOS,ZFRAC,GTOT,LOGFUG,GNORM,SOLIDS) C C SUBROUTINE NTOPI CALCULATES FUGACITY COEFFICIENTS, GIVEN PHASE COMPO- C SITION; IN ADDITION THE GRADIENT NORM (IN LN-K-FACTORS) IS CALCULATED C C INPUT: NCOMP: NO. OF COMPONENTS IN MIXTURE C ZFRAC: ARRAY OF MOLE FRACTIONS C SOLIDS: AMOUNTS OF SOLID PHASES C IPOS: VECTOR OF POINTER TO DEPENDENT PHASE C C OUTPUT : C GTOT: GIBBS ENERGY, RELATIVE TO INITIAL PHASE C LOGFUG: ARRAY OF LN FUGACITY COEFFICIENTS C GNORM: GRADIENT NORM C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION LOGFUG PARAMETER(MDIM=20,NFMAX=5) DIMENSION ZFRAC(MDIM,*) , LOGFUG(MDIM,*) , SOLIDS(*) , IPOS(*) DIMENSION X(MDIM) , Y(MDIM) LOGICAL EOSTYP COMMON /SOLCL / FUGAR(MDIM) , ISOLC(MDIM) , ISOLPC COMMON /FPROP / FUG(MDIM) , FUGX(MDIM,MDIM) COMMON /KEEP / NFAS , T , P , GNUL , ZFEED(MDIM) , FR(MDIM) , & XVL(MDIM,NFMAX) , SFAS(NFMAX) , ZVAL(NFMAX) , & IFLU(MDIM) , EOSTYP , ISTFAS , ISTAB GNORM = 0. GTOT = -GNUL C C SUM OVER PHASES C DO 200 J = 1 , NFAS GFAS = 0.D0 DO 50 I = 1 , NCOMP X(I) = ZFRAC(I,J) 50 CONTINUE CALL FUGMOD(NCOMP,1,IFLU(J),IC,T,P,ZVAL(J),X,FUG,FUGX) DO 100 I = 1 , NCOMP XI = X(I) XVL(I,J) = XI/ZFEED(I)*SFAS(J) IF( XI.GT.1.D-15 )THEN FAC = FUG(I) + DLOG(XI) GFAS = GFAS + XI*FAC ENDIF LOGFUG(I,J) = FUG(I) 100 CONTINUE GTOT = GTOT + GFAS*SFAS(J) 200 CONTINUE C C CHECK CURRENT DEVIATION C C C NORMALIZE FUGACTIY COEFFICIENTS C DO 300 I = 1 , NCOMP JPOS = IPOS(I) Y(I) = DLOG(ZFRAC(I,JPOS)+1.D-20) ARG = LOGFUG(I,JPOS) FUG(I) = ARG DO 250 J = 1 , NFAS EE = LOGFUG(I,J) - ARG LOGFUG(I,J) = EE IF( J.NE.JPOS .AND. ZFRAC(I,J).GT.1.D-15 )THEN DIFF = DLOG(ZFRAC(I,J)) + EE - Y(I) GNORM = GNORM + (ZFEED(I)*DIFF)**2 ENDIF 250 CONTINUE 300 CONTINUE C C LOOP SOLIDS C DO 400 I = 1 , NCOMP IF( ISOLC(I).EQ.0 )THEN LOGFUG(I,NFAS+1) = 0.D0 ELSE LOGFUG(I,NFAS+1) = FUGAR(I) - FUG(I) GTOT = GTOT + FUGAR(I)*SOLIDS(I) IF( SOLIDS(I).GT.0.D0 )GNORM = GNORM + & (ZFEED(I)*(LOGFUG(I,NFAS+1)-Y(I)))**2 ENDIF C ELSE IF( ISOLC(I).GT.0 .AND. Y(I).GT.FUGAR(I) )ISOLPC = 1 400 CONTINUE RETURN END **==FSPLIT.FOR SUBROUTINE FSPLIT(NCOM,NFAS,FEED,BETA,BETASL,FUGFLU,FUGSOL,XMBR, & MAXIT,ICODE,ISOLST,SOLIN) C C INPUT PARAMETERS ARE: C C NCOM: ACTUAL NUMBER OF COMPONENTS C NFAS: NUMBER OF MIXED PHASES C FEED: OVERALL MOLAR COMPOSITION OF FEED C BETA: INITIAL ESTIMATE OF PHASE AMOUNTS; AT LEAST ONE ELEMENT C MUST BE NONZERO C FUGFLU: ARRAY OF FLUID PHASE LN FUGACITY COEFFICIENTS C FUGSOL: VECTOR OF COMPONENT PURE PHASE FUGACITY COEFFICIENTS C SOLIN: LOGICAL VARIABLE, SET FOR SOLIDS FORMATION C C OUTPUT: C C BETASL: VECTOR OF PURE PHASE MOLAR AMOUNTS C XMBR: ARRAY OF FLUID PHASE MOLE FRACTIONS C MAXIT: LIMIT TO ITERATION COUNT C ICODE: STATUS INDICATOR ICODE=0: CONVERGED, ICODE=1: NOT CONVERGED C ISOLST: SOLIDS STATUS C IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C MAX.SIZE CAN BE MODIFIED BY CHANGING NFMAX AND MDIM C PARAMETER(NFMAX=5,MDIM=20) LOGICAL FIRST , SOLIN DIMENSION BETA(*) , BETASL(*) , FUGFLU(MDIM,*) , FUGSOL(*) , & FEED(*) DIMENSION COPFLU(MDIM,NFMAX) , ISOLST(*) , XMBR(MDIM,*) DIMENSION BETOLD(NFMAX) , DBETA(NFMAX) , HAR(MDIM,NFMAX) C C REORDER FUGACITY COEFFICIENTS AND CREATE INVERSE ARRAY C DO 200 I = 1 , NCOM C---GET PHASE WITH SMALLEST FUGACITY; THE FUGACITY COEFFICIENT FOR THIS C PHASE IS SET TO UNITY, AND THE REMAINING PHASES ARE ADJUSTED ACCOR- C DINGLY C IF( ISOLST(I).EQ.0 )THEN FGMIN = 100 ELSE FGMIN = FUGSOL(I) ENDIF BETASL(I) = 0.D0 DO 50 J = 1 , NFAS TT = FUGFLU(I,J) IF( TT.LT.FGMIN )FGMIN = TT 50 CONTINUE DO 100 J = 1 , NFAS TT = FGMIN - FUGFLU(I,J) IF( TT.GT.-50.D0 )THEN COPFLU(I,J) = DEXP(TT) ELSE COPFLU(I,J) = 0.D0 ENDIF 100 CONTINUE 200 CONTINUE C C PROBLEM IS SET UP; START OPTIMIZATION C CALL OBFUN(NCOM,NFAS,FEED,BETA,BETASL,COPFLU,ISOLST,OBVAL,SOLIN) XVAL = OBVAL ITER = 0 300 ITER = ITER + 1 C C CALCULATE CORRECTION VECTOR C CALL GRDHES(NCOM,NFAS,BETA,FEED,ISOLST,COPFLU,DBETA,GNORM,HAR, & SOLIN) C C CHECK CORRECTION FOR EXCEEDING FLUID PHASE BOUNDARY C CALL MAXUD(NFAS,BETA,DBETA,RFAK,LIMIT) FIRST = .TRUE. DO 400 J = 1 , NFAS BETOLD(J) = BETA(J) 400 CONTINUE 500 DO 600 J = 1 , NFAS BETA(J) = BETOLD(J) + RFAK*DBETA(J) 600 CONTINUE CALL OBFUN(NCOM,NFAS,FEED,BETA,BETASL,COPFLU,ISOLST,OBVAL,SOLIN) IF( OBVAL.GT.XVAL+1.D-10 )THEN C----CORRECTION FAILED FIRST = .FALSE. RFAK = RFAK/2 GOTO 500 ENDIF XVAL = OBVAL IF( LIMIT.NE.0 .AND. FIRST )BETA(LIMIT) = 0. IF( GNORM.GT.1.D-12 .AND. ITER.LT.MAXIT )GOTO 300 IF( ITER.EQ.MAXIT )THEN C---IT DID NOT CONVERGE; ICODE = 1 ELSE ICODE = 0 ENDIF C C CALCULATE PHASE COMPOSITION CORRESPONDING TO CURRENT STATE C DO 800 I = 1 , NCOM BSUM = BETASL(I) DO 650 J = 1 , NFAS BSUM = BSUM + BETA(J)*COPFLU(I,J) 650 CONTINUE DO 700 J = 1 , NFAS XMBR(I,J) = COPFLU(I,J)*FEED(I)/BSUM 700 CONTINUE 800 CONTINUE C C RENORMALIZE - TO BE ON THE SAFE SIDE C DO 1000 J = 1 , NFAS SUMMBR = 0.D0 DO 850 I = 1 , NCOM SUMMBR = SUMMBR + XMBR(I,J) 850 CONTINUE BETA(J) = BETA(J)*SUMMBR DO 900 I = 1 , NCOM XMBR(I,J) = XMBR(I,J)/SUMMBR 900 CONTINUE 1000 CONTINUE C C XMBR NOW CONTAINS NORMALIZED MOLE FRACTIONS, C WHICH,WHEN MULTIPLIED BY BETA SATISFIES THE C MATERIAL BALANCE C RETURN END **==OBFUN.FOR SUBROUTINE OBFUN(NCOM,NFAS,FEED,BETA,BETASL,COPFLU,ISOLST,OBVAL, & SOLIN) C C OBFUN CALCULATES OBJECTIVE FUNCTION FOR IDEAL SOLUTION CALCULATION C C INPUT: NCOM: NO. OF COMPONENTS C NFAS: NO. OF FLUID PHASES C FEED: OVERALL COMPOSITION C BETA: VECTOR OF FLUID PHASE AMOUNTS C BETASL: VECTOR OF SOLID PHASE AMOUNTS C COPFLU: LN FUG.COEFF. DIFFERENCES C ISOLST: SOLID STATUS INDICATOR (MAY BE MODIFIED) C SOLIN: SOLID FORMATION INDICATOR C C OUTPUT: OBVAL: VALUE OF OBJECTIVE FUNCTION C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20) DIMENSION BETA(*) , BETASL(*) , COPFLU(MDIM,*) , ISOLST(*) , & FEED(*) LOGICAL SOLIN , SOLOUT C C GIVEN THE BETA-S, CALCULATE SOLIDS AND OBJECTIVE FUNCTION C C START WITH SOLID DISPOSITION C SOLOUT = .NOT.SOLIN DO 100 I = 1 , NCOM IF( .NOT.(ISOLST(I).EQ.0 .OR. SOLOUT) )THEN C---FORMATION POSSIBLE SSM = 0. DO 20 J = 1 , NFAS SSM = SSM + BETA(J)*COPFLU(I,J) 20 CONTINUE DIF = FEED(I) - SSM IF( DIF.GT.0. )THEN C----PHASE FORMS BETASL(I) = DIF ISOLST(I) = 2 ELSE C----PHASE DID NOT FORM BETASL(I) = 0 ISOLST(I) = 1 ENDIF ENDIF 100 CONTINUE C C DISTRIBUTION CLEAR; EVALUATE OBJECTIVE FUNCTION C QFUN = 0. DO 200 I = 1 , NCOM IF( SOLIN .AND. ISOLST(I).EQ.2 )THEN TERM = DLOG(FEED(I)) - 1 ELSE TERM = 0. DO 120 J = 1 , NFAS TERM = TERM + BETA(J)*COPFLU(I,J) 120 CONTINUE TERM = DLOG(TERM) ENDIF QFUN = QFUN - FEED(I)*TERM 200 CONTINUE DO 300 J = 1 , NFAS TT = 1. DO 250 I = 1 , NCOM IF( ISOLST(I).EQ.2 .AND. SOLIN )TT = TT - COPFLU(I,J) 250 CONTINUE QFUN = QFUN + BETA(J)*TT 300 CONTINUE OBVAL = QFUN RETURN END **==GRDHES.FOR SUBROUTINE GRDHES(NCOM,NFAS,BETA,FEED,ISOLST,COPFLU,DBETA,GNORM, & HARRY,SOLIN) C C GRDHES CALCULATES - YES, YOU GUESSED IT - THE GRADIENT AND THE HESSIAN C MATRIX FOR THE IDEAL SOLUTION CALCULATION, AND CALCULATES THE CORREC- C TION TO THE PHASE AMOUNTS C C INPUT: NCOM: NO. OF COMPONENTS C NFAS: NO. OF FLUID PHASES C BETA: VECTOR OF FLUID PHASE AMOUNTS C FEED: OVERALL COMPOSITION C ISOLST: SOLID STATUS INDICATOR (MAY BE MODIFIED) C COPFLU: LN FUG.COEFF. DIFFERENCES C SOLIN: SOLID FORMATION INDICATOR C C OUTPUT: DBETA: CORRECTION VECTOR C GNORM: GRADIENT NORM C HARRY: HESSIAN MATRIX (WORK AREA) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(NFMAX=5,MDIM=20) LOGICAL SOLIN , SOLOUT DIMENSION INDEX(NFMAX) , IPIV(NFMAX) , GRDL(NFMAX) DIMENSION BETA(*) , COPFLU(MDIM,*) , ISOLST(*) , FEED(*) , & DBETA(*) DIMENSION GRAD(NFMAX) , HESS(NFMAX,NFMAX) , HARRY(MDIM,*) SOLOUT = .NOT.SOLIN C C ZERO INITIAL STATE C NEGINC = 1 C----INITIALIZE 100 DO 200 I = 1 , NFAS DBETA(I) = 0. DO 150 J = 1 , NFAS HESS(I,J) = 0. 150 CONTINUE 200 CONTINUE C-----CALCULATE ALL GRADIENT CONSTANT TERMS DO 300 J = 1 , NFAS GRD = 1 DO 250 I = 1 , NCOM C----CONTRIBUTION FROM PURE PHASES IF( ISOLST(I).EQ.2 .AND. SOLIN )GRD = GRD - COPFLU(I,J) 250 CONTINUE GRAD(J) = GRD 300 CONTINUE C C KAPPA-FACTORS C DO 400 I = 1 , NCOM IF( ISOLST(I).NE.2 .OR. SOLOUT )THEN TT = 0. DO 320 J = 1 , NFAS TT = TT + BETA(J)*COPFLU(I,J) 320 CONTINUE DO 340 J = 1 , NFAS HARRY(I,J) = COPFLU(I,J)/TT 340 CONTINUE ENDIF 400 CONTINUE C C GRADIENT CONTRIBUTION FROM FLUID PHASES C DO 500 I = 1 , NCOM IF( ISOLST(I).NE.2 .OR. SOLOUT )THEN DO 420 J = 1 , NFAS T1 = FEED(I)*HARRY(I,J) GRAD(J) = GRAD(J) - T1 420 CONTINUE ENDIF 500 CONTINUE MAXNEG = 0 VALNEG = 0 NTAL = 0 C C SELECT ACTIVE COMPONENTS, I.E. ALL PHASES WITH NONZERO BETAS C ONE PHASE WITH ZERO BETA-VALUE AND NEGATIVE GRADIENT CAN BE C INCLUDED C DO 600 J = 1 , NFAS C----SELECT ACTIVE PHASES IF( BETA(J).GT.0. )THEN NTAL = NTAL + 1 C----INDEX KEEPS TRACK ON ACTIVE PHASES INDEX(NTAL) = J GRDL(NTAL) = GRAD(J) C---AND CHECK INACTIVE FOR NEGATIVE GRADIENTS ELSEIF( GRAD(J).LT.VALNEG .AND. NEGINC.NE.0 )THEN VALNEG = GRAD(J) MAXNEG = J ENDIF 600 CONTINUE C C CHECK FOR NEGATIVE GRADIENT AMONG INACTIVE PHASES C IF( MAXNEG.NE.0 )THEN NTAL = NTAL + 1 INDEX(NTAL) = MAXNEG GRDL(NTAL) = GRAD(MAXNEG) ENDIF C------FORM APPROPRIATE HESSIAN DO 700 I = 1 , NCOM IF( ISOLST(I).NE.2 .OR. SOLOUT )THEN DO 620 J = 1 , NTAL JVAL = INDEX(J) T1 = HARRY(I,JVAL)*FEED(I) DO 610 K = 1 , J KVAL = INDEX(K) TT = T1*HARRY(I,KVAL) HESS(J,K) = HESS(J,K) + TT 610 CONTINUE 620 CONTINUE ENDIF 700 CONTINUE DO 800 J = 1 , NTAL DO 750 K = 1 , J HESS(K,J) = HESS(J,K) 750 CONTINUE C----ADD DIAGONAL TERM TO HANDLE SINGULAR HESSIAN HESS(J,J) = HESS(J,J) + 1.D-7 800 CONTINUE IF( NTAL.EQ.0 )RETURN CALL LU(NFMAX,NTAL,IPIV,HESS) CALL BACK(NFMAX,NTAL,IPIV,HESS,GRDL) IF( MAXNEG.NE.0 .AND. GRDL(NTAL).GT.0. )THEN C----CHANGE IN REACTIVATED PHASE GIVES NEGATIVE AMOUNT; RECALCU C----LATE WITHOUT THIS PHASE NEGINC = 0 GOTO 100 ENDIF C C SAVE CHANGES IN BETA C GNORM = 0. DO 900 J = 1 , NTAL DELB = GRDL(J) GNORM = GNORM + DELB**2 JVAL = INDEX(J) DBETA(JVAL) = -DELB 900 CONTINUE RETURN END **==MAXUD.FOR SUBROUTINE MAXUD(NFAS,BETA,DBETA,RFAK,LIMIT) C C MAXUD IS USED TO C DETERMINE REDUCTION FACTOR TO ENSURE NONNEGATIVITY OF BETA-VALUES C LIMITING CASE: EXACTLY ONE BETA BECOMES ZERO C C INPUT: NFAS: NO. OF PHASES C BETA: CURRENT PHASE FRACTIONS C DBETA: CORRECTION C C OUTPUT: RFAK: REDUCTION FACTOR C LIMIT: LIMITING ELEMENT IN DBETA C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION BETA(*) , DBETA(*) C C CHECK IF CORRECTION CAUSES NEGATIVE BETA C CALCULATE REDUCTION FACTOR, IF APPROPRIATE C C RFAK: CALCULATED STEPSIZE REDUCTION FACTOR C LIMIT: BETA-INDEX FOR PHASE WHICH BECOMES ZERO C RFAK = 1.D0 LIMIT = 0 DO 100 J = 1 , NFAS IF( DBETA(J).LT.0. )THEN BETN = BETA(J) + RFAK*DBETA(J) IF( BETN.LT.0. )THEN RFAK = -BETA(J)/DBETA(J) LIMIT = J ENDIF ENDIF 100 CONTINUE RETURN END **==PHASN.FOR SUBROUTINE PHASN(NARG,ICALC,FUN,GRAD,HESS,CMPFRC) C C PHASN EVALUATES FUNCTION, GRADIENT AND HESSIAN FOR PHASE SPLIT C CALCULATIONS. PHASN IS AN ARGUMENT OF THE MURRAY ROUTINE C C INPUT: C C NARG: NO. OF INDEPENDENT VARIABLES C ICALC: CALCULATION TYPE: C ICALC = 0: OBJ. FUNCTION ONLY C ICALC = 1: ALSO GRADIENT C ICALC = 2: ALSO HESSIAN C C OUTPUT: C C FUNC: OBJECTIVE FUNCTION C GRAD: GRADIENT C HESS: HESSIAN C CMPFRC: C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,NFMAX=5,NEQM=NFMAX*MDIM) DIMENSION INDEX(MDIM) , NCON(MDIM,NFMAX+1) , ENER(NFMAX+1) , & HESS(NARG,*) DIMENSION CMPFRC(*) , GRAD(*) , X(MDIM) , TM(MDIM) DIMENSION XMIN(MDIM) , DELEN(NFMAX+1) , FUGSAV(MDIM,NFMAX+1) DIMENSION IRPOS(NEQM) LOGICAL ACTIVE , NOTSOL , EOSTYP INTEGER DEPFAS COMMON /ERRDEB/ NDPRT COMMON /KEEP / NFAS , T , P , GNUL , ZZ(MDIM) , FR(MDIM) , & XVL(MDIM,NFMAX) , SFAS(NFMAX) , ZVAL(NFMAX) , & IFLU(MDIM) , EOSTYP , ISTFAS , ISTAB COMMON /FPROP / FUG(MDIM) , FUGX(MDIM,MDIM) COMMON /SOLCL / FUGAR(MDIM) , ISOL(MDIM) , ISTAT COMMON /OPTIM / VBEST(MDIM*(NFMAX+1)) , GMIN COMMON /AKTIV / ACTIVE(MDIM) , IPOS(NEQM) , NACT , EPSV SAVE INDEX , NCON , ENER , FUGSAV , XMIN SAVE IRPOS JD = 1 NCOMP = NARG/(NFAS-1) 100 NG = NFAS - 1 NFE = NFAS - ISTAT NARG = NCOMP*NG IF( ICALC.LT.2 )THEN C C CALLED FROM MURML TO EVALUATE CORRECTION C C REARRANGE VARIABLES TO NATURAL ORDER C CALL SWAP(NCOMP,NFAS,INDEX,CMPFRC) DO 200 I = 1 , NCOMP C---- DEPFAS IS POSITION OF DEPENDENT VARIABLE (CALCULATED BY MASSBALANCE) DEPFAS = INDEX(I) XLSUM = 1.D0 C--- GET AMOUNT OF DEPENDENT DO 120 J = 1 , NFAS IF( J.NE.DEPFAS )THEN EE = CMPFRC(I+(J-1)*NCOMP) XVL(I,J) = EE XLSUM = XLSUM - EE ENDIF 120 CONTINUE XVL(I,DEPFAS) = XLSUM C C CORRECT TRACE COMPONENTS, USING OLD K-FACTORS AND PHASE AMOUNTS C DO 140 J = 1 , NFAS IF( J.NE.DEPFAS )THEN IF( XVL(I,J).NE.0 .AND. NCON(I,J).EQ.-1 )THEN C---CORRECT TRACE RADIF = DEXP(FUGSAV(I,DEPFAS)-FUGSAV(I,J)) C---- DEPENDENT SOLID IF( DEPFAS.GT.NFE )THEN XVLN = SFAS(J)*RADIF/ZZ(I) ELSE C- - DEPENDENT FLUID XJ = XVL(I,DEPFAS)/SFAS(DEPFAS)*RADIF XVLN = XJ*SFAS(J) ENDIF DIF = XVL(I,J) - XVLN XVL(I,J) = XVLN XVL(I,DEPFAS) = XVL(I,DEPFAS) + DIF ENDIF ENDIF 140 CONTINUE SM = 0.D0 C---- CHECK FOR STEP MAKING CONCENTRATIONS NEGATIVE OR ZERO DO 160 J = 1 , NFAS SLIMT = 1.D-8 C - - SOLID MIGHT BE REMOVABLE IF( J.EQ.NFAS .AND. ISTAT.EQ.1 )SLIMT = 0. IF( XVL(I,J).LT.0.D0 )XVL(I,J) = SLIMT SM = SM + XVL(I,J) 160 CONTINUE DO 180 J = 1 , NFAS XVL(I,J) = XVL(I,J)/SM CMPFRC(I+(J-1)*NCOMP) = XVL(I,J) 180 CONTINUE C----- END CHECK 200 CONTINUE C C END FIRST DERIVATIVE MODIFICATION C ELSE C C SECOND DERIVATIVES NEEDED. C IDENTIFY TRACE COMPONENTS AND DEPENDENT VARIABLES C JD = 2 SMIN = 1. DO 250 J = 1 , NFAS IF( SFAS(J).LE.SMIN )THEN SMIN = SFAS(J) JMIN = J ENDIF 250 CONTINUE C C ARRAY NCON: 1: INDEPENDENT; 0: DEPENDENT; -1: TRACE C INDEX KEEPS TRACE OF DEPENDENT VARIABLES C DO 300 I = 1 , NCOMP C--- ACTIVE FLAGS WHEN ALL COMPONENTS ARE 'TRACE' COMPONENTS ACTIVE(I) = .FALSE. YMAX = 0 DO 260 J = 1 , NFAS IPLACE = I + (J-1)*NCOMP CMPFRC(IPLACE) = VBEST(IPLACE) XVL(I,J) = VBEST(IPLACE) XM = CMPFRC(IPLACE)*ZZ(I)/SFAS(J) C--- SAVE COMPOSITION OF PHASE PRESENT IN SMALLEST AMOUNT IF( J.EQ.JMIN )XMIN(I) = XM NCON(I,J) = 1 IF( XM.GT.EPSV )THEN ACTIVE(I) = .TRUE. NCON(I,J) = 1 ELSE NCON(I,J) = -1 ENDIF C--- IDENTIFY PHASE CONTAINING LARGEST PROPORTION OF COMPONENT IF( CMPFRC(IPLACE).GE.YMAX )THEN YMAX = CMPFRC(IPLACE) JMAX = J ENDIF 260 CONTINUE C---- INDEX(I): PHASE WHERE COMP. I IS PRESENT IN LARGEST AMOUNT INDEX(I) = JMAX C---- NCON SET ZERO FOR 'DEPENDENT' PHASE NCON(I,JMAX) = 0 IF( ACTIVE(I) )THEN DO 270 J = 1 , NFAS C--- SKIP LOOP, IF ACTIVE COMPONENTS ARE FOUND IF( NCON(I,J).GT.0 )GOTO 300 270 CONTINUE NCON(I,JMAX) = -1 ACTIVE(I) = .FALSE. ELSE NCON(I,JMAX) = -1 ENDIF 300 CONTINUE C C ORDER ACTUAL INDEPENDENT VARIABLES FOR MURML C CALL SWAP(NCOMP,NFAS,INDEX,CMPFRC) NACT = 0 IPOS(1) = 1 C--- GET LIST OF ACTIVE VARIABLES (DISCOUNT DEPENDENT AND TRACE) DO 350 J = 1 , NG JNUL = (J-1)*NCOMP DO 320 I = 1 , NCOMP IF( ACTIVE(I) )THEN C- - - - - POSITION OF VARIABLE IJ = I + JNUL C- -- - PHASE JJ = J C- - -- PHASE DEPENDENT; SWAP WITH LAST IF( J.EQ.INDEX(I) )JJ = NFAS IF( NCON(I,JJ).GE.0 )THEN NACT = NACT + 1 IPOS(NACT) = IJ ENDIF ENDIF 320 CONTINUE 350 CONTINUE IF( NACT.EQ.0 )NACT = 1 C C IRPOS-ARRAY LOCATES ACTIVE VARIABLES C IRPOS=0 MEANS INACTIVE C DO 400 I = 1 , NARG IRPOS(I) = 0 400 CONTINUE DO 450 I = 1 , NACT IJ = IPOS(I) IRPOS(IJ) = I 450 CONTINUE ENDIF C------INITIATE OBJECTIVE FUNCTION FUN = -GNUL IF( ICALC.NE.0 )THEN C------INITIATE GRADIENT DO 500 I = 1 , NACT GRAD(I) = 0. IF( ICALC.GE.2 )THEN C------INITIATE HESSIAN DO 460 J = 1 , NACT HESS(I,J) = 0. 460 CONTINUE ENDIF 500 CONTINUE ENDIF C C START CALCULATING CONTRIBUTION FROM INDIVIDUAL PHASES C DO 700 K = 1 , NFAS NOTSOL = (K.NE.NFAS .OR. ISTAT.NE.1) SFAS(K) = 0. C--- ENER IS THE ACTUAL PHASE GIBBS ENERGY ENER(K) = 0. DO 550 I = 1 , NCOMP X(I) = XVL(I,K)*ZZ(I) SFAS(K) = SFAS(K) + X(I) 550 CONTINUE IF( NOTSOL )CALL FUGMOD(NCOMP,JD,IFLU(K),IC,T,P,ZVAL(K),X,FUG, & FUGX) DELEN(K) = 0 DO 600 I = 1 , NCOMP IF( NOTSOL )THEN FUGSAV(I,K) = FUG(I) C--- ADDITION OF SMALL TERM TO GUARD AGAINST ZERO CONC. XE = X(I)/SFAS(K) + 1.D-20 TMI = FUG(I) + DLOG(XE) ELSE IF( ISOL(I).LT.2 )FUGSAV(I,K) = 100. IF( ISOL(I).EQ.2 )FUGSAV(I,K) = FUGAR(I) TMI = 0 IF( ISOL(I).EQ.2 )TMI = FUGAR(I) ENDIF TM(I) = TMI DELEN(K) = DELEN(K) + XMIN(I)*TMI ENER(K) = ENER(K) + TMI*X(I) 600 CONTINUE C C ADD CURRENT CONTRIBUTION TO OBJECTIVE FUNCTION C FUN = FUN + ENER(K) C C END GIBBS ENERGY, START GRADIENT C IF( ICALC.NE.0 )THEN DO 620 I = 1 , NCOMP S = TM(I)*ZZ(I) IPL = K DEPFAS = INDEX(I) IF( K.EQ.NFAS )IPL = DEPFAS IF( K.EQ.DEPFAS )THEN DO 605 M = 1 , NG IPS = I + (M-1)*NCOMP IR = IRPOS(IPS) IF( IR.GT.0 )GRAD(IR) = GRAD(IR) - S 605 CONTINUE ELSE IPS = I + (IPL-1)*NCOMP IR = IRPOS(IPS) IF( IR.GT.0 )GRAD(IR) = GRAD(IR) + S ENDIF 620 CONTINUE C C END GRADIENT, START HESSIAN MATRIX C IF( ICALC.GE.2 )THEN IF( NOTSOL )THEN DO 625 I = 1 , NCOMP S = ZZ(I)/SFAS(K) DO 622 J = 1 , I FUGX(I,J) = S*(FUGX(I,J)-1.)*ZZ(J) FUGX(J,I) = FUGX(I,J) 622 CONTINUE FUGX(I,I) = FUGX(I,I) + ZZ(I)/(XVL(I,K)+1.D-20) 625 CONTINUE DO 635 I = 1 , NCOMP DEPFAS = INDEX(I) IPL = K IF( K.EQ.NFAS )IPL = DEPFAS DO 630 J = 1 , NCOMP JV = INDEX(J) JPL = K IF( K.EQ.NFAS )JPL = JV S = FUGX(I,J) M11 = IPL M12 = IPL M21 = JPL M22 = JPL IF( DEPFAS.EQ.K )THEN S = -S M11 = 1 M12 = NG ENDIF IF( JV.EQ.K )THEN M21 = 1 M22 = NG S = -S ENDIF DO 628 M1 = M11 , M12 DO 626 M2 = M21 , M22 IPS = I + (M1-1)*NCOMP JPS = J + (M2-1)*NCOMP IR = IRPOS(IPS) IF( IR.NE.0 )THEN JR = IRPOS(JPS) IF( JR.NE.0 )HESS(IR,JR) = HESS(IR,JR) & + S ENDIF 626 CONTINUE 628 CONTINUE 630 CONTINUE 635 CONTINUE ENDIF ENDIF ENDIF 700 CONTINUE C C CHANGE BACK TO PERTURBED FORM FOR FIRST DERIVATIVES C IF( ICALC.EQ.1 )CALL SWAP(NCOMP,NFAS,INDEX,CMPFRC) IF( ICALC.GE.2 .AND. NFAS.GT.2 )THEN C C CHECK FOR REMOVAL OF A PHASE C DMIN = 1. DO 750 J = 1 , NFAS IF( J.NE.JMIN )THEN DELT = DELEN(J) - DELEN(JMIN) IF( DELT.LE.DMIN )THEN K = J DMIN = DELT ENDIF ENDIF 750 CONTINUE IF( DMIN.NE.1. )THEN J = JMIN DO 760 I = 1 , NCOMP X(I) = ZZ(I)*(XVL(I,J)+XVL(I,K)) 760 CONTINUE SFT = SFAS(J) + SFAS(K) CALL FUGMOD(NCOMP,1,IFLU(K),IC,T,P,ZNEW,X,FUG,FUGX) ENGO = ENER(J) + ENER(K) DO 780 I = 1 , NCOMP XE = X(I)/SFT + 1.D-20 ENGO = ENGO - X(I)*(FUG(I)+DLOG(XE)) 780 CONTINUE IF( ENGO.GE.-1.D-10 )THEN IF( J.EQ.NFAS )ISTAT = 0 C C COMBINE PHASES AND RESET FLAGS C IF( K.GT.J )THEN IFLUK = IFLU(J) IFLU(J) = IFLU(K) IFLU(K) = IFLUK JJJ = J J = K K = JJJ ENDIF C- - - - PHASE TO REMOVE NOW HAS THE HIGHER INDEX C C STORE BEST RESULT SO FAR IN VBEST C DO 790 I = 1 , NCOMP SS = XVL(I,K) + XVL(I,J) XVL(I,K) = SS XVL(I,J) = XVL(I,NFE) VBEST(I+(K-1)*NCOMP) = SS VBEST(I+(J-1)*NCOMP) = VBEST(I+(NFE-1)*NCOMP) C ZFRAC(I,NFAS) = XMIN(I) 790 CONTINUE IFLUJ = IFLU(J) IFLU(J) = IFLU(NFE) IFLU(NFE) = IFLUJ NFAS = NFAS - 1 IF( NDPRT.GT.0 )WRITE(NDPRT,9010) GOTO 100 ENDIF ENDIF ENDIF IF( FUN.LE.GMIN+1.D-10 )THEN C C NO REDUCTION, STORE BEST RESULT C GMIN = FUN DO 850 I = 1 , NCOMP DO 800 J = 1 , NFAS IPLACE = I + (J-1)*NCOMP VBEST(IPLACE) = XVL(I,J) 800 CONTINUE 850 CONTINUE ENDIF RETURN 9010 FORMAT(/,' NUMBER OF PHASES REDUCED ') END SUBROUTINE MAXFRC(NCOMP,NFAS,ARRAY,IPLACE) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20) C C SUBROUTINE RETURNS THE INDEX (IPLACE) C FOR THE LARGEST ELEMENT IN EACH ROW C OF THE MATRIX (NCOMP*NFAS) ARRAY C DIMENSION ARRAY(MDIM,*) , IPLACE(*) DO 100 I = 1 , NCOMP IVAL = 1 VALUE = ARRAY(I,1) DO 50 J = 2 , NFAS EE = ARRAY(I,J) IF( EE.GT.VALUE )THEN IVAL = J VALUE = EE ENDIF 50 CONTINUE IPLACE(I) = IVAL 100 CONTINUE RETURN END **==MURML.FOR SUBROUTINE MURML(N,IPR,NMAX,XLAM,GLIM,F,X,GD,G,W,FUNC) C C SECOND ORDER MINIMIZER USING MURRAY'S METHOD C C N: ACTUAL NO. OF INDEPENDENT VARIABLES (INPUT) C IPR: PRINT INDICATOR FOR TEST OUTPUT (INPUT) C NMAX: MAX NO. OF ITERATIONS (INPUT) C XLAM: UPPER BOUND ON STEPSIZE (INPUT) C GLIM: GRADIENT NORM CONVERGENCE CRITERION (INPUT) C X: VECTOR OF INDEPENDENT VARIABLES` (IN/OUT) C F: OBJECTIVE FUNCTION (OUTPUT) C GD: GRADIENT VECTOR C G: HESSIAN MATRIX C W: WORKSPACE FOR ROUTINE C FUNC: SUBROUTINE FOR EVALUATION OF OBJECTIVE FUNCTION C C CALLED FROM: MAIN C C CALLS: SPLIT,LINE,FUNC C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,NFMAX=5,NEQM=(NFMAX-1)*MDIM) EXTERNAL FUNC LOGICAL ACTIVE DIMENSION X(*) , GD(*) , G(N,*) , W(5,*) COMMON /AKTIV / ACTIVE(MDIM) , IVAR(NEQM+MDIM) , NACT , EPSV COMMON /ERRDEB/ NDPRT NTM = 0 C----DEL MUST BE APPROXIMATELY THE SQUARE ROOT OF MACHINE ACCURACY DEL = 1.D-7 100 NTM = NTM + 1 GPS = 0. C----INITIAL VALUE OF OBJECTIVE FUNCTION CALL FUNC(N,2,F,GD,G,X) DO 200 I = 1 , N W(4,I) = X(I) 200 CONTINUE DO 300 I = 1 , NACT W(5,I) = GD(I) GPS = GPS + GD(I)**2 300 CONTINUE C-----CHECK GRADIENT NORM IF( GPS.LT.GLIM )GOTO 1400 BETA = 0. C-----PICK SUBSET OF ACTIVE VARIABLES DO 400 I = 1 , NACT W(1,I) = 0. DO 350 J = 1 , I S = DABS(G(I,J)) IF( I.EQ.J )S = S*NACT IF( S.GT.BETA )BETA = S 350 CONTINUE 400 CONTINUE BETA = DSQRT(BETA/NACT) C----AUX. ROUTINE TO FACTORIZE CALL SPLIT(N,NACT,IDI,BETA,DEL,G,W) NTS = 0 500 NTS = NTS + 1 C C INNER ITERATIONS WITH CONSTANT HESSIAN C C-------CALCULATE STEP; RETURNED IN W(*,3) CALL LINE(N,NACT,GD,G,W) SMAX = 0. GP = 0. DP = 0. GPS = 0. DO 600 I = 1 , NACT S = W(3,I) IF( DABS(S).GT.SMAX )SMAX = DABS(S) GP = GP + S*GD(I) GPS = GPS + GD(I)**2 DP = DP + S*S*W(1,I) 600 CONTINUE FAC = 1. IF( SMAX.GT.XLAM )FAC = XLAM/SMAX C-----FAC IS REDUCTION FACTOR,WHEN MAX. STEP IS EXCEEDED ALFA = 1. 700 FF = ALFA*FAC DO 800 I = 1 , N X(I) = W(4,I) 800 CONTINUE DO 900 I = 1 , NACT IV = IVAR(I) X(IV) = X(IV) + FF*W(3,I) 900 CONTINUE IF( DABS(ALFA).LT.1.D-4 )ALFA = -ALFA IF( DABS(ALFA).GE.1.D-5 )THEN CALL FUNC(N,1,FNEW,GD,G,X) DELS = FNEW - F IF( DELS.GT.1.D-10 )THEN C----INCREASE IN OBJECTIVE FUNCTION; REDO STEP IF( NTS.GT.1 )GOTO 1200 ALFA = ALFA/3 GOTO 700 ENDIF ENDIF PRED = FF*GP - FF**2/2*(GP+DP) C----RATIO OF CALCUATED GAIN/PREDICTED GAIN CALPRE = DELS/(PRED+1.D-19) GNW = 0. DO 1000 I = 1 , NACT GNW = GNW + GD(I)**2 1000 CONTINUE RG = GNW/GPS GPS = GNW IF( GPS.LT.GLIM )GOTO 1400 F = FNEW DO 1100 I = 1 , N W(4,I) = X(I) 1100 CONTINUE IF( NTS.LT.5 )THEN C---------SUFFICIENT IMPROVEMENT TO CONTINUE ? IF( RG.LT..02 .AND. CALPRE.GT..3 )GOTO 500 ENDIF 1200 DO 1300 I = 1 , N X(I) = W(4,I) 1300 CONTINUE IF( IPR.NE.0 .AND. NDPRT.GT.0 )WRITE(NDPRT,9010)NTM , F , CALPRE , & GPS , ALFA IF( IPR.GT.1 .AND. NDPRT.GT.0 )WRITE(NDPRT,9020)(X(I),I=1,N) IF( IPR.GT.1 .AND. NDPRT.GT.0 )WRITE(NDPRT,9030) IF( NTM.LT.NMAX .AND. GPS.GT.GLIM )GOTO 100 1400 IF( NDPRT.GT.0 )WRITE(NDPRT,9040)NTM , GPS , F XLAM = GPS/GLIM IF( NDPRT.GT.0 )WRITE(NDPRT,9020)(X(I),I=1,N) RETURN 9010 FORMAT(' IT ',I3,' FNEW ',D12.5,' CP ',F8.3,' GNORM ',D9.2, & ' ALFA ',D10.2) 9020 FORMAT(/,' X-VECTOR: ',/,(4X,8F9.5)) 9030 FORMAT(' ') 9040 FORMAT(/,' ITEND ,',I4,' ITER, GNORM = ',D9.2,' OBJF ',D15.8/) END **==LINE.FOR SUBROUTINE LINE(N,NACT,GD,G,W) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C LINE IS A SERVICE ROUTINE FOR MURML C LINE CALCULATES THE BACKSUBSTITUTION C C N: TOTAL NO. OF INDEPENDENT VARIABLES C NACT: ACTUAL SIZE C GD: GRADIENT VECTOR (RHS OF EQN) C G: HESSIAN MATRIX C W: WORK AREA ; W(3,*) RETURNS CORRECTIONS C C CALLED FROM: MURML C DIMENSION GD(*) , G(N,*) , W(5,*) C----COPY GRADIENT ACCOREDING TO ACTIVE ELEMENTS (W5) DO 100 I = 1 , NACT W(5,I) = GD(I) 100 CONTINUE C----- PREPARED FOR MARQUARDT-LIKE VARIANT,BUT HERE XLAM=0 DO 200 I = 1 , NACT S = -W(5,I) DO 150 J = 1 , I - 1 S = S - G(I,J)*W(2,J) 150 CONTINUE W(2,I) = S/G(I,I) 200 CONTINUE DO 300 I = NACT , 1 , -1 S = W(2,I) DO 250 J = I + 1 , NACT S = S - G(J,I)*W(3,J) 250 CONTINUE W(3,I) = S/G(I,I) 300 CONTINUE RETURN END **==SPLIT.FOR SUBROUTINE SPLIT(N,NACT,IDI,BETA,DEL,G,W) C C SPLIT IS A SERVICE ROUTINE FOR MURML C SPLIT CALCULATES THE FACTORIZATION OF THE HESSIAN C C N: TOTAL NO. OF INDEPENDENT VARIABLES C NACT: ACTUAL SIZE C IDI: COUNTS DIAGONAL CORRECTION C BETA: TOLERANCE VARIABLE C DEL: TOLERANCE FOR POSITIVE DEFINITENESS C G: HESSIAN MATRIX C W: WORK AREA C C CALLED FROM: MURML C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION G(N,*) , W(5,*) IDI = 0 DO 200 I = 1 , NACT ID = 0 TV = G(I,I) SV = TV DO 50 J = 1 , I - 1 SV = SV - G(I,J)**2 50 CONTINUE IF( SV.LT.DEL*DEL )ID = 1 SVR = DSQRT(DABS(SV)) IF( SVR.LT.DEL )SVR = DEL XM = 0. DO 100 J = I + 1 , NACT S = G(J,I) DO 60 K = 1 , I - 1 S = S - G(I,K)*G(J,K) 60 CONTINUE S = S/SVR IF( DABS(S).GT.XM )XM = DABS(S) G(J,I) = S 100 CONTINUE IF( XM.GE.BETA )THEN ID = 1 XM = XM/BETA DO 120 J = I , NACT G(J,I) = G(J,I)/XM 120 CONTINUE SVR = SVR*XM ENDIF IF( ID.EQ.1 )W(1,I) = SVR**2 - SV G(I,I) = SVR IDI = IDI + ID 200 CONTINUE RETURN END **==LU.FOR SUBROUTINE LU(ND,N,IPIV,A) C C GAUSSIAN ELIMINATION ROUTINE; LU IS DECOMPOSITION STEP C C ND: ROW DIMENSION OF ARRAY TO BE FACTORIZED C N: ACTUAL SIZE OF MATRIX C IPIV: ROW PIVOT VECTOR C A: MATRIX TO FACTORIZE; RETURNS FACTORIZED MATRIX C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION IPIV(*) , A(ND,*) IPIV(N) = N DO 200 I = 1 , N - 1 X = DABS(A(I,I)) IPIV(I) = I DO 50 J = I + 1 , N Y = DABS(A(J,I)) IF( Y.GT.X )THEN X = Y IPIV(I) = J ENDIF 50 CONTINUE IF( IPIV(I).NE.I )THEN K = IPIV(I) DO 60 J = I , N X = A(I,J) A(I,J) = A(K,J) A(K,J) = X 60 CONTINUE ENDIF DO 100 J = I + 1 , N X = -A(J,I)/A(I,I) A(J,I) = X DO 80 K = I + 1 , N A(J,K) = A(J,K) + X*A(I,K) 80 CONTINUE 100 CONTINUE 200 CONTINUE RETURN END **==BACK.FOR SUBROUTINE BACK(ND,N,IPIV,A,V) C C GAUSSIAN ELIMINATION ROUTINE; BACK IS BACKSUBSTITUTION STEP C C ND: ROW DIMENSION OF ARRAY TO BE FACTORIZED C N: ACTUAL SIZE OF MATRIX C IPIV: ROW PIVOT VECTOR C A: FACTORIZED MATRIX C V: RHS=VECTOR; RETURNS RESULT C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION IPIV(*) , A(ND,*) , V(*) DO 100 I = 1 , N - 1 K = IPIV(I) IF( K.NE.I )THEN X = V(I) V(I) = V(K) V(K) = X ENDIF VI = V(I) DO 50 J = I + 1 , N V(J) = V(J) + A(J,I)*VI 50 CONTINUE 100 CONTINUE V(N) = V(N)/A(N,N) DO 200 I = N - 1 , 1 , -1 VI = V(I) DO 150 J = I + 1 , N VI = VI - A(I,J)*V(J) 150 CONTINUE V(I) = VI/A(I,I) 200 CONTINUE RETURN END **==SWAP.FOR SUBROUTINE SWAP(NCOMP,NFAS,INDEX,YVAL) C C SWAP IS USED FROM PHASN TO REARRANGE ELEMENTS IN YVAL C C NCOMP: NO. OF COMPONENTS C NFAS: NO. OF PHASES C INDEX: ARRAY GIVING LOCATION OF 'DEPENDENT' VARIABLE C YVAL: VECTOR OF INDEPENDENT- AND DEPENDENT VARIABLES C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION INDEX(*) , YVAL(*) DO 100 I = 1 , NCOMP IV = INDEX(I) IC1 = I + (IV-1)*NCOMP IC2 = I + (NFAS-1)*NCOMP S = YVAL(IC1) YVAL(IC1) = YVAL(IC2) YVAL(IC2) = S 100 CONTINUE RETURN END **==COPAR.FOR SUBROUTINE COPAR(NR,NC,NADIM,AR1,AR2) C C SUBROUTINE COPIES ARRAY AR1 INTO ARRAY AR2 C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION AR1(NADIM,*) , AR2(NADIM,*) DO 100 I = 1 , NR DO 50 J = 1 , NC AR2(I,J) = AR1(I,J) 50 CONTINUE 100 CONTINUE RETURN END **==ONEDIR.FOR SUBROUTINE ONEDIR(NCOMP,ZFRAC,RESV,YBEST) C C PURPOSE OF ROUTINE: TO LOCATE A NEAR-CRTICAL MINIMUM IN THE TPD. C C INPUT: C NCOMP: NO. OF COMPONENTS C ZFRAC: CURRENT MOLE FRACTIONS IN EACH PHASE C C OUTPUT: C RESV : MIN. OF OBJECTIVE FUNCTION C YBEST: LOG MOLE FRACTIONS AT MINIMUM C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,NFMAX=5) LOGICAL FIRST , POS , NEG DIMENSION YBEST(*) , ZFRAC(MDIM,*) DIMENSION Z(MDIM) , ZR(MDIM) , ZL(MDIM) , IP(MDIM) , YV(MDIM) DIMENSION YD(MDIM) , F1(MDIM) , FP(MDIM) DIMENSION YNEW(MDIM) , CORR(MDIM) , DIST(NFMAX) DIMENSION YVAL(MDIM) , GRD(MDIM) , GDG(MDIM) LOGICAL EOSTYP COMMON /ERRDEB/ NDPRT COMMON /KEEP / NFAS , T , P , GNUL , ZFEED(MDIM) , FR(MDIM) , & XVL(MDIM,NFMAX) , SFAS(NFMAX) , ZVAL(NFMAX) , & IFLU(MDIM) , EOSTYP , ISTFAS , ISTAB COMMON /FPROP / FUG(MDIM) , FUGX(MDIM,MDIM) SFACT = 1.30 RESV = 0. DO 1200 KKK = 1 , NFAS C - - -SKIP VAPOUR PHASE WITH HYBRID MODELS IF( IFLU(KKK).LT.0 )GOTO 1200 C C STARTING POINT: PHASES ALLREADY PRESENT C DO 50 I = 1 , NCOMP Z(I) = ZFRAC(I,KKK) + 1.D-25 ZR(I) = DSQRT(Z(I)) ZL(I) = DLOG(Z(I)) 50 CONTINUE C C SECOND DERIVATIVE MATRIX FOR PHASE; SCALE BY SQUARE ROOTS C CALL FUGMOD(NCOMP,2,0,IC,T,P,ZCOMP,Z,FUG,FUGX) DO 100 I = 1 , NCOMP CORR(I) = 0.D0 F1(I) = ZL(I) + FUG(I) RVAL = ZR(I) DO 60 J = I , NCOMP FUGX(I,J) = RVAL*FUGX(I,J)*ZR(J) FUGX(J,I) = FUGX(I,J) 60 CONTINUE 100 CONTINUE XMAX = 0.D0 DO 150 I = 1 , NCOMP XMT = DABS(FUGX(I,I)) IF( XMT.GT.XMAX )THEN IMAX = I XMAX = XMT ENDIF 150 CONTINUE C C EIGENVECTOR-EIGENVALUE ESTIMATE BY DIRECT ITERATION C SSQ = 0. DO 200 I = 1 , NCOMP YD(I) = FUGX(I,IMAX) SSQ = SSQ + YD(I)**2 200 CONTINUE SSQS = DSQRT(SSQ) DO 250 I = 1 , NCOMP YD(I) = YD(I)/SSQS 250 CONTINUE DO 350 M = 1 , 3 PR1 = 0. PR2 = 0. DO 280 I = 1 , NCOMP TERM = 0. DO 260 J = 1 , NCOMP TERM = TERM + FUGX(I,J)*YD(J) 260 CONTINUE PR1 = PR1 + TERM*YD(I) PR2 = PR2 + TERM*TERM YV(I) = TERM 280 CONTINUE PR2S = DSQRT(PR2) DO 300 I = 1 , NCOMP YD(I) = YV(I)/PR2S 300 CONTINUE 350 CONTINUE RAT = PR2/PR1 RAT1 = RAT + 1 IF( NDPRT.GT.0 )WRITE(NDPRT,9010)KKK , RAT1 C C EIGENVALUE TOO HIGH MEANS PHASE IS NOT NEAR CRITICAL C IF( RAT.GT.-.5 )GOTO 1200 C C OTHERWISE REFINE BY INVERSE ITERATION C DO 400 I = 1 , NCOMP FUGX(I,I) = FUGX(I,I) + 1 400 CONTINUE C---FACTORIZATION CALL LU(MDIM,NCOMP,IP,FUGX) RATO = 0. DO 500 K = 1 , 5 DO 420 I = 1 , NCOMP YV(I) = YD(I) 420 CONTINUE CALL BACK(MDIM,NCOMP,IP,FUGX,YV) SP1 = 0. SP2 = 0. DO 440 I = 1 , NCOMP YV(I) = YV(I) - YD(I) SP1 = SP1 + YV(I)*YD(I) SP2 = SP2 + YV(I)**2 440 CONTINUE RAT = SP2/SP1 RTEST = 1/(RAT+1) SP2 = DSQRT(SP2) DO 460 I = 1 , NCOMP YD(I) = YV(I)/SP2 460 CONTINUE IF( DABS(RTEST-RATO).LT.1.D-6 )GOTO 550 RATO = RTEST 500 CONTINUE 550 RAT = 1/(RAT+1) DO 600 I = 1 , NCOMP YV(I) = YD(I)/ZR(I)*.5D0 600 CONTINUE C C START LINESEARCH C POS = .FALSE. NEG = .FALSE. FUNMIN = 0. GYOLD = 1. NCN = 0 SMAX = 1.50D0 SNUL = .05D0 C C STARTING DISTANCE MAY MISS VERY CLOSE MINIMA C DS = SFACT*SNUL FIRST = .TRUE. IST = 0 S = SNUL LOC = 0 GOTO 700 650 GYOLD = GY IF( S.EQ.SMAX )GOTO 1150 S = S + DS IF( S.GT.SMAX )S = SMAX FIRST = .FALSE. IF( S.EQ.0 )GOTO 1150 LOC = 0 C-----TO CHECK WHEN CALCULATION AND REFINEMENT ARE OF OPPOSITE SIGN 700 LOC = LOC + 1 IF( FIRST )IST = IST + 1 DO 750 I = 1 , NCOMP SY = S*YV(I) IF( SY.GT.0. )THEN YNEW(I) = (ZR(I)*(1+SY)) ELSE FAC = (1-SY)*(1+SY*SY) YNEW(I) = ZR(I)/FAC ENDIF 750 CONTINUE DO 800 I = 1 , NCOMP IF( .NOT.FIRST )THEN C------CORR EVALUATES A QUADRATIC CORRECTION TO SEARCH DIRECTION TTT = S*S*CORR(I) YN = YNEW(I) + TTT YNEW(I) = YN C------REPLACE VERY SMALL QUANTITIES,USING DIRECT SUBSTITUTION IF( YN.LE.0. .OR. DABS(YN/TTT).LT..5 )YNEW(I) & = DEXP((F1(I)-FP(I))/2) ENDIF YVAL(I) = YNEW(I)**2 800 CONTINUE C C LOCATE CLOSEST PHASE C DO 850 L = 1 , NFAS DIS = 0. DO 820 I = 1 , NCOMP DIS = DIS + (YVAL(I)-ZFRAC(I,L))**2 820 CONTINUE DIST(L) = DIS 850 CONTINUE DO 900 L = 1 , NFAS IF( IFLU(L).GE.0 )THEN IF( DIST(L)/DIST(KKK).LT..5 )SMAX = S ENDIF 900 CONTINUE CALL FUGMOD(NCOMP,1,ISTAB,IC,T,P,ZCOMP,YVAL,FP,FUGX) C C NEW TPD IN FUN C XPR = 0. XPR2 = 0. FUN = 1. DO 950 I = 1 , NCOMP DDS = F1(I) - FP(I) IF( DDS.LT.-15. )THEN C C CORRECT BY DIRECT SUBSTITUTION C YNEW(I) = DEXP(DDS/2) YVAL(I) = YNEW(I)**2 TM = 0 ELSE TM = DLOG(YVAL(I)) - DDS ENDIF FUN = FUN + YVAL(I)*(TM-1) TM = DLOG(YVAL(I)) + FP(I) - F1(I) QQ = YNEW(I)*TM GRD(I) = QQ GDG(I) = QQ XPR = XPR + QQ*YD(I) XPR2= XPR2+ QQ*(YD(I) + 4.*S*CORR(I)) 950 CONTINUE GY = XPR2/S/RAT TN = (GY-1)*RAT/S IF( FIRST .AND. TN*S.GT.0. .AND. IST.LT.2 )THEN C C TEST FOR SIGN CHANGE (THIRD DERIV.) ON FIRST TRY C DO 960 I = 1 , NCOMP YV(I) = -YV(I) YD(I) = -YD(I) 960 CONTINUE LOC = 0 GOTO 700 ENDIF DO 1000 I = 1 , NCOMP GDG(I) = GDG(I) - XPR*YD(I) 1000 CONTINUE CALL BACK(MDIM,NCOMP,IP,FUGX,GDG) C C CORRECTION TO OBJECTIVE FUNCTION C GAIN = 0. DO 1050 I = 1 , NCOMP GAIN = GAIN + GDG(I)*(GRD(I)+XPR*YD(I)) 1050 CONTINUE GAIN = GAIN/2 DO 1100 I = 1 , NCOMP IF( FIRST )THEN CORR(I) = -GDG(I)/SNUL**2/2 ELSE CORR(I) = CORR(I) - GDG(I)/S/S/2 ENDIF 1100 CONTINUE FUNM = FUN - GAIN IF( LOC.EQ.1 .AND. FUNM*FUN.LT.0. )GOTO 700 IF( FUN.LT.FUNMIN )THEN FUNMIN = FUN DO 1120 I = 1 , NCOMP YBEST(I) = YVAL(I) 1120 CONTINUE ENDIF C------GY IS THE PROCECTED GRADIENT,WHICH DETERMINES THE SEARCH IF( RAT.LT.0. )GY = -GY IF( .NOT.(FIRST) )THEN IF( NEG )THEN IF( .NOT.(GY.LT.0. .AND. (.NOT.POS)) )THEN C C BRACKETED MINIMUM C POS = .TRUE. NCN = NCN + 1 IF( NCN.LE.3 )THEN IF( GY.LE.0. )THEN GYNEG = GY SNEG = S ELSE GYPOS = GY SPOS = S ENDIF S = SNEG - GYNEG*(SPOS-SNEG)/(GYPOS-GYNEG) DS = 0 GOTO 650 ENDIF GOTO 1150 ENDIF ELSEIF( GY.GT.0. )THEN C C TEST FOR INCREASING RATIO; TERMINATE IF FOUND C IF( GY.GT.GYOLD )GOTO 1150 C C NOTHING INTERESTING, BUT STILL HOPE; GO ON. C DS = SFACT*DS GOTO 650 ENDIF NEG = .TRUE. GYNEG = GY SNEG = S DS = SFACT*DS ENDIF GOTO 650 1150 IF( FUNMIN.LT.0. )THEN IF( NDPRT.GT.0 )WRITE(NDPRT,9020)FUNMIN RESV = FUNMIN IF( FUNMIN.LE.-1.D-7 )RETURN ENDIF 1200 CONTINUE RETURN 9010 FORMAT(/,' PHASE INVESTIGATED ',I3,' EIGENVALUEESTIM.',D12.4) 9020 FORMAT(/,' SEARCH CONCLUSIVE, FMIN = ',D10.2) END **==PURE.FOR C C SUBROUTINE PURE(STRTUP,NCOMP,SELECT,ZFRAC,Y,TPDMIN) C C STRTUP: FIRST CALL OR NOT - INDICATOR C NCOMP: NO. OF COMPONENTS IN MIXTURE (INPUT) C ZFRAC: CURRENT MOLE FRACTIONS IN EACH PHASE C Y: TRIAL PHASE COMPOSITION (BEST) (OUTPUT) C TPDMIN: MIN. VALUE OF TANGENT PLANE C DISTANCE IN SEARCH (OUTPUT) C C CALLED FROM: MAIN C C CALLS EOSPUR C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,NFMAX=5) LOGICAL STRTUP , SELECT(*) DIMENSION Y(*) , ZFRAC(MDIM,*) DIMENSION LIST(MDIM+2) , RTOLD(MDIM+2) , OLDY(MDIM+2,MDIM) DIMENSION YTEST(MDIM) , DS(MDIM) , DSG(MDIM) LOGICAL EOSTYP C C LOCAL VARIABLES LIST, OLDY AND RTOLD MUST SURVIVE CALL C COMMON /KEEP / NFAS , T , P , GNUL , ZFEED(MDIM) , FR(MDIM) , & XVL(MDIM,NFMAX) , SFAS(NFMAX) , ZVAL(NFMAX) , & IFLU(MDIM) , EOSTYP , ISTFAS , ISTAB COMMON /FPROP / FUG(MDIM) , FUGX(MDIM,MDIM) COMMON /ERRDEB/ NDPRT SAVE LIST , RTOLD , OLDY ISTAB = 0 NA = 2 IF( .NOT.EOSTYP )NA = 3 TPDMIN = 1. JPUR = 0 JMIN = 0 KMIN = 1 NT = NCOMP + 1 C C LIST IS THE SEARCH STATUS INDICATOR C DO 100 I = 1 , NT IF( STRTUP )LIST(I) = 1 IF( LIST(I).NE.0 )LIST(I) = 1 100 CONTINUE IF( NFAS.NE.1 .AND. EOSTYP )THEN NT = NCOMP + 2 ZMAX = 0. JMAX = 0 C C GET PHASE PRESENT IN LARGEST AMOUNT C DO 150 J = 1 , NFAS IF( ZVAL(J).GT.ZMAX )THEN ZMAX = ZVAL(J) JMAX = J ENDIF 150 CONTINUE C C CALCULATE EXTRA 'MIDDLE' PHASE FOR SEARCH C USING THE LIKELY VAPOUR COMPOSITION C DO 200 I = 1 , NCOMP YTEST(I) = ZFRAC(I,JMAX) 200 CONTINUE C C SET TYPE AS LIQUID, AND INCREASE PRESSURE C PINC = 1.25D0 * P CALL FUGMOD(NCOMP,1,1,IC,T,PINC,Z,YTEST,FUG,FUGX) DO 205 I=1,NCOMP OLDY(NT-1,I) = FR(I) - FUG(I) -.25D0 205 CONTINUE LIST(NT-1) = 1 ENDIF LIST(NT) = 1 IF( STRTUP )THEN C----CALL PURE PHASE PROPERTIES, IF NOT CALLED BEFIORE CALL EOSPUR(NCOMP,SELECT,JPUR,OLDY) STRTUP = .FALSE. ENDIF DO 300 I = 1 , NCOMP C-----USE ONLY BASIS COMPONENTS SELECTED IN SELECT IF( .NOT.SELECT(I) )LIST(I) = 0 OLDY(NT,I) = FR(I) 300 CONTINUE IF( JPUR.GT.0 )THEN C----INSTABILITY IN PURE COMPONENT SEARCH DO 350 I = 1 , NCOMP Y(I) = 1.D-10 350 CONTINUE Y(JPUR) = 1. IF( EOSTYP )THEN ISTAB = 0 ELSE ISTAB = 1 ENDIF TPDMIN = -1. RETURN ENDIF C C LOOP THROUGH OTHER COMPONENTS NA TIMES C NTLOOP = NT IF( .NOT.EOSTYP )THEN DO 400 J = 1 , NFAS IF( IFLU(J).LT.0 )THEN NTLOOP = NT - 1 ISTAB = 1 ENDIF 400 CONTINUE ENDIF KVALUE = 0 DO 500 KK = 1 , NA KMIN = 3 C----SET BIG VALUE RTMIN = 1.D20 DO 450 K = 1 , NTLOOP C----LIST(K) EQUAL TO ZERO: COMP. NOT USED; EQUAL TO 3: SEARCH USELESS IF( LIST(K).NE.0 .AND. LIST(K).LT.3 )THEN C C PHASE TYPE SELECTION FOR THERMO.CALL C MFASE = 1 IF( K.EQ.NT )MFASE = -1 SUM = 0. DO 410 I = 1 , NCOMP Z = OLDY(K,I) C----LOG CONC. STORED IN OLDY IF( Z.LT.-60. )Z = -60. IF( Z.GT.10. )Z = 10. Z = DEXP(Z) YTEST(I) = Z SUM = SUM + Z 410 CONTINUE IF( SUM.GE..1 )THEN DLSUM = DLOG(SUM) CALL FUGMOD(NCOMP,1,MFASE,IC,T,P,Z,YTEST,FUG,FUGX) DMIN = NCOMP DO 415 J = 1 , NFAS DIST = 0. DO 412 I = 1 , NCOMP Z = ZFRAC(I,J) IF( J.EQ.1 )YTEST(I) = YTEST(I)/SUM DS(I) = YTEST(I) - Z IF( DABS(DS(I)).LT.1.D-30 )DS(I) = 0. DIST = DIST + DS(I)**2 412 CONTINUE C----FIND CLOSER PHASE IF( DIST.LT.DMIN )THEN DMIN = DIST MINVL = J DO 414 I = 1 , NCOMP DSG(I) = DS(I) 414 CONTINUE ENDIF 415 CONTINUE VAL = 0. GRAD = 0. C-----INVESTIGATE PROGRESS TOWARDS PHASE DO 420 I = 1 , NCOMP DNEW = FR(I) - FUG(I) GD = OLDY(K,I) - DNEW - DLSUM OLDY(K,I) = DNEW VAL = VAL + YTEST(I)*GD GRAD = GRAD + GD*DSG(I) 420 CONTINUE RT = 1 + VAL IF( GRAD.GT.0. )RT = 2*VAL/GRAD IF (GRAD .LT. 0.) RT = 0.5 + VAL IF( VAL.LT.0. )RT = VAL IF( DMIN.LT.2.D-3 .AND. RT.GT..8 )LIST(K) = 0 C---SCAN RESULT AND MODIFY LIST IN CASE OF TRIVIAL SOLUTION APPROACH IF( KK.NE.1 )THEN IF( RT.GE.0. )THEN IF( KK.GT.1 .AND. RT.GE.0.D0 )THEN IF( LIST(K).EQ.0 )GOTO 450 IF( RT.GE.RTOLD(K) )THEN LIST(K) = LIST(K) + 1 IF( RT.GT.1. )LIST(K) = LIST(K) + 1 ENDIF ENDIF ENDIF ENDIF RTOLD(K) = RT IF( LIST(K).LT.3 .AND. RT.LT.RTMIN )THEN C----IMPROVEMENT KVALUE = K KMIN = LIST(K) RTMIN = RT JMIN = MINVL DO 422 I = 1 , NCOMP Y(I) = YTEST(I) 422 CONTINUE ENDIF ENDIF ENDIF 450 CONTINUE C----CHECK IF WE HAD ANY IMROVEMENT IF( RTMIN.EQ.1.D20 )RETURN C----OR EXIT, IF WE FIND INSTABILITY IF( RTMIN.LT.0. )GOTO 600 500 CONTINUE C C FINSIHED HERE C 600 TPDMIN = 0. IF( RTMIN.LT.0. )TPDMIN = RTMIN IF( TPDMIN.LT.0.D0 )GOTO 1200 C C USE A TEST POINT, REFLECTED IN TRIVIAL SOLUTION PNT C COMPARE RESULT FROM THIS POINT TO THE RESULT FROM C THE SELECTED POINT AND CHOOSE BEST, IF NEITHER GIVE C CONCLUSIVE EVIDENCE C DO 700 I = 1 , NCOMP DS(I) = ZFRAC(I,JMIN) Z = DS(I)/Y(I) IF( Z.GT.2. )Z = 2. IF( Z.LT..5 )Z = .5 YTEST(I) = DS(I)*Z 700 CONTINUE DO 800 K = 1 , NA CALL FUGMOD(NCOMP,1,MFASE,IC,T,P,Z,YTEST,FUG,FUGX) IF( K.LT.NA )THEN DO 720 I = 1 , NCOMP Z = FR(I) - FUG(I) IF( Z.LT.-60. )Z = -60. YTEST(I) = DEXP(Z) 720 CONTINUE ENDIF 800 CONTINUE SUM = 0. DO 900 I = 1 , NCOMP SUM = SUM + YTEST(I) 900 CONTINUE VAL = 0. GRAD = 0. DO 1000 I = 1 , NCOMP YTEST(I) = YTEST(I)/SUM GD = DLOG(YTEST(I)) + FUG(I) - FR(I) VAL = VAL + YTEST(I)*GD GRAD = GRAD + GD*(YTEST(I)-DS(I)) 1000 CONTINUE IF( VAL.GT.0. )THEN RT = 1 + VAL IF( GRAD.GT.0. )RT = 2*VAL/GRAD IF( KMIN.NE.0 .AND. KMIN.NE.3 )THEN RTMIN = (1+KMIN)*RTMIN/2 IF( RT.GT.RTMIN )GOTO 1200 ENDIF ENDIF IF( VAL.LT.0. )TPDMIN = VAL DO 1100 I = 1 , NCOMP Y(I) = YTEST(I) 1100 CONTINUE 1200 IF( EOSTYP )THEN ISTAB = 0 ELSEIF( KVALUE.LT.NT )THEN ISTAB = 1 ELSE ISTAB = -1 ENDIF IF( NDPRT.GT.0 )WRITE(NDPRT,9010)(Y(I),I=1,NCOMP) RETURN 9010 FORMAT(/,' COMPOSITION FOR STABILITY TEST ',/,(3X,8F9.6)) END **==EOSPUR.FOR SUBROUTINE EOSPUR(NCOMP,SELECT,JPUR,OLDY) C C NCOMP: NO. OF COMPONENTS IN MIXTURE (INPUT) C JPUR: INDICATOR FOR PURE PHASE INSTABILITY (OUTPUT) C OLDY: LOGARITHM OF CONCENTRATIONS FOR INITI- C ATION OF STABILITY ANALYSIS (OUTPUT) C C CALLED FROM: PURE C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,NFMAX=5) LOGICAL SELECT(*) , EOSTYP DIMENSION OLDY(MDIM+2,*) C----LOCAL VARIABLES DIMENSION X(MDIM) COMMON /FPROP / FUG(MDIM) , FUGX(MDIM,MDIM) COMMON /KEEP / NFAS , T , P , GNUL , ZFEED(MDIM) , FR(MDIM) , & XVL(MDIM,NFMAX) , SFAS(NFMAX) , ZVAL(NFMAX) , & IFLU(MDIM) , EOSTYP , ISTFAS , ISTAB FMAX = 0. PFACT = 1.25 PLOG = DLOG(PFACT) DO 100 I = 1 , NCOMP IF( SELECT(I) )THEN DO 20 J = 1 , NCOMP X(J) = 0. 20 CONTINUE X(I) = 1. C C----GET PROPERTIES FOR PURE PHASE I C FULL ROUTINE USED, BUT A SPECIAL VERSION COULD BE C USED (NO NEED FOR MIXING RULES! ) C C SPECIFY LIQUID C CALL FUGMOD(NCOMP,1,1,IC,T,P,ZZZ,X,FUG,FUGX) C C CHECK FOR LIQUID FOUND; IF NOT, INCREASE PRESSURE C IF (IC .GT. 0) THEN DO 40 J = 1 , NCOMP OLDY(I,J) = FR(J) - FUG(J) 40 CONTINUE IF( OLDY(I,I).GT.FMAX )THEN C---THIS IMPLIES INSTABILITY FMAX = OLDY(I,I) JPUR = I ENDIF ELSE CALL FUGMOD(NCOMP,1,1,IC,T,PFACT*P,ZZZ,X,FUG,FUGX) DO 45 J = 1 , NCOMP OLDY(I,J) = FR(J) - FUG(J) - PLOG 45 CONTINUE ENDIF ENDIF 100 CONTINUE RETURN END **==FSTAB.FOR SUBROUTINE FSTAB(NCOMP,IDERIV,FUN,GRAD,XMAT,Y) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C MODIFIED OCTOBER 1989. C C C FSTAB IS THE OBJECTIVE FUNCTION ROUTINE FOR STABILITY ANALYSIS C USING THE SECOND ORDER MURRAY METHOD C C NCOMP: NO. OF COMPONENTS IN MIXTURE C IDERIV: DERIVATIVE SPECIFIER: 1: FIRST ONLY; 2: ALSO SECOND C FUN : VALUE OF OBJECTIVE FUNCTION C GRAD: GRADIENT VECTOR C XMAT: HESSIAN MATRIX C Y: VECTOR OF INDEPENDENT VARIABLES - HERE 2 TIMES SQRT(Y) C PARAMETER(MDIM=20,NFMAX=5,NEQM=(NFMAX-1)*MDIM) LOGICAL ACTIVE , EOSTYP DIMENSION GRAD(*) , XMAT(NCOMP,*) , Y(*) , YEX(MDIM) DIMENSION IRPOS(MDIM) COMMON /KEEP / NFAS , T , P , GNUL , ZFEED(MDIM) , FR(MDIM) , & XVL(MDIM,NFMAX) , SFAS(NFMAX) , ZVAL(NFMAX) , & IFLU(MDIM) , EOSTYP , ISTFAS , ISTAB COMMON /FPROP / FUG(MDIM) , FUGX(MDIM,MDIM) COMMON /AKTIV / ACTIVE(MDIM) , IVAR(NEQM+MDIM) , NACT , EPSV SAVE IRPOS DO 100 I = 1 , NCOMP C-----TRACE COMPONENTS BY DIRECT SUBSTITUTION IF( .NOT.ACTIVE(I) )THEN YI = FR(I) - FUG(I) Y(I) = 2.D0*DEXP(YI/2.D0) ENDIF 100 CONTINUE IF( IDERIV.EQ.2 )THEN NACT = 0 DO 150 I = 1 , NCOMP IRPOS(I) = 0 IF( Y(I).LT.1.D-4 )THEN ACTIVE(I) = .FALSE. ELSE NACT = NACT + 1 IVAR(NACT) = I ACTIVE(I) = .TRUE. ENDIF 150 CONTINUE DO 200 I = 1 , NACT IR = IVAR(I) IRPOS(IR) = I 200 CONTINUE ENDIF SUM = 0. DO 300 I = 1 , NCOMP ALFA = Y(I) IF( ALFA.GT.0.D0 )THEN YEX(I) = (ALFA/2)**2 ELSE YEX(I) = 2.5D-9 Y(I) = 1.D-4 ENDIF SUM = SUM + YEX(I) 300 CONTINUE JD = 1 IF( IDERIV.EQ.2 )JD = 2 CALL FUGMOD(NCOMP,JD,ISTAB,IC,T,P,Z,YEX,FUG,FUGX) IF( .NOT.EOSTYP )THEN ISTFAS = ISTAB ELSEIF( IC.GT.0 )THEN ISTFAS = 1 ELSE ISTFAS = -1 ENDIF FUN = 1. DO 400 I = 1 , NCOMP S = DLOG(YEX(I)) + FUG(I) - FR(I) IF( IDERIV.GT.0 )THEN IR = IRPOS(I) IF( IR.GT.0 )GRAD(IR) = S*Y(I)/2 ENDIF FUN = FUN + YEX(I)*(S-1.) 400 CONTINUE IF( IDERIV.EQ.2 )THEN DO 450 I = 1 , NCOMP IR = IRPOS(I) IF( IR.GT.0 )THEN S = Y(I)/SUM/4 DO 410 J = 1 , I JR = IRPOS(J) IF( JR.GT.0 )THEN XMAT(IR,JR) = S*Y(J)*FUGX(I,J) XMAT(JR,IR) = XMAT(IR,JR) ENDIF 410 CONTINUE C CORRECT XMAT(IR,IR) = XMAT(IR,IR) + 1.D0 + GRAD(I)/Y(I) XMAT(IR,IR) = XMAT(IR,IR) + 1.D0 + GRAD(IR)/Y(I) ENDIF 450 CONTINUE ENDIF RETURN END **==FUGMOD.FOR C C FUGMOD IS THE THERMODYNAMIC INTERFACE ROUTINE C USED FOR CALLING AN EQUATION OF STATE AS WELL C AS A HYBRID MODEL C C SUBROUTINE FUGMOD(NCOMP,INDIC,MTYP,IC,T,P,Z,X,FUG,FUGX) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C PRESENT EOS INTERFACES TO SRK/PR OR SPLIT MODEL ROUTINES C C INPUT: C NCOMP: NO. OF COMPONENTS IN MIXTURE (NOT USED) C INDIC: LEVEL OF DERIVATIVES: 1: FUG ONLY; 3: ALSO FUGX C MTYP: PHASE TYPE: 1: LIQ; -1: VAP; 0: MIN GIBBS ENERGY C IC : RETURN PHASE INDICATOR; 1: OK; -1: NOT OK; C T : TEMPERATURE C P : PRESSURE C OUTPUT: C Z : COMPRESSIBILITY FACTOR C X : MIXTURE COMPOSITION (MOLES) C FUG : LOG FUGACITY COEFFICIENTS C FUGX: SCALED MOLE DERIVATIVES OF FUG AT CONSTANT T,P C PARAMETER(MDIM=20) DIMENSION X(*) , FUG(*) , FUGX(MDIM,*) DIMENSION FUGT(MDIM) , FUGP(MDIM) , AUX(11) COMMON /STYR / N , NH , NDER , NTEMP COMMON /TYPEOS/ MODTYP IF( MODTYP.EQ.0 )THEN C C EOS C N = NCOMP NTEMP = 0 NDER = 1 IF( INDIC.GT.1 )NDER = 2 CALL CUBGEN(MTYP,IC,T,P,Z,X,FUG,FUGT,FUGP,FUGX,AUX) ELSE C C HYBRID MODEL; APPROPRIATE MODEL MUST BE SPECIFIED C CALL SPLITM(INDIC,MTYP,IC,T,P,Z,X,FUG,FUGT,FUGP,FUGX,AUX) ENDIF RETURN END **==EOSINP.FOR C C END OF GENERAL PROGRAM ROUTINES; C FOLLOWING ROUTINES: SRK/PR EQUATION OF STATE ROUTINES C C C C INPUT ROUTINE FOR EOS=MODEL C SUBROUTINE EOSINP(NC,SELECT,IEQ,LIST,INTYP,INFIL) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C INPUT: NC: NO. OF COMPONENTS IN MIXTURE C IEQ: EOS-SELECTION; 0=SRK, 1=PR C LIST: COMPONENT SELECTION LIST (FOR INTYP=0) C INTYP: INPUT INDICATOR: C INTYP = 0: BUILT-IN TABLE OF COMPONENTS C INTYP > 0: INPUT FROM FILE C INFIL: UNIT NO. FOR INPUT FILE C C OUTPUT: SELECT: ARRAY OF SELECTORS FOR TRIAL PHASE SELECTION C IN STABILITY ANALYSIS. C SELECT=0: COMPONENT NOT USED FOR TRIAL PHASE SELECTION C SELECT=1: COMPONENT USED FOR TRIAL PHASE SELECTION C C MAXBNK IS THE NUMBER OF ELEMENTS IN THE DATA BANK C MDIM IS THE MAX. NO OF ELEMENTS THAT CAN BE EXTRACTED C PARAMETER(MAXBNK=15) PARAMETER(MDIM=20,MC2=((MDIM)*(MDIM+1))/2) DIMENSION LIST(*) LOGICAL SELECT(*) DIMENSION ISOLID(MAXBNK) DIMENSION CSAPAR(MAXBNK) , CSBPAR(MAXBNK) , CSVPAR(MAXBNK) CHARACTER*4 NAME(MAXBNK) CHARACTER*8 NEW (MDIM) REAL T(MAXBNK) , P(MAXBNK) , OMEGA(MAXBNK) , CKV(MAXBNK,MAXBNK) REAL O(MDIM,MDIM) DIMENSION NHYCIF(MAXBNK) C C COMMON AREAS SET FOR THERMODYNAMIC PROPERTIES C COMMON /CUB / C , C1 , C2 , CRAT COMMON /CRIT / TCR(MDIM) , PCR(MDIM) , OMG(MDIM) COMMON /STYR / NCOMP , NH , NDER , NTEMP COMMON /OVER / AC(MDIM) , Q(MDIM) , TSQR(MDIM) COMMON /PAR / CK(MC2) , BC(MDIM) , AC0(MDIM) , AC1(MDIM) , & AC2(MDIM) COMMON /SOLPAR/ ISOLV(MDIM) , CSA(MDIM) , CSB(MDIM) , CSV(MDIM) C-----IDENTIFY HYDROCARBONS DATA NHYCIF/11*0 , 4*1/ DATA NAME/' C1 ' , ' C2 ' , ' C3 ' , ' IC4' , ' C4 ' , ' IC5' , & ' C5 ' , ' C6 ' , ' C7 ' , ' C8 ' , ' C9 ' , ' H20' , & ' N2 ' , ' C02' , ' H2S'/ C------CRITICAL TEMPERATURES (K) AND PRESSURES (ATM) (REID ET AL.) DATA T/190.6 , 305.4 , 369.8 , 408.1 , 425.2 , 460.4 , 469.6 , & 507.4 , 540.2 , 568.8 , 594.6 , 647.3 , 126.2 , 304.2 , & 373.2/ DATA P/45.4 , 48.2 , 41.9 , 36. , 37.5 , 33.4 , 33.3 , 29.3 , & 27. , 24.5 , 22.8 , 218. , 33.5 , 72.8 , 88.2/ C------ACENTRIC FACTORS (REID ET AL.) DATA OMEGA/.008 , .098 , .152 , .176 , .193 , .227 , .251 , .296 , & .351 , .394 , .440 , .344 , .04 , .225 , .1/ C------NONZERO KIJ FOR SRK-EQUATION (REID ET AL., HEIDEMANN ET AL.) DATA CKV/11*0. , .45 , .02 , .12 , .08 , 11*0. , .45 , .06 , .15 , & .07 , 11*0. , .53 , .08 , .15 , .07 , 11*0. , .52 , .08 , & .15 , .06 , 11*0. , .52 , .08 , .15 , .06 , 11*0. , .5 , & .08 , .15 , .06 , 11*0. , .5 , .08 , .15 , .06 , 11*0. , .5 , & .08 , .15 , .05 , 11*0. , .5 , .08 , .15 , .04 , 11*0. , .5 , & .08 , .15 , .04 , 11*0. , .5 , .08 , .15 , .03 , 2*.45 , & .53 , 2*.52 , 6*.5 , 4*0. , .02 , .06 , 9*.08 , 4*0. , .12 , & 10*.15 , 3*0. , .12 , .08 , 2*.07 , 4*.06 , .05 , 2*.04 , & .03 , 2*0. , .12 , 0./ C C SOLID PHASE PARAMETERS SET FOR NO SOLIDS C DATA ISOLID/15*0/ DATA CSAPAR/15*0./ DATA CSBPAR/15*0./ DATA CSVPAR/15*0./ IPRT = 1 C = IEQ CALL CONST(C,C1,C2,CONA,CONB) C C CRAT IS THE V/B RATIO AT THE PURE COMPONENT CRITICAL POINT C CRAT = (1./CONB-C)/3.D0 NCOMP = NC DO 100 I = 1 , NCOMP IF( INTYP.EQ.0 )THEN L = LIST(I) ISEL = NHYCIF(L) ISOLV(I) = ISOLID(I) CSA(I) = CSAPAR(L) CSB(I) = CSBPAR(L) CSV(I) = CSVPAR(L) NEW(I) = NAME(L) TCR(I) = T(L) PCR(I) = P(L) OM = OMEGA(L) ELSE READ(INFIL,*) NEW(I) , TCR(I) , PCR(I) , OM , ISEL , & ISOLV(I) ,CSA(I) , CSB(I) , CSV(I) ENDIF SELECT(I) = (ISEL.NE.0) TSQR(I) = 1./DSQRT(TCR(I)) BC(I) = CONB*TCR(I)/PCR(I) AC(I) = CONA*TCR(I)/DSQRT(PCR(I)) OMG(I) = OM C-----ICEQ=0: SRK Q(I) = .48 + OM*(1.574-.176*OM) C-----IEQ=1: PENG-ROBINSON IF( IEQ.EQ.1 )Q(I) = .37464 + OM*(1.54226-.26992*OM) 100 CONTINUE IF( INTYP.EQ.0 )THEN CRMAX = 0. CRMIN = 1.D5 LMAX = 0 LMIN = 0 C------CHOOSE LIGHT AND HEAVY HYDROCARB. FROM CRIT. TEMP DO 150 I = 1 , NCOMP IF( .NOT.SELECT(I) )THEN TTRY = TCR(I) IF( TTRY.GE.CRMAX )THEN CRMAX = TTRY LMAX = I ENDIF IF( TTRY.LE.CRMIN )THEN CRMIN = TTRY LMIN = I ENDIF ENDIF 150 CONTINUE IF( LMIN.NE.0 )SELECT(LMIN) = .TRUE. IF( LMAX.NE.0 )SELECT(LMAX) = .TRUE. ENDIF IF( IPRT.GE.0 )THEN WRITE(22,9010) DO 200 I = 1 , NCOMP WRITE(22,9020)NEW(I) , TCR(I) , PCR(I) ,OMG(I), SELECT(I) 200 CONTINUE ENDIF IF( INTYP.EQ.0 )THEN DO 250 I = 1 , NCOMP L1 = LIST(I) DO 220 J = I , NCOMP L2 = LIST(J) O(I,J) = CKV(L1,L2) O(J,I) = O(I,J) 220 CONTINUE 250 CONTINUE ELSE O(1,1) = 0. DO 300 I = 2 , NCOMP O(I,I) = 0. READ(INFIL,*) (O(I,J),J=1,I-1) DO 260 J = 1 , I-1 O(J,I) = O(I,J) 260 CONTINUE 300 CONTINUE ENDIF NAD = 0 DO 400 I = 1 , NCOMP DO 350 J = I , NCOMP NAD = NAD + 1 CK(NAD) = O(I,J) 350 CONTINUE 400 CONTINUE IF( IPRT.LT.0 )RETURN WRITE(22,9030)(NEW(K)(1:4),K=1,NCOMP-1) DO 500 J = 2 , NCOMP WRITE(22,9040)NEW(J)(1:4),(O(J,I),I=1,J-1) 500 CONTINUE RETURN 9010 FORMAT(' COMP TCRIT PCRIT OMEGA SELECT * ',/) 9020 FORMAT(1X,A8,2F9.2,1X,F8.4,4X,L2) 9030 FORMAT(/,' BINARY INTERACTION COEFFICIENTS ',//,(T9,12(A4,3X))) 9040 FORMAT(1X,A4,(T8,12F7.3)) END **==TEMSET.FOR SUBROUTINE TEMSET(T) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,MC2=(MDIM*(MDIM+1))/2) COMMON /STYR / NCOMP , NH , NDER , NTEMP COMMON /OVER / AC(MDIM) , Q(MDIM) , TSQR(MDIM) COMMON /PAR / CK(MC2) , BC(MDIM) , AC0(MDIM) , AC1(MDIM) , & AC2(MDIM) SAVE TOLD DATA TOLD/0.D0/ IF( T.EQ.TOLD )RETURN C C COEFFICIENTS FOR TEMPERATURE DEPENDENT COEFFICIENTS ONLY HAVE C TO BE RECALCULATED WHEN T IS CHANGED. C TOLD = T SQT = DSQRT(T) T2R = 1.D0/T/2 T2RF = -3*T2R DO 100 I = 1 , NCOMP ALF = TSQR(I) Q1 = AC(I)*(1+Q(I))/SQT AC0(I) = Q1 - AC(I)*Q(I)*ALF AC1(I) = -Q1*T2R AC2(I) = T2RF*AC1(I) 100 CONTINUE RETURN END **==CUBGEN.FOR SUBROUTINE CUBGEN(MT,IC,T,P,Z,XX,FG,FT,FP,FX,AUX) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,MC2=(MDIM*(MDIM+1))/2) C-----DUMMY VARIABLES DIMENSION XX(*) , FG(*) , FT(*) , FP(*) , FX(MDIM,*) , AUX(*) C-----LOCAL VARIABLES DIMENSION X(MDIM) , PD(MDIM) , AD1(MDIM) , BD1(MDIM) , ADT(MDIM) , & AD2(MC2) COMMON /STYR / NCOMP , NH , NDER , NTEMP COMMON /CUB / C , C1 , C2 , CRAT CALL TEMSET(T) SUM = 0. DO 100 I = 1 , NCOMP SUM = SUM + XX(I) 100 CONTINUE DO 200 I = 1 , NCOMP X(I) = XX(I)/SUM 200 CONTINUE CALL ANEW(X,A,AT,ATT,B,AD1,AD2,ADT,BD1) APT = A*P/T BPT = B*P/T CALL CUBIC(MT,APT,BPT,Z,Z2,GDER) C-------SECOND Z AND ENERGY DIFF. STORED IN AUX(1),AUX(2) AUX(1) = Z2 AUX(2) = GDER V = Z*T/P VRAT = V/B C C V/B IS COMPARED WITH (V/B) AT THE PSEUDOCRITICAL POINT C TO IDENTIFY THE PHASE C IC = 1 TL = CRAT/VRAT - 1.0 IF( MT.EQ.0 .AND. TL.GT.0. )IC = 2 IF( MT.EQ.0 .AND. TL.LT.0. )IC = -2 IF( MT*TL.LT.0. )IC = -1 C------V = VOLUME/R S1 = 1./(V+C1*B) S2 = 1./(V+C2*B) P1 = P/T + A*S1*S2 PA = -S1*S2 PN = P1 FAC = C1*S1 + C2*S2 P2 = A*PA PB = P1*P1 - FAC*P2 C------ DERIVATIVE IS D(P/T) / D(V/R) DPDV = -P1*P1 - P2*(S1+S2) AUX(8) = DPDV*T XL1 = DLOG(V*P1) FN = XL1 XL2 = DLOG(S1/S2)/(C2-C1) FA = -XL2/B F2 = -A*FA FF = FN - F2 GB = -V*PA F2B = (A*GB-F2)/B FB = P1 - F2B FNB = P1 FAB = -F2B/A GBB = -GB*FAC F2BB = (A*GBB-2*F2B)/B FBB = P1*P1 - F2BB XLZ = DLOG(Z) FNP = FN - XLZ IF( NTEMP.NE.0 )THEN DFT = FA*AT HE = -T*DFT + Z - 1 HE = HE*T C------EXCESS ENTHALPY/R STORED IN AUX(3) AUX(3) = HE SE = -T*DFT - FF + XLZ C------EXCESS ENTROPY/R STORED IN AUX(4) AUX(4) = SE PTR = PA*AT DPDT = P/T + T*PTR IF( NTEMP.GE.2 )THEN FTT = FA*ATT CP = -T*(T*FTT+2*DFT) - DPDT**2/DPDV - 1 C------EXCESS HEAT CAPACITY/R IN AUX(5) AUX(5) = CP DVDT = -DPDT/DPDV/T C---- AUX(6) IS PRESSURE DERIVATIVE OF RESID. ENTROPY AUX(6) = -DVDT + 1/P C-----AUX(7) IS PRESSURE DERIVATIVE OF ENTHALPY AUX(7) = -T*DVDT + V AUX(9) = DVDT ENDIF ENDIF IF( NDER.NE.0 )THEN DO 250 I = 1 , NCOMP FG(I) = FNP + FA*AD1(I) + FB*BD1(I) 250 CONTINUE IF( NDER.GE.2 .OR. NTEMP.NE.0 )THEN DO 260 I = 1 , NCOMP PD(I) = (PN+PA*AD1(I)+PB*BD1(I))/DPDV 260 CONTINUE IF( NDER.GE.2 )THEN NADR = 0 DO 270 I = 1 , NCOMP PDI = PD(I) CC = 1 + FNB*BD1(I) + PN*PDI CA = FAB*BD1(I) + PA*PDI CB = FNB + FAB*AD1(I) + FBB*BD1(I) + PB*PDI DO 265 J = I , NCOMP NADR = NADR + 1 FX(I,J) = CC + CA*AD1(J) + CB*BD1(J) + FA*AD2(NADR) FX(J,I) = FX(I,J) 265 CONTINUE 270 CONTINUE ENDIF IF( NTEMP.NE.0 )THEN CC = 1./T CB = FAB*AT DPDTT = DPDT/T DO 280 I = 1 , NCOMP FP(I) = -PD(I)/T - 1./P FT(I) = CC + CB*BD1(I) + FA*ADT(I) + DPDTT*PD(I) 280 CONTINUE ENDIF ENDIF ENDIF RETURN END **==ANEW.FOR SUBROUTINE ANEW(X,A0,AT,ATT,B0,AD1,AD2,ADT,BD1) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20,MC2=(MDIM*(MDIM+1))/2) DIMENSION X(*) , AD1(*) , AD2(*) , ADT(*) , BD1(*) DIMENSION AX0(MDIM) , AX1(MDIM) , AS1(MDIM) , AST1(MDIM) , & AX2(MDIM) COMMON /STYR / NCOMP , NH , NDER , NTEMP COMMON /PAR / CK(MC2) , BC(MDIM) , AC0(MDIM) , AC1(MDIM) , & AC2(MDIM) S0 = 0. S1 = 0. S2 = 0. DO 100 I = 1 , NCOMP XI = X(I) AS1(I) = 0. AX0(I) = XI*AC0(I) IF( NTEMP.NE.0 )THEN C------NTEMP=0: NO TEMP. DERIVATIVES AST1(I) = 0. AX1(I) = XI*AC1(I) S1 = S1 + AX1(I) IF( NTEMP.NE.1 )THEN C------NTEMP=1: NO SEC. TEMP DERIVATIVES AX2(I) = XI*AC2(I) S2 = S2 + AX2(I) ENDIF ENDIF S0 = S0 + AX0(I) 100 CONTINUE A0 = S0*S0 NADR = 0 GSM = S0*S2 + S1*S1 DO 200 I = 1 , NCOMP FAC = 2*AC0(I) AD1(I) = S0*FAC IF( NTEMP.NE.0 )ADT(I) = FAC*S1 + 2*S0*AC1(I) IF( NDER.GE.2 )THEN DO 120 J = I , NCOMP NADR = NADR + 1 AD2(NADR) = FAC*AC0(J) 120 CONTINUE ENDIF 200 CONTINUE NADR = 0 DO 300 I = 1 , NCOMP C-----DOUBLE LOOP FOR A-PARAMETER IF( I.NE.NCOMP )THEN K1 = I + 1 NADR = NADR + 1 XI = X(I) AX0I = AX0(I) IF( NTEMP.NE.0. )AX1I = AX1(I) FAC = 2*AC0(I) DO 220 K = K1 , NCOMP NADR = NADR + 1 EF = CK(NADR) IF( EF.NE.0. )THEN AX0K = AX0(K) AS1(I) = AS1(I) + AX0K*EF AS1(K) = AS1(K) + AX0I*EF IF( NTEMP.NE.0 )THEN AST1(I) = AST1(I) + AX1(K)*EF AST1(K) = AST1(K) + AX1I*EF IF( NTEMP.NE.1 )GSM = GSM - & EF*(AX0K*AX2(I)+AX0I*AX2(K)+2*AX1I*AX1(K)) ENDIF IF( NDER.GE.2 )THEN C------NDER=1: NO SECOND COMPOSITION DERIVATIVES FAC1 = AC0(K)*FAC*EF AD2(NADR) = AD2(NADR) - FAC1 ENDIF ENDIF 220 CONTINUE ENDIF 300 CONTINUE A0 = 0 AT = 0. IF( NTEMP.GT.1 )ATT = 2*GSM DO 400 I = 1 , NCOMP FAC = 2*AC0(I) AD1(I) = AD1(I) - FAC*AS1(I) A0 = A0 + AD1(I)*X(I)/2 IF( NTEMP.NE.0 )THEN ADT(I) = ADT(I) - 2*AC1(I)*AS1(I) - FAC*AST1(I) AT = AT + X(I)*ADT(I) ENDIF 400 CONTINUE AT = AT/2 B0 = 0. C------ HERE USED: SINGLE SUM IN B DO 500 I = 1 , NCOMP B0 = B0 + X(I)*BC(I) BD1(I) = BC(I) 500 CONTINUE RETURN END **==CUBIC.FOR SUBROUTINE CUBIC(MTYP,A,B,Z,ZV2,DELG) C------Z IS DESIRED ROOT; ZV2 IS THE OTHER ROOT; DELG IS G-DIFFERENCE IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /CUB / C , C1 , C2 , CRAT ZV2 = 0. DELG = 0. BC = B*C X2 = 1 - BC X1 = A - B*(1+B+C+2*BC) X0 = B*(A-BC*(1+B)) Z = X2/3 F = ((Z-X2)*Z+X1)*Z - X0 Z = B IF( F.LT.0. .OR. B.GT.X2/3 )Z = Z + 1 100 DF2 = 3*Z - X2 DF1 = Z*(DF2-X2) + X1 F = ((Z-X2)*Z+X1)*Z - X0 DZ = F/DF1 DZ = DZ*(1+DZ*DF2/DF1) Z = Z - DZ IF( DABS(DZ).GT.1.D-10 )GOTO 100 IF( MTYP*DF2.GE.0. )THEN E1 = Z - X2 E0 = Z*E1 + X1 D = E1*E1 - 4*E0 IF( D.GE.0. )THEN Z1 = (DABS(E1)+DSQRT(D))/2 IF( E1.GT.0. )Z1 = -Z1 IF( Z.GT.Z1 )Z1 = E0/Z1 DF2 = 3*Z1 - X2 DF1 = Z1*(DF2-X2) + X1 F = ((Z1-X2)*Z1+X1)*Z1 - X0 DZ = F/DF1 DZ = DZ*(1+DZ*DF2/DF1) Z1 = Z1 - DZ IF( Z1.GE.B )THEN IF( MTYP.EQ.0 )THEN C-----CALCULATE EXCESS GIBBS ENERGY DIFFERENCE FOR MULTIPLE SOLUTIONS F = DLOG((Z-B)/(Z1-B)) + A/B/(C2-C1) & *DLOG((Z+C2*B)*(Z1+C1*B)/(Z+C1*B)/(Z1+C2*B)) & - (Z-Z1) DELG = DABS(F) ZV2 = Z1 IF( F.LT.0. )THEN ZV2 = Z Z = Z1 ENDIF ELSEIF( MTYP.EQ.1 )THEN IF( Z1.LT.Z )Z = Z1 ELSE IF( Z1.GT.Z )Z = Z1 ENDIF ENDIF ENDIF ENDIF RETURN END **==CONST.FOR SUBROUTINE CONST(C,C1,C2,AC,BC) IMPLICIT DOUBLE PRECISION(A-H,O-Z) U = 1 + C W = -C D = U*U - 4*W C2 = (U+DSQRT(D))/2 C1 = W/C2 S1 = 1 + ((U+W)*(U+3)-W)/2 S2 = 1 + U + W R = DSQRT(S1*S1-S2*S2*S2) A = (S1+R)**(1.D0/3) B = S2/A X = A + B + 1 BC = 1./(3*X+U-1) AC = BC*DSQRT(X*(X*X-3*W)-U*W) RETURN END **==SOLFUG.FOR C C SOLFUG IS AN EXAMPLE ROUTINE FOR CALCULATION OF SOLID FUGACITIES, C WHICH CAN REPLACED BY ANY APPROPRIATE ALTERNATIVE C C THE LN (FUGACITY COEFFICIENT) IS ASSUMED TO BE GIVEN BY AN C EXPRESSION OF THE FORM C C LN (FG) = A + B/T + V*P/T - LN(P) C C BUT ANY ALTERNATIVE FORM COULD BE USED AS WELL C C SUBROUTINE SOLFUG(N,T,P,ISOLC,FUGAR) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20) DIMENSION ISOLC(*) , FUGAR(*) COMMON /SOLPAR/ ISOLV(MDIM) , CSA(MDIM) , CSB(MDIM) , CSV(MDIM) C C THIS ROUTINE MUST CALCULATE LN FUGACITY COEFFICIENTS FOR SOLID C PHASE FORMERS AND RETURN VALUE IN FUGAR C C FOR NON-SOLID FORMERS, RETURN A VALUE OF 100 C DO 100 I = 1 , N IF( ISOLV(I).EQ.0 )THEN ISOLC(I) = 0 FUGAR(I) = 100.D0 ELSE ISOLC(I) = 1 FUGLOG = CSA(I) + CSB(I)/T + CSV(I)*P/T FUGAR(I) = FUGLOG - DLOG(P) ENDIF 100 CONTINUE RETURN END **==SPLITM.FOR C C DUMMY SUBROUTINE FOR HYBRID MODEL C SUBROUTINE SPLITM(INDIC,MTYP,IC,T,P,Z,X,FUG,FUGT,FUGP,FUGX,AUX) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MDIM=20) DIMENSION X(*) , FUG(*) , FUGT(*) , FUGP(*) , FUGX(MDIM,*) , & AUX(11) C C C ENTER BODY OF ROUTINE HERE C C RETURN END **==HYBINP.FOR C C DUMMY SUBROUTINE FOR INITIALIZATION OF HYBRID MODEL C SUBROUTINE HYBINP(NCOMP) RETURN END
Write, Run & Share Fortran code online using OneCompiler's Fortran online compiler for free. It's one of the robust, feature-rich online compilers for Fortran language, running on the latest version 7. Getting started with the OneCompiler's Fortran compiler is simple and pretty fast. The editor shows sample boilerplate code when you choose language as Fortran and start coding.
OneCompiler's Fortran online editor supports stdin and users can give inputs to programs using the STDIN textbox under the I/O tab. Following is a sample Fortran program which takes name as input and prints hello message with your name.
program hello
character :: name*30
read *, name
print *, "Hello ", name
end program hello
Fortran language was initially developed for scientific calculations by IBM in 1957. It has a number of in-built functions to perform mathematical calculations and is ideal for applications which has more mathematical calculations.
| Data type | Description | Usage |
|---|---|---|
| Integer | To store integer variables | integer :: x |
| Real | To store float values | real :: x |
| Complex | To store complex numbers | complex :: x,y |
| Logical | To store boolean values True or false | logical :: x=.True. , logical :: x = .FALSE. |
| Character | To store characters and strings | character :: x |
Variable is a name given to the storage area in order to manipulate them in our programs.
data type :: variable_name
Array is a collection of similar data which is stored in continuous memory addresses.
data-type, dimension (x,y) :: array-name
integer, dimension(3,3) :: cube
Do is used to execute a set of statement(s) iteratively when a given condition is true and the loop variable must be an integer.
do i = start, stop [,step]
! code
end do
Do-While is used to execute a set of statement(s) iteratively when a given condition is true.
do while (condition)
!Code
end do
If is used to execute a set of statements based on a condition.
if (logical-expression) then
!Code
end if
If is used to execute a set of statements based on a condition and execute another set of statements present in else block, if condition specified in If block fails.
if (logical-expression) then
!code when the condition is true
else
!code when the condition fails
end if
Case is similar to switch in C language.
[name:] select case (regular-expression)
case (value1)
! code for value 1
... case (value2)
! code for value 2
...
case default
! default code
...
end select [name]