c ---------------------------------------------------------- c Indicating the Region of Mixture c ---------------------------------------------------------- INTEGER FUNCTION IPHASE(TT,PP,Z) INTEGER KPA,MESS,KSTAN,KAS,J $ ,FMF043 REAL T,P,Z,TT,PP CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56),TEMPR,PRESSR,DBMOL,DBZ c$$$ $ ,PBUBR,PDEWR c$$$ $ ,FMF039,FMF040 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C CALL FMF000(J,CONST,COEFF) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF C PRNAME='IPHASE' c ------ Transforming the unit of properties ----------------- T=TT P=PP IF ((KPA.EQ.0).OR.(KPA.GE.4)) THEN P=PP*1.0E-5 ELSEIF (KPA.EQ.1) THEN T=TT+273.15 ELSEIF (KPA.EQ.3) THEN P=PP*1.0E-5 T=TT+273.15 ENDIF C c ------ Checking the range of Temperature and Pressure ------ c ------ and fraction ------ IF ((T.LT.APPLIM(1)).AND.(T.GT.APPLIM(2))) THEN CALL FMF049(2,PRNAME) IPHASE=-2 RETURN ELSEIF ((P.LT.APPLIM(3)).OR.(P.GT.APPLIM(4))) THEN CALL FMF049(2,PRNAME) IPHASE=-2 RETURN ELSEIF (Z.LT.0.0D0.OR.Z.GT.1.0D0) THEN CALL FMF049(2,PRNAME) IPHASE=-2 RETURN ENDIF C c ============ Double precision ============================== c c --------- Reading properties of NH3 and H2O --------------- C DBZ=DBLE(Z) IF (KAS.EQ.1) THEN DBMOL=DBZ*CONST(21)/((1.0D0-DBZ)*CONST(11)+DBZ*CONST(21)) ELSE DBMOL=DBZ ENDIF C c ------- real properties => reduced properties ------------- TEMPR=DBLE(T)/CONST(1) PRESSR=DBLE(P)/CONST(2) C c$$$ PBUBR=FMF039(CONST,COEFF,TEMPR,DBMOL) c$$$ IF (PBUBR.LE.-1.0D15) THEN c$$$ CALL FMF049(2,PRNAME) c$$$ IPHASE=-2 c$$$ RETURN c$$$ ELSEIF (PBUBR.LT.0.0D0) THEN c$$$ CALL FMF049(1,PRNAME) c$$$ IPHASE=-1 c$$$ RETURN c$$$ ENDIF c$$$ PDEWR=FMF040(CONST,COEFF,TEMPR,DBMOL) c$$$ IF (PDEWR.LE.-1.0D15) THEN c$$$ CALL FMF049(2,PRNAME) c$$$ IPHASE=-2 c$$$ RETURN c$$$ ELSEIF (PDEWR.LT.0.0D0) THEN c$$$ CALL FMF049(1,PRNAME) c$$$ IPHASE=-1 c$$$ RETURN c$$$ ENDIF c$$$ IF (PRESSR.GE.PBUBR) THEN c$$$ IPHASE=1 c$$$ RETURN c$$$ ELSEIF ((PRESSR.LT.PBUBR).AND.(PRESSR.GT.PDEWR)) THEN c$$$ IPHASE=2 c$$$ RETURN c$$$ ELSEIF (PRESSR.LE.PDEWR) THEN c$$$ IPHASE=3 c$$$ RETURN c$$$ ENDIF IPHASE=FMF043(CONST,COEFF,TEMPR,PRESSR,DBMOL) RETURN END c --------------------------------------------------------- c The Identification of Ammonia and Water c --------------------------------------------------------- CHARACTER*40 FUNCTION IDENTM(I,A) INTEGER I,KPA,MESS,KSTAN,KAS COMMON /UNIT/ KPA,MESS,KSTAN,KAS CHARACTER*1 A c IF (I.EQ.1) THEN IF (A.EQ.'C') THEN IDENTM='NH3' ELSEIF (A.EQ.'S') THEN IDENTM='AMMONIA' ELSEIF (A.EQ.'V') THEN IDENTM='12.1' ELSE IDENTM='Out of range' ENDIF ELSEIF (I.EQ.2) THEN IF (A.EQ.'C') THEN IDENTM='H2O' ELSEIF (A.EQ.'S') THEN IDENTM='WATER' ELSEIF (A.EQ.'V') THEN IDENTM='12.1' ELSE IDENTM='Out of range' ENDIF ELSE IDENTM='Out of range' ENDIF RETURN END c -------------------------------------------------------------- c Subroutine subpur c calculate the properties of pure substance c Temperature,Pressure => Volume,Enthalpy,Entropy c i=1 : Ammonia c i=2 : Water c -------------------------------------------------------------- SUBROUTINE SUBPUR(I,J,TT,PP,V,H,S) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL T,P,V,H,S,TT,PP CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56),TEMPR,PRESSR $ ,PSATR,DBVR,DBHR,DBSR,DBV,DBH,DBS,TEMPJ $ ,FMF023,FMF024,FMF025,FMF026,FMF027,FMF028,FMF037,FMF038 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C J=0 PRNAME='SUBPUR' c c --------- Reading properties of NH3 and H2O --------------- CALL FMF000(J,CONST,COEFF) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF C c ------ Transforming the unit of properties ----------------- T=TT P=PP IF ((KPA.EQ.0).OR.(KPA.GE.4)) THEN P=PP*1.0E-5 ELSEIF (KPA.EQ.1) THEN T=TT+273.15 ELSEIF (KPA.EQ.3) THEN P=PP*1.0E-5 T=TT+273.15 ENDIF C c ------ Checking the range of Temperature and Pressure ------ c ------ and Fraction ------ IF ((T.LT.APPLIM(1)).AND.(T.GT.APPLIM(2))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF ((P.LT.APPLIM(3)).OR.(P.GT.APPLIM(4))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c ============ Double precision ============================== C c ------- real properties => reduced properties ------------- TEMPR=DBLE(T)/CONST(1) PRESSR=DBLE(P)/CONST(2) C c ----------- NH3 (ammonia) -------------------- IF (I.EQ.1) THEN TEMPJ=DBLE(T) IF (TEMPJ.GT.CONST(12)) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF c -------- judging phase ------------------- PSATR=FMF038(CONST,COEFF,TEMPR) IF (PSATR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PSATR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF C IF (PRESSR.GT.PSATR) THEN IF (PRESSR.GE.PSATR) THEN DBVR=FMF024(COEFF,TEMPR,PRESSR,1.0D0) DBHR=FMF026(COEFF,TEMPR,PRESSR,1.0D0) DBSR=FMF028(COEFF,TEMPR,PRESSR,1.0D0) ELSEIF (PRESSR.LT.PSATR) THEN DBVR=FMF023(COEFF,TEMPR,PRESSR,1.0D0) DBHR=FMF025(COEFF,TEMPR,PRESSR,1.0D0) DBSR=FMF027(COEFF,TEMPR,PRESSR,1.0D0) c ELSE c WRITE(*,*) ' SATURATED STATE' c J=-2 c RETURN ENDIF c ------------ H2O (water) ------------------- ELSEIF (I.EQ.2) THEN c ----- judging phase ----------- PSATR=FMF037(CONST,COEFF,TEMPR) IF (PSATR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PSATR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF C IF (PRESSR.GT.PSATR) THEN IF (PRESSR.GE.PSATR) THEN DBVR=FMF024(COEFF,TEMPR,PRESSR,0.0D0) DBHR=FMF026(COEFF,TEMPR,PRESSR,0.0D0) DBSR=FMF028(COEFF,TEMPR,PRESSR,0.0D0) ELSEIF (PRESSR.LT.PSATR) THEN DBVR=FMF023(COEFF,TEMPR,PRESSR,0.0D0) DBHR=FMF025(COEFF,TEMPR,PRESSR,0.0D0) DBSR=FMF027(COEFF,TEMPR,PRESSR,0.0D0) c ELSE c WRITE(*,*) ' SATURATED STATE' c J=-2 c RETURN ENDIF ELSE J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c ----------- reduced properties => real properties ---------- c ----------- kJ => J ----------- DBV=DBVR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBH=(DBHR*CONST(1)*CONST(3)+CONST(3+I))*1.0D3 DBS=(DBSR*CONST(3)+CONST(5+I))*1.0D3 C c ----------- Transforming the unit of properties ------------ IF (KSTAN.EQ.1) THEN DBH=DBH+CONST(10*I+9) DBS=DBS+CONST(10*I+10) ENDIF IF (KAS.EQ.1) THEN DBV=DBV/CONST(10*I+1) DBH=DBH/CONST(10*I+1) DBS=DBS/CONST(10*I+1) ENDIF C c ============= real ================= c V=SNGL(DBV) H=SNGL(DBH) S=SNGL(DBS) RETURN END c ------------------------------------------------------------------- c Subroutine subpst c calculate the properties at saturated state of pure substance c Temperature => Pressure,Volume,Enthalpy,Entropy c i=1 : Ammonia c i=2 : Water c -------------------------------------------------------------------- SUBROUTINE SUBPST(I,J,TT,PS,VL,VV,HL,HV,SL,SV) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL T,PS,VL,VV,HL,HV,SL,SV,TT CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56) $ ,TEMPR,PSATR,DBVLR,DBVVR,DBHLR,DBHVR,DBSLR,DBSVR $ ,DBVL,DBVV,DBHL,DBHV,DBSL,DBSV,PSAT,NH3LIM,H2OLIM $ ,FMF023,FMF024,FMF025,FMF026,FMF027,FMF028,FMF037,FMF038 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C J=0 PRNAME='SUBPST' c c ------------ Reading properties of H2O and NH3 -------- CALL FMF000(J,CONST,COEFF) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF C c ------- Transforming the unit of properties ---------- T=TT IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ENDIF C c ------------ Checking the range of temperature -------- IF ((T.LT.APPLIM(1)).OR.(T.GT.APPLIM(2))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c =========== Double precision ================ C c ----------- real property => reduced property --------- TEMPR=DBLE(T)/CONST(1) c c ----------- NH3(ammonia) ---------------- IF (I.EQ.1) THEN NH3LIM=DBLE(APPLIM(5))/CONST(1) IF (T.GT.NH3LIM) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF PSATR=FMF038(CONST,COEFF,TEMPR) IF (PSATR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PSATR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF DBVLR=FMF023(COEFF,TEMPR,PSATR,1.0D0) DBVVR=FMF024(COEFF,TEMPR,PSATR,1.0D0) DBHLR=FMF025(COEFF,TEMPR,PSATR,1.0D0) DBHVR=FMF026(COEFF,TEMPR,PSATR,1.0D0) DBSLR=FMF027(COEFF,TEMPR,PSATR,1.0D0) DBSVR=FMF028(COEFF,TEMPR,PSATR,1.0D0) C c ------------- H2O(water) ----------------------- ELSEIF (I.EQ.2) THEN H2OLIM=DBLE(APPLIM(7))/CONST(1) IF (T.GT.H2OLIM) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF PSATR=FMF037(CONST,COEFF,TEMPR) IF (PSATR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PSATR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF DBVLR=FMF023(COEFF,TEMPR,PSATR,0.0D0) DBVVR=FMF024(COEFF,TEMPR,PSATR,0.0D0) DBHLR=FMF025(COEFF,TEMPR,PSATR,0.0D0) DBHVR=FMF026(COEFF,TEMPR,PSATR,0.0D0) DBSLR=FMF027(COEFF,TEMPR,PSATR,0.0D0) DBSVR=FMF028(COEFF,TEMPR,PSATR,0.0D0) ELSE J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c -------- reduced properties => real properties ------- PSAT=PSATR*CONST(2) DBVL=DBVLR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBVV=DBVVR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBHL=(DBHLR*CONST(1)*CONST(3)+CONST(3+I))*1.0D3 DBHV=(DBHVR*CONST(1)*CONST(3)+CONST(3+I))*1.0D3 DBSL=(DBSLR*CONST(3)+CONST(5+I))*1.0D3 DBSV=(DBSVR*CONST(3)+CONST(5+I))*1.0D3 C c -------- Transforming the unit of properties ------------- IF ((KPA.EQ.0).OR.(KPA.GE.3)) THEN PS=SNGL(PSAT*1.0D5) ELSE PS=SNGL(PSAT) ENDIF IF (KSTAN.EQ.1) THEN DBHL=DBHL+CONST(10*I+9) DBHV=DBHV+CONST(10*I+9) DBSL=DBSL+CONST(10*I+10) DBSV=DBSV+CONST(10*I+10) ENDIF IF (KAS.EQ.1) THEN VL=SNGL(DBVL/CONST(10*I+1)) VV=SNGL(DBVV/CONST(10*I+1)) HL=SNGL(DBHL/CONST(10*I+1)) HV=SNGL(DBHV/CONST(10*I+1)) SL=SNGL(DBSL/CONST(10*I+1)) SV=SNGL(DBSV/CONST(10*I+1)) ELSE VL=SNGL(DBVL) VV=SNGL(DBVV) HL=SNGL(DBHL) HV=SNGL(DBHV) SL=SNGL(DBSL) SV=SNGL(DBSV) ENDIF RETURN END c ------------------------------------------------------------------- c Temperature,overall composition => Pressure at bubble point c ------------------------------------------------------------------- SUBROUTINE SUBPB(J,TT,P,Z,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,Z,V,H,S,TT CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56),PRESSR,PRESS $ ,TEMPR,DBVLR,DBHLR,DBSLR,DBVL,DBHL,DBSL,DBX $ ,FMF023,FMF025,FMF027,FMF039 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C J=0 PRNAME='SUBPB ' c c ------------ Reading properties of H2O and NH3 -------- CALL FMF000(J,CONST,COEFF) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF C C c ------- Transforming the unit of properties ---------- T=TT IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ENDIF c c ------------ Checking the range of temperature -------- IF ((T.LT.APPLIM(1)).OR.(T.GT.APPLIM(2))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (Z.LT.0.0E0.OR.Z.GT.1.0E0) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c =========== Double precision ================ c ----------- real property => reduced property --------- TEMPR=DBLE(T)/CONST(1) DBX=DBLE(Z) IF (KAS.EQ.1) THEN DBX=DBX*CONST(21)/((1.0D0-DBX)*CONST(11)+DBX*CONST(21)) ENDIF C C PRESSR=FMF039(CONST,COEFF,TEMPR,DBX) IF (PRESSR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PRESSR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF DBVLR=FMF023(COEFF,TEMPR,PRESSR,DBX) DBHLR=FMF025(COEFF,TEMPR,PRESSR,DBX) DBSLR=FMF027(COEFF,TEMPR,PRESSR,DBX) C c --------- reduced properties => real properties ----------- c --------- kJ => J ----------------------------------------- PRESS=PRESSR*CONST(2) DBVL=DBVLR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBHL=(DBHLR*CONST(1)*CONST(3)+DBX*CONST(4)+(1.0D0-DBX)*CONST(5)) $ *1.0D3 DBSL=(DBSLR*CONST(3)+DBX*CONST(6)+(1.0D0-DBX)*CONST(7))*1.0D3 C c -------- Transforming the unit of properties ------------- IF (KPA.NE.1.AND.KPA.NE.2) THEN PRESS=PRESS*1.0D5 ENDIF IF (KSTAN.EQ.1) THEN DBHL=DBHL+DBX*CONST(19)+(1.0D0-DBX)*CONST(29) DBSL=DBSL+DBX*CONST(20)+(1.0D0-DBX)*CONST(30) ENDIF IF (KAS.EQ.1) THEN DBVL=DBVL/(DBX*CONST(11)+(1.0D0-DBX)*CONST(21)) DBHL=DBHL/(DBX*CONST(11)+(1.0D0-DBX)*CONST(21)) DBSL=DBSL/(DBX*CONST(11)+(1.0D0-DBX)*CONST(21)) ENDIF c c ============= Real ===================== C P=SNGL(PRESS) V=SNGL(DBVL) H=SNGL(DBHL) S=SNGL(DBSL) RETURN END c ------------------------------------------------------------------- c Temperature,overall composition => Pressure at dew point c ------------------------------------------------------------------- SUBROUTINE SUBPD(J,TT,P,Z,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,Z,V,H,S,TT CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56),PRESSR,PRESS $ ,TEMPR,DBVVR,DBHVR,DBSVR,DBVV,DBHV,DBSV,DBY $ ,FMF024,FMF026,FMF028,FMF040 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C J=0 PRNAME='SUBPD ' c c ------------ Reading properties of H2O and NH3 -------- CALL FMF000(J,CONST,COEFF) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF C c ------- Transforming the unit of properties ---------- T=TT IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ENDIF C c ------------ Checking the range of temperature -------- IF ((T.LT.APPLIM(1)).OR.(T.GT.APPLIM(2))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (Z.LT.0.0E0.OR.Z.GT.1.0E0) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c =========== Double precision ================ C c ----------- real property => reduced property --------- TEMPR=DBLE(T)/CONST(1) DBY=DBLE(Z) IF (KAS.EQ.1) THEN DBY=DBY*CONST(21)/((1.0D0-DBY)*CONST(11)+DBY*CONST(21)) ENDIF C C PRESSR=FMF040(CONST,COEFF,TEMPR,DBY) IF (PRESSR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PRESSR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF DBVVR=FMF024(COEFF,TEMPR,PRESSR,DBY) DBHVR=FMF026(COEFF,TEMPR,PRESSR,DBY) DBSVR=FMF028(COEFF,TEMPR,PRESSR,DBY) C c --------- reduced properties => real properties ----------- c --------- kJ => J ----------------------------------------- PRESS=PRESSR*CONST(2) DBVV=DBVVR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBHV=(DBHVR*CONST(1)*CONST(3)+DBY*CONST(4)+(1.0D0-DBY)*CONST(5)) $ *1.0D3 DBSV=(DBSVR*CONST(3)+DBY*CONST(6)+(1.0D0-DBY)*CONST(7))*1.0D3 C c -------- Transforming the unit of properties ------------- IF (KPA.NE.1.AND.KPA.NE.2) THEN PRESS=PRESS*1.0D5 ENDIF IF (KSTAN.EQ.1) THEN DBHV=DBHV+DBY*CONST(19)+(1.0D0-DBY)*CONST(29) DBSV=DBSV+DBY*CONST(20)+(1.0D0-DBY)*CONST(30) ENDIF IF (KAS.EQ.1) THEN DBVV=DBVV/(DBY*CONST(11)+(1.0D0-DBY)*CONST(21)) DBHV=DBHV/(DBY*CONST(11)+(1.0D0-DBY)*CONST(21)) DBSV=DBSV/(DBY*CONST(11)+(1.0D0-DBY)*CONST(21)) ENDIF C c ============= Real ===================== C P=SNGL(PRESS) V=SNGL(DBVV) H=SNGL(DBHV) S=SNGL(DBSV) RETURN END c ----------------------------------------------------------------- c Subroutine subxy c calculate the composition of gas and liquid region c ----------------------------------------------------------------- SUBROUTINE SUBXY(J,TT,PP,X,Y,VL,VV,HL,HV,SL,SV) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,X,Y,VL,VV,HL,HV,SL,SV,TT,PP CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56),TEMPR,PRESSR $ ,TSTMNR,TSTMXR,PSTAR,PSTWR,PSTCR $ ,DBX,DBY,DBVLR,DBVVR,DBHLR,DBHVR,DBSLR,DBSVR $ ,DBVL,DBVV,DBHL,DBHV,DBSL,DBSV,MXTMPR $ ,FMF023,FMF024,FMF025,FMF026,FMF027,FMF028,FMF037,FMF038 C $ ,FMF041,FMF042 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C J=0 PRNAME='SUBXY ' c c --------- Reading properties of NH3 and H2O --------------- CALL FMF000(J,CONST,COEFF) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF C c ------ Transforming the unit of properties ----------------- T=TT P=PP IF ((KPA.EQ.0).OR.(KPA.GE.4)) THEN P=PP*1.0E-5 ELSEIF (KPA.EQ.1) THEN T=TT+273.15 ELSEIF (KPA.EQ.3) THEN P=PP*1.0E-5 T=TT+273.15 ENDIF C c ------ Checking the range of Temperature and Pressure ------ IF ((T.LT.APPLIM(1)).AND.(T.GT.APPLIM(2))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF ((P.LT.APPLIM(3)).OR.(P.GT.APPLIM(4))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c ============ Double precision ============================== C c ------- real properties => reduced properties ------------- TEMPR=DBLE(T)/CONST(1) PRESSR=DBLE(P)/CONST(2) MXTMPR=DBLE(APPLIM(2))/CONST(1) C c ------- judging phase ---------------- C TSTMNR=FMF042(CONST,COEFF,MXTMPR) C IF (TSTMNR.LE.-1.0D15) THEN C J=-2 C CALL FMF049(J,PRNAME) C RETURN C ELSEIF (TSTMNR.LT.0.0D0) THEN C J=-1 C CALL FMF049(J,PRNAME) C RETURN C ENDIF C TSTMXR=FMF041(CONST,COEFF,MXTMPR) C IF (TSTMXR.LE.-1.0D15) THEN C J=-2 C CALL FMF049(J,PRNAME) C RETURN C ELSEIF (TSTMXR.LT.0.0D0) THEN C J=-1 C CALL FMF049(J,PRNAME) C RETURN C ENDIF TSTMNR=DBLE(APPLIM(5)*1.0E-2) TSTMXR=DBLE(APPLIM(7)*1.0E-2) IF (TEMPR.LT.TSTMNR) THEN PSTAR=FMF038(CONST,COEFF,TEMPR) IF (PSTAR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PSTAR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF PSTWR=FMF037(CONST,COEFF,TEMPR) IF (PSTWR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PSTWR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF IF ((PRESSR.GT.PSTAR).OR.(PRESSR.LT.PSTWR)) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PRESSR.EQ.PSTAR) THEN DBX=1.0D0 DBY=1.0D0 DBVLR=FMF023(COEFF,TEMPR,PRESSR,DBX) DBVVR=FMF024(COEFF,TEMPR,PRESSR,DBY) DBHLR=FMF025(COEFF,TEMPR,PRESSR,DBX) DBHVR=FMF026(COEFF,TEMPR,PRESSR,DBY) DBSLR=FMF027(COEFF,TEMPR,PRESSR,DBX) DBSVR=FMF028(COEFF,TEMPR,PRESSR,DBY) GOTO 1000 ELSEIF (PRESSR.EQ.PSTWR) THEN DBX=0.0D0 DBY=0.0D0 DBVLR=FMF023(COEFF,TEMPR,PRESSR,DBX) DBVVR=FMF024(COEFF,TEMPR,PRESSR,DBY) DBHLR=FMF025(COEFF,TEMPR,PRESSR,DBX) DBHVR=FMF026(COEFF,TEMPR,PRESSR,DBY) DBSLR=FMF027(COEFF,TEMPR,PRESSR,DBX) DBSVR=FMF028(COEFF,TEMPR,PRESSR,DBY) GOTO 1000 ENDIF ELSEIF (TEMPR.EQ.TSTMNR) THEN PSTAR=DBLE(APPLIM(6))/CONST(2) PSTWR=FMF037(CONST,COEFF,TEMPR) IF (PSTWR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PSTWR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF IF (PRESSR.GT.PSTAR.OR.PRESSR.LT.PSTWR) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PRESSR.EQ.PSTAR) THEN DBX=1.0D0 DBY=1.0D0 DBVLR=FMF023(COEFF,TEMPR,PRESSR,DBX) DBVVR=FMF024(COEFF,TEMPR,PRESSR,DBY) DBHLR=FMF025(COEFF,TEMPR,PRESSR,DBX) DBHVR=FMF026(COEFF,TEMPR,PRESSR,DBY) DBSLR=FMF027(COEFF,TEMPR,PRESSR,DBX) DBSVR=FMF028(COEFF,TEMPR,PRESSR,DBY) GOTO 1000 ELSEIF (PRESSR.EQ.PSTWR) THEN DBX=0.0D0 DBY=0.0D0 DBVLR=FMF023(COEFF,TEMPR,PRESSR,DBX) DBVVR=FMF024(COEFF,TEMPR,PRESSR,DBY) DBHLR=FMF025(COEFF,TEMPR,PRESSR,DBX) DBHVR=FMF026(COEFF,TEMPR,PRESSR,DBY) DBSLR=FMF027(COEFF,TEMPR,PRESSR,DBX) DBSVR=FMF028(COEFF,TEMPR,PRESSR,DBY) GOTO 1000 ENDIF ELSEIF ((TEMPR.GT.TSTMNR).AND.(TEMPR.LT.TSTMXR)) THEN PSTWR=FMF037(CONST,COEFF,TEMPR) IF (PSTWR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PSTWR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF PSTCR=(DBLE(APPLIM(6))+(DBLE(APPLIM(8))-DBLE(APPLIM(6))) $ *(TEMPR-TSTMNR)/(TSTMXR-TSTMNR))*1.0E-1 IF (PRESSR.LT.PSTWR) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PRESSR.GE.PSTCR) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PRESSR.EQ.PSTWR) THEN DBX=0.0D0 DBY=0.0D0 DBVLR=FMF023(COEFF,TEMPR,PRESSR,DBX) DBVVR=FMF024(COEFF,TEMPR,PRESSR,DBY) DBHLR=FMF025(COEFF,TEMPR,PRESSR,DBX) DBHVR=FMF026(COEFF,TEMPR,PRESSR,DBY) DBSLR=FMF027(COEFF,TEMPR,PRESSR,DBX) DBSVR=FMF028(COEFF,TEMPR,PRESSR,DBY) GOTO 1000 ENDIF ELSEIF (TEMPR.EQ.TSTMXR) THEN DBX=0.0D0 DBY=0.0D0 DBVLR=FMF023(COEFF,TEMPR,PRESSR,DBX) DBVVR=FMF024(COEFF,TEMPR,PRESSR,DBY) DBHLR=FMF025(COEFF,TEMPR,PRESSR,DBX) DBHVR=FMF026(COEFF,TEMPR,PRESSR,DBY) DBSLR=FMF027(COEFF,TEMPR,PRESSR,DBX) DBSVR=FMF028(COEFF,TEMPR,PRESSR,DBY) GOTO 1000 ELSEIF (TEMPR.GT.TSTMXR) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c ---------- calculating properties ------------- CALL FMF050(J,CONST,COEFF,TEMPR,PRESSR,DBX,DBY) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF DBVLR=FMF023(COEFF,TEMPR,PRESSR,DBX) DBVVR=FMF024(COEFF,TEMPR,PRESSR,DBY) DBHLR=FMF025(COEFF,TEMPR,PRESSR,DBX) DBHVR=FMF026(COEFF,TEMPR,PRESSR,DBY) DBSLR=FMF027(COEFF,TEMPR,PRESSR,DBX) DBSVR=FMF028(COEFF,TEMPR,PRESSR,DBY) 1000 CONTINUE c c --------- reduced properties => real properties ----------- c --------- kJ => J ----------------------------------------- DBVL=DBVLR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBVV=DBVVR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBHL=(DBHLR*CONST(1)*CONST(3)+DBX*CONST(4)+(1.0D0-DBX)*CONST(5)) $ *1.0D3 DBHV=(DBHVR*CONST(1)*CONST(3)+DBY*CONST(4)+(1.0D0-DBY)*CONST(5)) $ *1.0D3 DBSL=(DBSLR*CONST(3)+DBX*CONST(6)+(1.0D0-DBX)*CONST(7))*1.0D3 DBSV=(DBSVR*CONST(3)+DBY*CONST(6)+(1.0D0-DBY)*CONST(7))*1.0D3 C c -------- Transforming the unit of properties ------------- IF (KSTAN.EQ.1) THEN DBHL=DBHL+DBX*CONST(19)+(1.0D0-DBX)*CONST(29) DBHV=DBHV+DBY*CONST(19)+(1.0D0-DBY)*CONST(29) DBSL=DBSL+DBX*CONST(20)+(1.0D0-DBX)*CONST(30) DBSV=DBSV+DBY*CONST(20)+(1.0D0-DBY)*CONST(30) ENDIF IF (KAS.EQ.1) THEN DBVL=DBVL/(DBX*CONST(11)+(1.0D0-DBX)*CONST(21)) DBVV=DBVV/(DBY*CONST(11)+(1.0D0-DBY)*CONST(21)) DBHL=DBHL/(DBX*CONST(11)+(1.0D0-DBX)*CONST(21)) DBHV=DBHV/(DBY*CONST(11)+(1.0D0-DBY)*CONST(21)) DBSL=DBSL/(DBX*CONST(11)+(1.0D0-DBX)*CONST(21)) DBSV=DBSV/(DBY*CONST(11)+(1.0D0-DBY)*CONST(21)) DBX=DBX*CONST(11)/((1.0D0-DBX)*CONST(21)+DBX*CONST(11)) DBY=DBY*CONST(11)/((1.0D0-DBY)*CONST(21)+DBY*CONST(11)) ENDIF C c ============= Real ===================== c X=SNGL(DBX) Y=SNGL(DBY) VL=SNGL(DBVL) VV=SNGL(DBVV) HL=SNGL(DBHL) HV=SNGL(DBHV) SL=SNGL(DBSL) SV=SNGL(DBSV) RETURN END c ------------------------------------------------------------- c Subroutine submix c calculate the properties of ammonia-water mixture c i=1 : T,p,z => v,h,s c i=2 : p,z,h => T,v,s c i=3 : p,z,s => T,v,h c i=4 : p,z,v => T,h,s c ------------------------------------------------------------- SUBROUTINE SUBMIX(I,J,TT,PP,Z,V,H,S) INTEGER I,J,KPA,MESS,KSTAN,KAS,PHASE $ ,FMF043 REAL T,P,Z,V,H,S,TT,PP CHARACTER*6 PRNAME DOUBLE PRECISION DBHR,DBSR,DBVR,TEMPR,PRESSR,DBMOL,DBZ $ ,CONST(1:30),COEFF(1:56),DBH,DBS,DBV,TEMP,DBHLR,DBHVR $ ,DBSLR,DBSVR,DBVLR,DBVVR,QUAL,MOLX,MOLY $ ,FMF023,FMF024,FMF025,FMF026,FMF027,FMF028 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C J=0 PRNAME='SUBMIX' c c --------- Reading properties of NH3 and H2O --------------- CALL FMF000(J,CONST,COEFF) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF C c ************************************************************ c ------ i=1 ; (z,t,p) => (v,h,s) ---------------------------- IF (I.EQ.1) THEN c ************************************************************ c c ------ Transforming the unit of properties ----------------- T=TT P=PP IF ((KPA.EQ.0).OR.(KPA.GE.4)) THEN P=PP*1.0E-5 ELSEIF (KPA.EQ.1) THEN T=TT+273.15 ELSEIF (KPA.EQ.3) THEN P=PP*1.0E-5 T=TT+273.15 ENDIF C c ------ Checking the range of Temperature and Pressure ------ IF ((T.LT.APPLIM(1)).AND.(T.GT.APPLIM(2))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF ((P.LT.APPLIM(3)).OR.(P.GT.APPLIM(4))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (Z.LT.0.0E0.OR.Z.GT.1.0E0) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c ============ Double precision ============================== DBZ=DBLE(Z) IF (KAS.EQ.1) THEN DBMOL=DBZ*CONST(21)/((1.0D0-DBZ)*CONST(11)+DBZ*CONST(21)) ELSE DBMOL=DBZ ENDIF C c ------- real properties => reduced properties ------------- TEMPR=DBLE(T)/CONST(1) PRESSR=DBLE(P)/CONST(2) C c ------- judging the phase of mixture ------------- PHASE=FMF043(CONST,COEFF,TEMPR,PRESSR,DBMOL) IF (PHASE.EQ.-1) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ELSEIF (PHASE.EQ.-2) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (PHASE.EQ.1) THEN DBVR=FMF023(COEFF,TEMPR,PRESSR,DBMOL) DBHR=FMF025(COEFF,TEMPR,PRESSR,DBMOL) DBSR=FMF027(COEFF,TEMPR,PRESSR,DBMOL) ELSEIF (PHASE.EQ.2) THEN CALL FMF050(J,CONST,COEFF,TEMPR,PRESSR,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF QUAL=(DBMOL-MOLX)/(MOLY-MOLX) DBVLR=FMF023(COEFF,TEMPR,PRESSR,MOLX) DBVVR=FMF024(COEFF,TEMPR,PRESSR,MOLY) DBHLR=FMF025(COEFF,TEMPR,PRESSR,MOLX) DBHVR=FMF026(COEFF,TEMPR,PRESSR,MOLY) DBSLR=FMF027(COEFF,TEMPR,PRESSR,MOLX) DBSVR=FMF028(COEFF,TEMPR,PRESSR,MOLY) DBVR=(1.0D0-QUAL)*DBVLR+QUAL*DBVVR DBHR=(1.0D0-QUAL)*DBHLR+QUAL*DBHVR DBSR=(1.0D0-QUAL)*DBSLR+QUAL*DBSVR ELSEIF (PHASE.EQ.3) THEN DBVR=FMF024(COEFF,TEMPR,PRESSR,DBMOL) DBHR=FMF026(COEFF,TEMPR,PRESSR,DBMOL) DBSR=FMF028(COEFF,TEMPR,PRESSR,DBMOL) ENDIF C c --------- reduced properties => real properties ---------------- DBV=DBVR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBH=(DBHR*CONST(1)*CONST(3)+DBMOL*CONST(4) $ +(1.0D0-DBMOL)*CONST(5))*1.0D3 DBS=(DBSR*CONST(3)+DBMOL*CONST(6)+(1.0D0-DBMOL)*CONST(7)) $ *1.0D3 C c ---------- Transfoming the unit of properties ------------ IF (KSTAN.EQ.1) THEN DBH=DBH+DBMOL*CONST(19)+(1.0D0-DBMOL)*CONST(29) DBS=DBS+DBMOL*CONST(20)+(1.0D0-DBMOL)*CONST(30) ENDIF IF (KAS.EQ.1) THEN DBV=DBV/(DBMOL*CONST(11)+(1.0D0-DBMOL)*CONST(21)) DBH=DBH/(DBMOL*CONST(11)+(1.0D0-DBMOL)*CONST(21)) DBS=DBS/(DBMOL*CONST(11)+(1.0D0-DBMOL)*CONST(21)) ENDIF C c --------- real ------------------- V=SNGL(DBV) H=SNGL(DBH) S=SNGL(DBS) RETURN C c ************************************************************ c ------ i=2 ; (z,p,h) => (t,v,s) ---------------------------- ELSEIF (I.EQ.2) THEN c ************************************************************ c c ------ Transforming the unit of properties ----------------- P=PP IF ((KPA.EQ.0).OR.(KPA.GE.3)) THEN P=PP*1.0E-5 ENDIF C c ------ Checking the range of Temperature and Pressure ------ IF ((P.LT.APPLIM(3)).OR.(P.GT.APPLIM(4))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (Z.LT.0.0E0.OR.Z.GT.1.0E0) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c ============ Double precision ============================== c c --------- Reading properties of NH3 and H2O --------------- DBZ=DBLE(Z) C c ------ Transforming the unit of properties ----------------- IF (KAS.EQ.1) THEN DBMOL=DBZ*CONST(21)/((1.0D0-DBZ)*CONST(11)+DBZ*CONST(21)) DBH=DBLE(H)*(CONST(11)*DBMOL+(1.0D0-DBMOL)*CONST(21)) ELSE DBMOL=DBZ DBH=DBLE(H) ENDIF IF (KSTAN.EQ.1) THEN DBH=(DBH-(DBMOL*CONST(19)+(1.0D0-DBMOL)*CONST(29)))*1.0D-3 ELSE DBH=DBH*1.0D-3 ENDIF DBH=DBH-(DBMOL*CONST(4)+(1.0D0-DBMOL)*CONST(5)) C c ------- real properties => reduced properties ------------- PRESSR=DBLE(P)/CONST(2) DBHR=DBH/(CONST(3)*CONST(1)) CALL FMF046(J,CONST,COEFF,DBMOL,PRESSR,DBHR,TEMPR,PHASE) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF IF (PHASE.EQ.1) THEN DBVR=FMF023(COEFF,TEMPR,PRESSR,DBMOL) DBSR=FMF027(COEFF,TEMPR,PRESSR,DBMOL) ELSEIF (PHASE.EQ.2) THEN CALL FMF050(J,CONST,COEFF,TEMPR,PRESSR,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF QUAL=(DBMOL-MOLX)/(MOLY-MOLX) DBVLR=FMF023(COEFF,TEMPR,PRESSR,MOLX) DBVVR=FMF024(COEFF,TEMPR,PRESSR,MOLY) DBSLR=FMF027(COEFF,TEMPR,PRESSR,MOLX) DBSVR=FMF028(COEFF,TEMPR,PRESSR,MOLY) DBVR=(1.0D0-QUAL)*DBVLR+QUAL*DBVVR DBSR=(1.0D0-QUAL)*DBSLR+QUAL*DBSVR ELSEIF (PHASE.EQ.3) THEN DBVR=FMF024(COEFF,TEMPR,PRESSR,DBMOL) DBSR=FMF028(COEFF,TEMPR,PRESSR,DBMOL) ENDIF C c --------- reduced properties => real properties ---------------- DBV=DBVR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBS=(DBSR*CONST(3)+DBMOL*CONST(6)+(1.0D0-DBMOL)*CONST(7)) $ *1.0D3 TEMP=TEMPR*CONST(1) C c ---------- Transfoming the unit of properties ------------ IF (KSTAN.EQ.1) THEN DBS=DBS+DBMOL*CONST(20)+(1.0D0-DBMOL)*CONST(30) ENDIF IF (KAS.EQ.1) THEN DBV=DBV/(DBMOL*CONST(11)+(1.0D0-DBMOL)*CONST(21)) DBS=DBS/(DBMOL*CONST(11)+(1.0D0-DBMOL)*CONST(21)) ENDIF IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TEMP=TEMP-273.15 ENDIF C c --------- real ------------------- V=SNGL(DBV) TT=SNGL(TEMP) S=SNGL(DBS) RETURN C c ************************************************************ c ------ i=3 ; (z,p,s) => (t,v,h) ---------------------------- ELSEIF (I.EQ.3) THEN c ************************************************************ c c ------ Transforming the unit of properties ----------------- P=PP IF ((KPA.EQ.0).OR.(KPA.GE.3)) THEN P=PP*1.0E-5 ENDIF C c ------ Checking the range of Temperature and Pressure ------ IF ((P.LT.APPLIM(3)).OR.(P.GT.APPLIM(4))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (Z.LT.0.0E0.OR.Z.GT.1.0E0) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF c c ============ Double precision ============================== c c --------- Reading properties of NH3 and H2O --------------- DBZ=DBLE(Z) C c ------ Transforming the unit of properties ----------------- IF (KAS.EQ.1) THEN DBMOL=DBZ*CONST(21)/((1.0D0-DBZ)*CONST(11)+DBZ*CONST(21)) DBS=DBLE(S)*(CONST(11)*DBMOL+(1.0D0-DBMOL)*CONST(21)) ELSE DBMOL=DBZ DBS=DBLE(S) ENDIF IF (KSTAN.EQ.1) THEN DBS=(DBS-(DBMOL*CONST(20)+(1.0D0-DBMOL)*CONST(30)))*1.0D-3 ELSE DBS=DBS*1.0D-3 ENDIF DBS=DBS-(DBMOL*CONST(4)+(1.0D0-DBMOL)*CONST(5)) c c ------- real properties => reduced properties ------------- PRESSR=DBLE(P)/CONST(2) DBSR=DBS/CONST(3) CALL FMF047(J,CONST,COEFF,DBMOL,PRESSR,DBSR,TEMPR,PHASE) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF IF (PHASE.EQ.1) THEN DBVR=FMF023(COEFF,TEMPR,PRESSR,DBMOL) DBHR=FMF025(COEFF,TEMPR,PRESSR,DBMOL) ELSEIF (PHASE.EQ.2) THEN CALL FMF050(J,CONST,COEFF,TEMPR,PRESSR,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF QUAL=(DBMOL-MOLX)/(MOLY-MOLX) DBVLR=FMF023(COEFF,TEMPR,PRESSR,MOLX) DBVVR=FMF024(COEFF,TEMPR,PRESSR,MOLY) DBHLR=FMF025(COEFF,TEMPR,PRESSR,MOLX) DBHVR=FMF026(COEFF,TEMPR,PRESSR,MOLY) DBVR=(1.0D0-QUAL)*DBVLR+QUAL*DBVVR DBHR=(1.0D0-QUAL)*DBHLR+QUAL*DBHVR ELSEIF (PHASE.EQ.3) THEN DBVR=FMF024(COEFF,TEMPR,PRESSR,DBMOL) DBHR=FMF026(COEFF,TEMPR,PRESSR,DBMOL) ENDIF C c --------- reduced properties => real properties ---------------- DBV=DBVR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBH=(DBHR*CONST(1)*CONST(3)+DBMOL*CONST(4) $ +(1.0D0-DBMOL)*CONST(5))*1.0D3 TEMP=TEMPR*CONST(1) C c ---------- Transfoming the unit of properties ------------ IF (KSTAN.EQ.1) THEN DBH=DBH+DBMOL*CONST(19)+(1.0D0-DBMOL)*CONST(29) ENDIF IF (KAS.EQ.1) THEN DBV=DBV/(DBMOL*CONST(11)+(1.0D0-DBMOL)*CONST(21)) DBH=DBH/(DBMOL*CONST(11)+(1.0D0-DBMOL)*CONST(21)) ENDIF IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TEMP=TEMP-273.15 ENDIF C c --------- real ------------------- V=SNGL(DBV) TT=SNGL(TEMP) H=SNGL(DBH) RETURN C c ************************************************************ c ------ i=4 ; (z,p,v) => (t,h,s) ---------------------------- ELSEIF (I.EQ.4) THEN c ************************************************************ c c ------ Transforming the unit of properties ----------------- P=PP IF ((KPA.EQ.0).OR.(KPA.GE.3)) THEN P=PP*1.0E-5 ENDIF C c ------ Checking the range of Temperature and Pressure ------ IF ((P.LT.APPLIM(3)).OR.(P.GT.APPLIM(4))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (Z.LT.0.0E0.OR.Z.GT.1.0E0) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c ============ Double precision ============================== c c --------- Reading properties of NH3 and H2O --------------- DBZ=DBLE(Z) C c ------ Transforming the unit of properties ----------------- IF (KAS.EQ.1) THEN DBMOL=DBZ*CONST(21)/((1.0D0-DBZ)*CONST(11)+DBZ*CONST(21)) DBV=DBLE(V)*(CONST(11)*DBMOL+(1.0D0-DBMOL)*CONST(21)) ELSE DBMOL=DBZ DBV=DBLE(V) ENDIF C c ------- real properties => reduced properties ------------- PRESSR=DBLE(P)/CONST(2) DBVR=DBV*CONST(2)*1.0D2/(CONST(3)*CONST(1)) C c ------- judging the phase of mixture ------------- CALL FMF048(J,CONST,COEFF,DBMOL,PRESSR,DBVR,TEMPR,PHASE) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF IF (PHASE.EQ.1) THEN DBHR=FMF025(COEFF,TEMPR,PRESSR,DBMOL) DBSR=FMF027(COEFF,TEMPR,PRESSR,DBMOL) ELSEIF (PHASE.EQ.2) THEN CALL FMF050(J,CONST,COEFF,TEMPR,PRESSR,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF QUAL=(DBMOL-MOLX)/(MOLY-MOLX) DBHLR=FMF025(COEFF,TEMPR,PRESSR,MOLX) DBHVR=FMF026(COEFF,TEMPR,PRESSR,MOLY) DBSLR=FMF027(COEFF,TEMPR,PRESSR,MOLX) DBSVR=FMF028(COEFF,TEMPR,PRESSR,MOLY) DBHR=(1.0D0-QUAL)*DBHLR+QUAL*DBHVR DBSR=(1.0D0-QUAL)*DBSLR+QUAL*DBSVR ELSEIF (PHASE.EQ.3) THEN DBHR=FMF026(COEFF,TEMPR,PRESSR,DBMOL) DBSR=FMF028(COEFF,TEMPR,PRESSR,DBMOL) ENDIF C c --------- reduced properties => real properties ---------------- DBH=(DBHR*CONST(1)*CONST(3)+DBMOL*CONST(4) $ +(1.0D0-DBMOL)*CONST(5))*1.0D3 DBS=(DBSR*CONST(3)+DBMOL*CONST(6)+(1.0D0-DBMOL)*CONST(7)) $ *1.0D3 TEMP=TEMPR*CONST(1) C c ---------- Transfoming the unit of properties ------------ IF (KSTAN.EQ.1) THEN DBH=DBH+DBMOL*CONST(19)+(1.0D0-DBMOL)*CONST(29) DBS=DBS+DBMOL*CONST(20)+(1.0D0-DBMOL)*CONST(30) ENDIF IF (KAS.EQ.1) THEN DBH=DBH/(DBMOL*CONST(11)+(1.0D0-DBMOL)*CONST(21)) DBS=DBS/(DBMOL*CONST(11)+(1.0D0-DBMOL)*CONST(21)) ENDIF IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TEMP=TEMP-273.15 ENDIF C c --------- real ------------------- TT=SNGL(TEMP) H=SNGL(DBH) S=SNGL(DBS) RETURN ELSE J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF END c -------------------------------------------------- c Fundamental properties c -------------------------------------------------- REAL FUNCTION FCM(I,A) INTEGER I,J,KPA,MESS,KSTAN,KAS INTEGER COMP CHARACTER*1 A CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56) COMMON /UNIT/ KPA,MESS,KSTAN,KAS c PRNAME=' FCM ' CALL FMF000(J,CONST,COEFF) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FCM=-1.0E10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FCM=-1.0E20 RETURN ENDIF c IF (I.EQ.1) THEN COMP=10 ELSEIF (I.EQ.2) THEN COMP=20 ELSE CALL FMF049(-2,PRNAME) FCM=-1.0E20 RETURN ENDIF c IF (A.EQ.'M') THEN FCM=SNGL(CONST(COMP+1)) RETURN ELSEIF (A.EQ.'R') THEN FCM=SNGL(CONST(3)/CONST(COMP+1))*1.0E3 RETURN ELSEIF (A.EQ.'T') THEN IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN FCM=SNGL(CONST(COMP+2))-273.15 ELSE FCM=SNGL(CONST(COMP+2)) ENDIF RETURN ELSEIF (A.EQ.'P') THEN IF ((KPA.EQ.1).OR.(KPA.EQ.2)) THEN FCM=SNGL(CONST(COMP+3)) ELSE FCM=SNGL(CONST(COMP+3))*1.0E5 ENDIF RETURN ELSEIF (A.EQ.'V') THEN IF (KAS.EQ.1) THEN FCM=SNGL(CONST(COMP+4)/CONST(COMP+1)) ELSE FCM=SNGL(CONST(COMP+4)) ENDIF RETURN ELSEIF (A.EQ.'E') THEN FCM=SNGL(CONST(COMP+8)) RETURN ELSE CALL FMF049(-2,PRNAME) FCM=-1.0E20 RETURN ENDIF END c ------------------------------------------------------------- c Temperature => Saturated pressure of pure substance c i=1 : Ammonia c i=2 : Water c ------------------------------------------------------------- REAL FUNCTION PSTM(I,TT) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL T,TT CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56),PSATR,PSAT,TEMPR $ ,FMF037,FMF038,NH3LIM,H2OLIM REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C PRNAME=' PSTM ' J=0 c --------- Reading Properties of H2O and NH3 ---------- CALL FMF000(J,CONST,COEFF) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) PSTM=-1.0E10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) PSTM=-1.0E20 RETURN ENDIF C c --------- Transforming The Unit of properties -------- T=TT IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ENDIF C c --------- Checking range of Temperature -------------- IF ((T.LT.APPLIM(1)).OR.(T.GT.APPLIM(2))) THEN J=-2 CALL FMF049(J,PRNAME) PSTM=-1.0E20 RETURN ENDIF C c ---------- real properties => reduced properties ------- TEMPR=DBLE(T)/CONST(1) C c ---------- NH3(ammonia) ----------- IF (I.EQ.1) THEN NH3LIM=DBLE(APPLIM(5))/CONST(1) C IF (T.GT.NH3LIM) THEN ----- modified by R. Akasaka at 26 Jan. 2000 IF (TEMPR.GT.NH3LIM) THEN J=-2 CALL FMF049(J,PRNAME) PSTM=-1.0E20 RETURN ENDIF PSATR=FMF038(CONST,COEFF,TEMPR) IF (PSATR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) PSTM=-1.0E20 RETURN ELSEIF (PSATR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) PSTM=-1.0E10 RETURN ENDIF c ---------- H2O(water) ------------- ELSEIF (I.EQ.2) THEN H2OLIM=DBLE(APPLIM(7))/CONST(1) C IF (T.GT.H2OLIM) THEN ----- modified by R. Akasaka at 26 Jan. 2000 IF (TEMPR.GT.H2OLIM) THEN J=-2 CALL FMF049(J,PRNAME) PSTM=-1.0E20 RETURN ENDIF PSATR=FMF037(CONST,COEFF,TEMPR) IF (PSATR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) PSTM=-1.0E20 RETURN ELSEIF (PSATR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) PSTM=-1.0E10 RETURN ENDIF c ---------- Others ----------------- ELSE J=-2 CALL FMF049(J,PRNAME) PSTM=-1.0E20 RETURN ENDIF C c ---------- reduced properties => real properties -------- PSAT=PSATR*CONST(2) C c ---------- Transforming the unit of property ------------ IF ((KPA.EQ.1).OR.(KPA.EQ.2)) THEN PSTM=SNGL(PSAT) ELSE PSTM=SNGL(PSAT*1.0D5) ENDIF RETURN END c ------------------------------------------------------------- c Pressure => Saturated Temperature of pure substance c i=1 : Ammonia c i=2 : Water c ------------------------------------------------------------- REAL FUNCTION TSPM(I,PP) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL P,PP CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56),TSATR,TSAT,PRESSR $ ,FMF041,FMF042,NH3LIM,H2OLIM REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C PRNAME=' TSPM ' C c --------- Reading Properties of H2O and NH3 ---------- CALL FMF000(J,CONST,COEFF) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) TSPM=-1.0E10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) TSPM=-1.0E20 RETURN ENDIF c --------- Transforming The Unit of properties -------- P=PP IF ((KPA.EQ.0).OR.(KPA.GE.3)) THEN P=PP*1.0E-5 ENDIF c c --------- Checking range of Pressure -------------- IF ((P.LT.APPLIM(3)).OR.(P.GT.APPLIM(4))) THEN J=-2 CALL FMF049(J,PRNAME) TSPM=-1.0E20 RETURN ENDIF C c ---------- real properties => reduced properties ------- PRESSR=DBLE(P)/CONST(2) c c ---------- NH3(ammonia) ----------- IF (I.EQ.1) THEN NH3LIM=DBLE(APPLIM(6))/CONST(2) C IF (P.GT.NH3LIM) THEN ----- modified by R. Akasaka at 14 Jan. 2000 IF (PRESSR.GT.NH3LIM) THEN J=-2 CALL FMF049(J,PRNAME) TSPM=-1.0D20 RETURN ENDIF TSATR=FMF042(CONST,COEFF,PRESSR) IF (TSATR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) TSPM=-1.0E20 RETURN ELSEIF (TSATR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) TSPM=-1.0E10 RETURN ENDIF c ---------- H2O(water) ------------- ELSEIF (I.EQ.2) THEN H2OLIM=DBLE(APPLIM(8))/CONST(2) C IF (P.GT.H2OLIM) THEN ----- modified by R. Akasaka at 14 Jan. 2000 IF (PRESSR.GT.H2OLIM) THEN J=-2 CALL FMF049(J,PRNAME) TSPM=-1.0D20 RETURN ENDIF TSATR=FMF041(CONST,COEFF,PRESSR) IF (TSATR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) TSPM=-1.0E20 RETURN ELSEIF (TSATR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) TSPM=-1.0E10 RETURN ENDIF c ---------- Others ----------------- ELSE J=-2 CALL FMF049(J,PRNAME) TSPM=-1.0E20 RETURN ENDIF C c ---------- reduced properties => real properties -------- TSAT=TSATR*CONST(1) c c ---------- Transforming the unit of property ------------ IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TSPM=SNGL(TSAT-273.15D0) ELSE TSPM=SNGL(TSAT) ENDIF RETURN END c --------------------------------------------- c Seting the common variables KPA and MESS c --------------------------------------------- SUBROUTINE KPAMES(KPAC,MESSC) INTEGER KPAC,MESSC,KPA,MESS,KSTAN,KAS COMMON /UNIT/ KPA,MESS,KSTAN,KAS C KPA=KPAC MESS=MESSC RETURN END c --------------------------------------------------- c Seting the common variables KSTAN and KAS c --------------------------------------------------- SUBROUTINE STNKAS(STAND,KASC) INTEGER KPA,MESS,KSTAN,KAS,STAND,KASC COMMON /UNIT/ KPA,MESS,KSTAN,KAS C KSTAN=STAND KAS=KASC RETURN END c ------------------------------------------------------------------- c Temperature,overall composition => Pressure at bubble point c ------------------------------------------------------------------- SUBROUTINE SUBTB(J,T,PP,Z,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,Z,V,H,S,PP CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56),PRESSR,TEMP $ ,TEMPR,DBVLR,DBHLR,DBSLR,DBVL,DBHL,DBSL,DBX $ ,FMF023,FMF025,FMF027,FMF044 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C J=0 PRNAME='SUBTB ' c c ------------ Reading properties of H2O and NH3 -------- CALL FMF000(J,CONST,COEFF) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF C c ------- Transforming the unit of properties ---------- P=PP IF ((KPA.NE.1).AND.(KPA.NE.2)) THEN P=PP*1.0E-5 ENDIF C c ------------ Checking the range of temperature -------- IF ((P.LT.APPLIM(3)).OR.(P.GT.APPLIM(4))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (Z.LT.0.0E0.OR.Z.GT.1.0E0) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF C c =========== Double precision ================ C c ----------- real property => reduced property --------- PRESSR=DBLE(P)/CONST(2) DBX=DBLE(Z) IF (KAS.EQ.1) THEN DBX=DBX*CONST(21)/((1.0D0-DBX)*CONST(11)+DBX*CONST(21)) ENDIF C C TEMPR=FMF044(CONST,COEFF,PRESSR,DBX) IF (TEMPR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (TEMPR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF DBVLR=FMF023(COEFF,TEMPR,PRESSR,DBX) DBHLR=FMF025(COEFF,TEMPR,PRESSR,DBX) DBSLR=FMF027(COEFF,TEMPR,PRESSR,DBX) c c --------- reduced properties => real properties ----------- c --------- kJ => J ----------------------------------------- TEMP=TEMPR*CONST(1) DBVL=DBVLR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBHL=(DBHLR*CONST(1)*CONST(3)+DBX*CONST(4)+(1.0D0-DBX)*CONST(5)) $ *1.0D3 DBSL=(DBSLR*CONST(3)+DBX*CONST(6)+(1.0D0-DBX)*CONST(7))*1.0D3 C c -------- Transforming the unit of properties ------------- IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TEMP=TEMP-273.15D0 ENDIF IF (KSTAN.EQ.1) THEN DBHL=DBHL+DBX*CONST(19)+(1.0D0-DBX)*CONST(29) DBSL=DBSL+DBX*CONST(20)+(1.0D0-DBX)*CONST(30) ENDIF IF (KAS.EQ.1) THEN DBVL=DBVL/(DBX*CONST(11)+(1.0D0-DBX)*CONST(21)) DBHL=DBHL/(DBX*CONST(11)+(1.0D0-DBX)*CONST(21)) DBSL=DBSL/(DBX*CONST(11)+(1.0D0-DBX)*CONST(21)) ENDIF c c ============= Real ===================== c T=SNGL(TEMP) V=SNGL(DBVL) H=SNGL(DBHL) S=SNGL(DBSL) RETURN END c ------------------------------------------------------------------- c Temperature,overall composition => Pressure at bubble point c ------------------------------------------------------------------- SUBROUTINE SUBTD(J,T,PP,Z,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,Z,V,H,S,PP CHARACTER*6 PRNAME DOUBLE PRECISION CONST(1:30),COEFF(1:56),PRESSR,TEMP $ ,TEMPR,DBVLR,DBHLR,DBSLR,DBVL,DBHL,DBSL,DBY $ ,FMF024,FMF026,FMF028,FMF045 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /UNIT/ KPA,MESS,KSTAN,KAS COMMON /FMFC/ ERRLIM,APPLIM C J=0 PRNAME='SUBTD ' c c ------------ Reading properties of H2O and NH3 -------- CALL FMF000(J,CONST,COEFF) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF C c ------- Transforming the unit of properties ---------- P=PP IF ((KPA.NE.1).AND.(KPA.NE.2)) THEN P=PP*1.0E-5 ENDIF C c ------------ Checking the range of temperature -------- IF ((P.LT.APPLIM(3)).OR.(P.GT.APPLIM(4))) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (Z.LT.0.0E0.OR.Z.GT.1.0E0) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF c c =========== Double precision ================ C c ----------- real property => reduced property --------- PRESSR=DBLE(P)/CONST(2) DBY=DBLE(Z) IF (KAS.EQ.1) THEN DBY=DBY*CONST(21)/((1.0D0-DBY)*CONST(11)+DBY*CONST(21)) ENDIF C C TEMPR=FMF045(CONST,COEFF,PRESSR,DBY) IF (TEMPR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (TEMPR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF DBVLR=FMF024(COEFF,TEMPR,PRESSR,DBY) DBHLR=FMF026(COEFF,TEMPR,PRESSR,DBY) DBSLR=FMF028(COEFF,TEMPR,PRESSR,DBY) c c --------- reduced properties => real properties ----------- c --------- kJ => J ----------------------------------------- TEMP=TEMPR*CONST(1) DBVL=DBVLR*CONST(1)*CONST(3)*0.01D0/CONST(2) DBHL=(DBHLR*CONST(1)*CONST(3)+DBY*CONST(4)+(1.0D0-DBY)*CONST(5)) $ *1.0D3 DBSL=(DBSLR*CONST(3)+DBY*CONST(6)+(1.0D0-DBY)*CONST(7))*1.0D3 C c -------- Transforming the unit of properties ------------- IF ((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TEMP=TEMP-273.15D0 ENDIF IF (KSTAN.EQ.1) THEN DBHL=DBHL+DBY*CONST(19)+(1.0D0-DBY)*CONST(29) DBSL=DBSL+DBY*CONST(20)+(1.0D0-DBY)*CONST(30) ENDIF IF (KAS.EQ.1) THEN DBVL=DBVL/(DBY*CONST(11)+(1.0D0-DBY)*CONST(21)) DBHL=DBHL/(DBY*CONST(11)+(1.0D0-DBY)*CONST(21)) DBSL=DBSL/(DBY*CONST(11)+(1.0D0-DBY)*CONST(21)) ENDIF c c ============= Real ===================== c T=SNGL(TEMP) V=SNGL(DBVL) H=SNGL(DBHL) S=SNGL(DBSL) RETURN END c ----------------------------------------------- c This program transfers the unit of fraction c from kmol to kg. c ----------------------------------------------- REAL FUNCTION AKG(FRMOL) REAL FRMOL DOUBLE PRECISION DMOL,FMF051,DMASS,CONST(1:30),COEFF(1:56) CHARACTER*6 PRNAME INTEGER J c J=0 PRNAME=' AKG ' CALL FMF000(J,CONST,COEFF) IF (J.EQ.-1) THEN AKG=-1.0E10 CALL FMF049(J,PRNAME) RETURN ELSEIF (J.EQ.-2) THEN AKG=-1.0E20 CALL FMF049(J,PRNAME) RETURN ENDIF IF (FRMOL.LT.0.0.OR.FRMOL.GT.1.0) THEN J=-2 AKG=-1.0E20 CALL FMF049(J,PRNAME) RETURN ENDIF DMOL=DBLE(FRMOL) DMASS=FMF051(CONST,DMOL) AKG=SNGL(DMASS) RETURN END c ----------------------------------------------- c This program transfers the unit of fraction c from kg to kmol. c ----------------------------------------------- REAL FUNCTION AKMOL(FRMASS) REAL FRMASS DOUBLE PRECISION DMASS,FMF052,DMOL,CONST(1:30),COEFF(1:56) CHARACTER*6 PRNAME INTEGER J c J=0 PRNAME='AKMOL ' CALL FMF000(J,CONST,COEFF) IF (J.EQ.-1) THEN AKMOL=-1.0E10 CALL FMF049(J,PRNAME) RETURN ELSEIF (J.EQ.-2) THEN AKMOL=-1.0E20 CALL FMF049(J,PRNAME) RETURN ENDIF IF (FRMASS.LT.0.0.OR.FRMASS.GT.1.0) THEN J=-2 AKMOL=-1.0E20 CALL FMF049(J,PRNAME) RETURN ENDIF DMASS=DBLE(FRMASS) DMOL=FMF052(CONST,DMASS) AKMOL=SNGL(DMOL) RETURN END c -------------------------------------------------------------- c Set coefficients for equation of state c -------------------------------------------------------------- SUBROUTINE FMF000(J,CONST,COEFF) INTEGER I,J DOUBLE PRECISION ANH3(1:4),AH2O(1:4),BNH3(1:3),BH2O(1:3) $ ,CNH3(1:4),CH2O(1:4),DNH3(1:3),DH2O(1:3),R0NH3(1:6) $ ,R0H2O(1:6),EXCES(1:16),CONS(1:30),COEFF(1:56),CONST(1:30) CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM C J=0 PRNAME='FMF000' C ----- Condition of Error ----- C ----- log domain error ----- ERRLIM(1)=1.0D-10 c ----- divide by 0 ----- ERRLIM(2)=1.0D-20 c ----- others ----- DO 40 I=3,10 ERRLIM(I)=0.0D0 40 CONTINUE c c ----- Applicable Limit ----- c ----- Minimum Temperature ----- APPLIM(1)=2.3E2 c ----- Maximum Temperature ----- c APPLIM(2)=6.0E2 APPLIM(2)=9.0E2 c ----- Minimum Pressure ----- APPLIM(3)=2.0E-1 c ----- Maximum Pressure ----- APPLIM(4)=1.1E2 c ----- Tsat-max pure NH3 ----- APPLIM(5)=402.33 c ----- Psat-max pure NH3 ----- APPLIM(6)=109.994 c ----- Tsat-max pure H2O ----- c APPLIM(7)=594.8877 APPLIM(7)=647.3 c ----- Psat-max pure H2O ----- APPLIM(8)=193.201 c ----------------------------- APPLIM(9)=0.0E0 APPLIM(10)=0.0E0 C C --------- the coefficients of equation of state ----------- DATA (ANH3(I),I=1,4) /3.971423D-2,-1.790557D-5,-1.308905D-2 $ ,3.752836D-3/ DATA (AH2O(I),I=1,4) /2.748796D-2,-1.016665D-5,-4.452025D-3 $ ,8.389246D-4/ DATA (BNH3(I),I=1,3) /1.634519D1,-6.508119D0,1.448937/ DATA (BH2O(I),I=1,3) /1.214557D1,-1.898065D0,2.911966D-1/ DATA (CNH3(I),I=1,4) /-1.049377D-2,-8.288224D0,-6.647257D2 $ ,-3.045352D3/ DATA (CH2O(I),I=1,4) /2.136131D-2,-3.169291D1,-4.634611D4 $ ,0.0D0/ DATA (DNH3(I),I=1,3) /3.673647D0,9.989629D-2,3.617622D-2/ DATA (DH2O(I),I=1,3) /4.019170D0,-5.175550D-2,1.951939D-2/ DATA (R0NH3(I),I=1,6) /4.878573D0,26.468879D0,1.644773D0 $ ,8.339026D0,3.2252D0,2.0000D0/ DATA (R0H2O(I),I=1,6) /21.821141D0,60.965058D0,5.733498D0 $ ,13.453430D0,5.0705D0,3.0000D0/ DATA (EXCES(I),I=1,16) /-4.1733398D1,2.414D-2,6.702285D0 $ ,-1.1475D-2,6.3608967D1,-6.2490768D1,1.761064D0,8.626D-3 $ ,3.87983D-1,-4.772D-3,-4.648107D0,8.36376D-1 $ ,-3.553627D0,9.04D-4,2.4361723D1,-2.0736547D1/ C DATA (EXCES(I),I=1,16) /-4.626129D1,2.060225D-2,7.292369D0 C $ ,-1.032613D-2,8.074824D1,-8.461214D1,2.452882D1 C $ ,9.598767D-3,-1.475383D0,-5.038107D-3,-9.640398D1 C $ ,1.226973D2,-7.582637D0,6.012445D-4,5.487018D1 C $ ,-7.667596D1/ C C --------- thermodynamic properties -------------------- C - i=1-10 C - 1:the value of the constants TB [K] C - 2:the value of the constants PB [bar] C - 3:the value of the constants R [kj/kmol K] C - 4:the value of the difference from standard state in 'h_{NH3}' C - [kj/kmol] C - 5:the value of the difference from standard state in 'h_{H2O}' C - [kj/kmol] C - 6:the value of the difference from standard state in 's_{NH3}' C - [kj/kmol K] C - 7:the value of the difference from standard state in 's_{H2O}' C - [kj/kmol K] C - 8: C - 9: C - 10: C - i=11-20 NH3(ammonia) C - 1:Molecular weight [kg/kmol] C - 2:Critical temperature [K] C - 3:Critical pressure [bar] C - 4:Critical volume [m3/kmol] C - 5-7:the coefficients of Antoine's equation C - 8:Acentric factor [-] C - 9:the value of the difference from standard state in 'h' C - (at kas=1) [j/kmol] C - 10:the value of the difference from standard state in 's' C - (at kas=1) [j/kmol K] C --------------------------------------------------------------- DATA (CONS(I),I=1,30) /1.0D2,1.0D1,8.314D0,7*0.0D0 $ ,1.703026D1,4.0565D2,1.1278D2,0.07247D0 $ ,1.032794455D1,2.13250D3,-3.298D1,0.2520D0 $ ,3406066.33003218611231D0,16835.5543199961605D0 $ ,1.80153D1,6.4713D2,2.2055D2,0.05595D0 $ ,1.168344455D1,3.81644D3,-4.613D1,0.3449D0 $ ,3609461.855855938028D0,17987.45651817819271D0/ c DO 10 I=1,4 COEFF(I)=ANH3(I) COEFF(I+4)=AH2O(I) COEFF(I+14)=CNH3(I) COEFF(I+18)=CH2O(I) COEFF(I+28)=EXCES(I) COEFF(I+32)=EXCES(I+4) COEFF(I+36)=EXCES(I+8) COEFF(I+40)=EXCES(I+12) 10 CONTINUE DO 20 I=1,3 COEFF(I+8)=BNH3(I) COEFF(I+11)=BH2O(I) COEFF(I+22)=DNH3(I) COEFF(I+25)=DH2O(I) COEFF(I+44)=R0NH3(I) COEFF(I+47)=R0NH3(I+3) COEFF(I+50)=R0H2O(I) COEFF(I+53)=R0H2O(I+3) 20 CONTINUE DO 30 I=1,30 CONST(I)=CONS(I) 30 CONTINUE RETURN END C --------------------------------------------------------- C TEMP_R,PRESSURE_R ===> GIBBS^L_{R,H_2O} C --------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF001(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,GIBB1,GIBB2,GIBB3,GIBB4 $ ,D3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF001' J=0 IF (TEMPR.LT.ERRLIM(1)) THEN J=-3 CALL FMF049(J,PRNAME) FMF001=-1.0D20 RETURN ENDIF D3=3.3333333333333D-1 GIBB1=COEFF(51)-TEMPR*COEFF(53) GIBB2=COEFF(12)*(TEMPR-COEFF(55)) $ +COEFF(13)*0.5D0*(TEMPR**2-COEFF(55)**2) $ +COEFF(14)*(TEMPR**3-COEFF(55)**3)*D3 C GIBB3=-1.0D0*COEFF(12)*TEMPR*(DLOG(TEMPR)-DLOG(COEFF(55))) GIBB3=-1.0D0*COEFF(12)*TEMPR*DLOG(TEMPR/COEFF(55)) $ -COEFF(13)*TEMPR*(TEMPR-COEFF(55)) $ -COEFF(14)*0.5D0*TEMPR*(TEMPR**2-COEFF(55)**2) GIBB4=(COEFF(5)+COEFF(7)*TEMPR+COEFF(8)*TEMPR**2) $ *(PRESSR-COEFF(56))+COEFF(6)*(PRESSR**2-COEFF(56)**2)*0.5D0 FMF001=GIBB1+GIBB2+GIBB3+GIBB4 RETURN END c --------------------------------------------------------- c Temp_R,Pressure_R ===> gibbs^l_{R,NH3} c --------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF002(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,GIBB1,GIBB2 $ ,GIBB3,GIBB4,D3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF002' J=0 IF (TEMPR.LT.ERRLIM(1)) THEN J=-3 CALL FMF049(J,PRNAME) FMF002=-1.0D20 RETURN ENDIF D3=3.3333333333333D-1 GIBB1=COEFF(45)-TEMPR*COEFF(47) GIBB2=COEFF(9)*(TEMPR-COEFF(49)) $ +COEFF(10)*0.5D0*(TEMPR**2-COEFF(49)**2) $ +COEFF(11)*(TEMPR**3-COEFF(49)**3)*D3 C GIBB3=-1.0D0*COEFF(9)*TEMPR*(DLOG(TEMPR)-DLOG(COEFF(49))) GIBB3=-1.0D0*COEFF(9)*TEMPR*DLOG(TEMPR/COEFF(49)) $ -COEFF(10)*TEMPR*(TEMPR-COEFF(49)) $ -COEFF(11)*0.5D0*TEMPR*(TEMPR**2-COEFF(49)**2) GIBB4=(COEFF(1)+COEFF(3)*TEMPR+COEFF(4)*TEMPR**2) $ *(PRESSR-COEFF(50))+COEFF(2)*(PRESSR**2-COEFF(50)**2)*0.5D0 FMF002=GIBB1+GIBB2+GIBB3+GIBB4 RETURN END c -------------------------------------------------------- c Temp_R,Pressure_R ===> g^g_{R,H2O} c -------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF003(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,GIBB1,GIBB2,GIBB3 $ ,GIBB4,GIBB5,GIBB6,D3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF003' J=0 IF (TEMPR.LT.ERRLIM(1).OR.PRESSR.LT.ERRLIM(1)) THEN J=-3 CALL FMF049(J,PRNAME) FMF003=-1.0D20 RETURN ELSEIF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF003=1.0D20 RETURN ENDIF D3=3.333333333333D-1 GIBB1=COEFF(52)-TEMPR*COEFF(54) GIBB2=COEFF(26)*(TEMPR-COEFF(55)) $ +0.5D0*COEFF(27)*(TEMPR**2-COEFF(55)**2) $ +COEFF(28)*(TEMPR**3-COEFF(55)**3)*D3 C GIBB3=-1.0D0*COEFF(26)*TEMPR*(DLOG(TEMPR)-DLOG(COEFF(55))) GIBB3=-1.0D0*COEFF(26)*TEMPR*DLOG(TEMPR/COEFF(55)) $ -COEFF(27)*TEMPR*(TEMPR-COEFF(55)) $ -0.5D0*COEFF(28)*TEMPR*(TEMPR**2-COEFF(55)**2) C GIBB4=TEMPR*(DLOG(PRESSR)-DLOG(COEFF(56))) GIBB4=TEMPR*DLOG(PRESSR/COEFF(56)) GIBB5=COEFF(19)*(PRESSR-COEFF(56)) $ +COEFF(20)*(PRESSR/TEMPR**3-4.0D0*COEFF(56)/COEFF(55)**3 $ +3.0D0*COEFF(56)*TEMPR/COEFF(55)**4) GIBB6=COEFF(21)*(PRESSR/TEMPR**11-1.2D1*COEFF(56)/COEFF(55)**11 $ +1.1D1*COEFF(56)*TEMPR/COEFF(55)**12) $ +COEFF(22)*D3*(PRESSR**3/TEMPR**11-1.2D1*COEFF(56)**3 $ /COEFF(55)**11+1.1D1*COEFF(56)**3*TEMPR/COEFF(55)**12) FMF003=GIBB1+GIBB2+GIBB3+GIBB4+GIBB5+GIBB6 RETURN END c -------------------------------------------------------- c Temp_R,Pressure_R ===> g^g_{R,NH3} c ------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF004(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,GIBB1,GIBB2,GIBB3 $ ,GIBB4,GIBB5,GIBB6,D3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF004' J=0 IF (TEMPR.LT.ERRLIM(1).OR.PRESSR.LT.ERRLIM(1)) THEN J=-3 CALL FMF049(J,PRNAME) FMF004=-1.0D20 RETURN ELSEIF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF004=-1.0D20 RETURN ENDIF D3=3.3333333333333D-1 GIBB1=COEFF(46)-TEMPR*COEFF(48) GIBB2=COEFF(23)*(TEMPR-COEFF(49)) $ +0.5D0*COEFF(24)*(TEMPR**2-COEFF(49)**2) $ +COEFF(25)*(TEMPR**3-COEFF(49)**3)*D3 C GIBB3=-1.0D0*COEFF(23)*TEMPR*(DLOG(TEMPR)-DLOG(COEFF(49))) GIBB3=-1.0D0*COEFF(23)*TEMPR*DLOG(TEMPR/COEFF(49)) $ -COEFF(24)*TEMPR*(TEMPR-COEFF(49)) $ -0.5D0*COEFF(25)*TEMPR*(TEMPR**2-COEFF(49)**2) C GIBB4=TEMPR*(DLOG(PRESSR)-DLOG(COEFF(50))) GIBB4=TEMPR*DLOG(PRESSR/COEFF(50)) GIBB5=COEFF(15)*(PRESSR-COEFF(50)) $ +COEFF(16)*(PRESSR/TEMPR**3-4.0D0*COEFF(50)/COEFF(49)**3 $ +3.0D0*COEFF(50)*TEMPR/COEFF(49)**4) GIBB6=COEFF(17)*(PRESSR/TEMPR**11-1.2D1*COEFF(50)/COEFF(49)**11 $ +1.1D1*COEFF(50)*TEMPR/COEFF(49)**12) $ +COEFF(18)*D3*(PRESSR**3/TEMPR**11-1.2D1*COEFF(50)**3 $ /COEFF(49)**11+1.1D1*COEFF(50)**3*TEMPR/COEFF(49)**12) FMF004=GIBB1+GIBB2+GIBB3+GIBB4+GIBB5+GIBB6 RETURN END c ---------------------------------------------------------- c Temp_R,Pressure_R,x ===> g^l_R c ---------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF005(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,GIBH2O,GIBNH3 $ ,GIBB1,GIBB2,GIBB3 $ ,FMF001,FMF002 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF005' J=0 IF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF005=-1.0D20 RETURN ENDIF GIBH2O=FMF001(COEFF,TEMPR,PRESSR) GIBNH3=FMF002(COEFF,TEMPR,PRESSR) GIBB1=(1.0D0-MOLFR)*GIBH2O+MOLFR*GIBNH3 GIBB3=(COEFF(29)+COEFF(30)*PRESSR+(COEFF(31)+COEFF(32)*PRESSR) $ *TEMPR+COEFF(33)/TEMPR+COEFF(34)/TEMPR**2 $ +(COEFF(35)+COEFF(36)*PRESSR+(COEFF(37)+COEFF(38)*PRESSR) $ *TEMPR+COEFF(39)/TEMPR+COEFF(40)/TEMPR**2) $ *(2.0D0*MOLFR-1.0D0) $ +(COEFF(41)+COEFF(42)*PRESSR+COEFF(43)/TEMPR+COEFF(44) $ /TEMPR**2)*(2.0D0*MOLFR-1.0D0)**2)*MOLFR*(1.0D0-MOLFR) IF ((MOLFR.EQ.1.0D0).OR.(MOLFR.EQ.0.0D0)) THEN FMF005=GIBB1+GIBB3 ELSE GIBB2=TEMPR*((1.0D0-MOLFR)*DLOG(1.0D0-MOLFR)+MOLFR $ *DLOG(MOLFR)) FMF005=GIBB1+GIBB2+GIBB3 ENDIF RETURN END c ------------------------------------------------------------ c Temp_R,Pressure_R,y ===> g^g_R c ------------------------------------------------------------ DOUBLE PRECISION FUNCTION FMF006(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,GIBH2O,GIBNH3 $ ,GIBB1,GIBB2 $ ,FMF003,FMF004 c GIBH2O=FMF003(COEFF,TEMPR,PRESSR) GIBNH3=FMF004(COEFF,TEMPR,PRESSR) GIBB1=(1.0D0-MOLFR)*GIBH2O+MOLFR*GIBNH3 IF ((MOLFR.EQ.1.0D0).OR.(MOLFR.EQ.0.0D0)) THEN FMF006=GIBB1 ELSE GIBB2=TEMPR*((1.0D0-MOLFR)*DLOG(1.0D0-MOLFR) $ +MOLFR*DLOG(MOLFR)) FMF006=GIBB1+GIBB2 ENDIF RETURN END c ------------------------------------------------------------- c Temp_R,Pressure_R ===> (dg^l_{R,H2O}/dp_R) c ------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF007(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR c FMF007=COEFF(5)+COEFF(6)*PRESSR+COEFF(7)*TEMPR+COEFF(8)*TEMPR**2 RETURN END c ------------------------------------------------------------- c Temp_R,Pressure_R ===> (dg^l_{R,NH3}/dp_R) c ------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF008(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR c FMF008=COEFF(1)+COEFF(2)*PRESSR+COEFF(3)*TEMPR+COEFF(4)*TEMPR**2 RETURN END c -------------------------------------------------------------- c Temp_R,Pressure_R,x ===> (dg_{ER}/dp_R) c -------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF009(COEFF,TEMPR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,MOLFR c FMF009=MOLFR*(1.0D0-MOLFR)*(COEFF(30)+COEFF(32)*TEMPR $ +(COEFF(36)+COEFF(38)*TEMPR)*(2.0D0*MOLFR-1.0D0) $ +COEFF(42)*(2.0D0*MOLFR-1.0D0)**2) RETURN END c ---------------------------------------------------------------- c Temp_R,Pressure_R ===> (dg^g_{R,H2O}/dp_R) c ----------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF010(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF010' J=0 IF (DABS(PRESSR).LT.ERRLIM(2).OR.DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF010=-1.0D20 RETURN ENDIF FMF010=TEMPR/PRESSR+COEFF(19)+COEFF(20)/TEMPR**3 $ +COEFF(21)/TEMPR**11+COEFF(22)*PRESSR**2/TEMPR**11 RETURN END c ---------------------------------------------------------------- c Temp_R,Pressure_R ===> (dg^g_{R,NH3}/dp_R) c ----------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF011(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF011' J=0 IF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF011=-1.0D20 RETURN ENDIF FMF011=TEMPR/PRESSR+COEFF(15)+COEFF(16)/TEMPR**3 $ +COEFF(17)/TEMPR**11+COEFF(18)*PRESSR**2/TEMPR**11 RETURN END c --------------------------------------------------------------- c Temp_R,Pressure_R ===> (d(g^l_{R,H2O}/T_R)/dT_R) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF012(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,DG1,DG2,DG3,D3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF012' J=0 IF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF012=-1.0D20 RETURN ENDIF D3=3.333333333333D-1 DG1=-1.0D0*COEFF(51)/TEMPR**2+COEFF(12)*COEFF(55)/TEMPR**2 $ +0.5D0*COEFF(13)*(1.0D0+COEFF(55)**2/TEMPR**2) $ +COEFF(14)*(2.0D0*TEMPR+COEFF(55)**3/TEMPR**2)*D3 DG2=-1.0D0*COEFF(12)/TEMPR-COEFF(13)-COEFF(14)*TEMPR DG3=-1.0D0*COEFF(5)*(PRESSR-COEFF(56))/TEMPR**2 $ -COEFF(6)*(PRESSR**2-COEFF(56)**2)/(2.0D0*TEMPR**2) $ +COEFF(8)*(PRESSR-COEFF(56)) FMF012=DG1+DG2+DG3 RETURN END c --------------------------------------------------------------- c Temp_R,Pressure_R ===> (d(g^l_{R,NH3}/T_R)/dT_R) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF013(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,DG1,DG2,DG3,D3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF013' J=0 IF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF013=-1.0D20 RETURN ENDIF D3=3.3333333333333D-1 DG1=-1.0D0*COEFF(45)/TEMPR**2+COEFF(9)*COEFF(49)/TEMPR**2 $ +0.5D0*COEFF(10)*(1.0D0+COEFF(49)**2/TEMPR**2) $ +COEFF(11)*(2.0D0*TEMPR+COEFF(49)**3/TEMPR**2)*D3 DG2=-1.0D0*COEFF(9)/TEMPR-COEFF(10)-COEFF(11)*TEMPR DG3=-1.0D0*COEFF(1)*(PRESSR-COEFF(50))/TEMPR**2 $ -COEFF(2)*(PRESSR**2-COEFF(50)**2)/(2.0D0*TEMPR**2) $ +COEFF(4)*(PRESSR-COEFF(50)) FMF013=DG1+DG2+DG3 RETURN END c ---------------------------------------------------------------- c Temp_R,Pressure_R ===> (d(g^g_{R,H2O}/T_R)/dT_R) c ---------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF014(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,DG1,DG2,DG3,D3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF014' J=0 IF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF014=-1.0D20 RETURN ENDIF D3=3.3333333333333D-1 DG1=-1.0D0*COEFF(52)/TEMPR**2+COEFF(26)*COEFF(55)/TEMPR**2 $ +0.5D0*COEFF(27)*(1.0D0+COEFF(55)**2/TEMPR**2) $ +COEFF(28)*(2.0D0*TEMPR+COEFF(55)**3/TEMPR**2)*D3 DG2=-1.0D0*COEFF(26)/TEMPR-COEFF(27)-COEFF(28)*TEMPR DG3=-1.0D0*COEFF(19)*(PRESSR-COEFF(56))/TEMPR**2 $ +COEFF(20)*(-4.0D0*PRESSR/TEMPR**5+4.0D0*COEFF(56) $ /(TEMPR**2*COEFF(55)**3))+COEFF(21)*(-1.2D1*PRESSR $ /TEMPR**13+1.2D1*COEFF(56)/(TEMPR**2*COEFF(55)**11)) $ +COEFF(22)*(-4.0D0*PRESSR**3/TEMPR**13+4.0D0*COEFF(56)**3 $ /(TEMPR**2*COEFF(55)**11)) FMF014=DG1+DG2+DG3 RETURN END c ---------------------------------------------------------------- c Temp_R,Pressure_R ===> (d(g^g_{R,NH3}/T_R)/dT_R) c ---------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF015(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,DG1,DG2,DG3,D3 INTEGER J REAL APPLIM(1:10) CHARACTER*6 PRNAME DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF015' J=0 IF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF015=-1.0D20 RETURN ENDIF D3=3.3333333333333D-1 DG1=-1.0D0*COEFF(46)/TEMPR**2+COEFF(23)*COEFF(49)/TEMPR**2 $ +0.5D0*COEFF(24)*(1.0D0+COEFF(49)**2/TEMPR**2) $ +COEFF(25)*(2.0D0*TEMPR+COEFF(49)**3/TEMPR**2)*D3 DG2=-1.0D0*COEFF(23)/TEMPR-COEFF(24)-COEFF(25)*TEMPR DG3=-1.0D0*COEFF(15)*(PRESSR-COEFF(50))/TEMPR**2 $ +COEFF(16)*(-4.0D0*PRESSR/TEMPR**5+4.0D0*COEFF(50) $ /(TEMPR**2*COEFF(49)**3))+COEFF(17)*(-1.2D1*PRESSR $ /TEMPR**13+1.2D1*COEFF(50)/(TEMPR**2*COEFF(49)**11)) $ +COEFF(18)*(-4.0D0*PRESSR**3/TEMPR**13+4.0D0*COEFF(50)**3 $ /(TEMPR**2*COEFF(49)**11)) FMF015=DG1+DG2+DG3 RETURN END c ------------------------------------------------------------ c Temp_R,Pressure_R,x ===> (d(g_{ER}/T_R)/dT_R) c ------------------------------------------------------------ DOUBLE PRECISION FUNCTION FMF016(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,F1,F2,F3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF016' J=0 IF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF016=-1.0D20 RETURN ENDIF F1=-1.0D0*COEFF(29)/TEMPR**2-COEFF(30)*PRESSR/TEMPR**2 $ -2.0D0*COEFF(33)/TEMPR**3-3.0D0*COEFF(34)/TEMPR**4 F2=-1.0D0*COEFF(35)/TEMPR**2-COEFF(36)*PRESSR/TEMPR**2 $ -2.0D0*COEFF(39)/TEMPR**3-3.0D0*COEFF(40)/TEMPR**4 F3=-1.0D0*COEFF(41)/TEMPR**2-COEFF(42)*PRESSR/TEMPR**2 $ -2.0D0*COEFF(43)/TEMPR**3-3.0D0*COEFF(44)/TEMPR**4 FMF016=MOLFR*(1.0D0-MOLFR)*(F1+F2*(2.0D0*MOLFR-1.0D0) $ +F3*(2.0D0*MOLFR-1.0D0)**2) RETURN END c --------------------------------------------------------------- c Temp_R,Pressure_R ===> (dg^l_{R,H2O}/dT_R) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF017(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,DG1,DG2,DG3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF017' J=0 IF (TEMPR.LT.ERRLIM(1)) THEN J=-3 CALL FMF049(J,PRNAME) FMF017=-1.0D20 RETURN ENDIF DG1=-1.0D0*COEFF(53)+COEFF(12)+COEFF(13)*TEMPR+COEFF(14) $ *TEMPR**2 DG2=-1.0D0*COEFF(12)*(DLOG(TEMPR/COEFF(55))+1.0D0) $ -COEFF(13)*(2.0D0*TEMPR-COEFF(55)) $ -0.5D0*COEFF(14)*(3.0D0*TEMPR**2-COEFF(55)**2) DG3=(COEFF(7)+2.0D0*COEFF(8)*TEMPR)*(PRESSR-COEFF(56)) FMF017=DG1+DG2+DG3 RETURN END c --------------------------------------------------------------- c Temp_R,Pressure_R ===> (dg^l_{R,NH3}/dT_R) c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF018(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,DG1,DG2,DG3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF018' J=0 IF (TEMPR.LT.ERRLIM(1)) THEN J=-3 CALL FMF049(J,PRNAME) FMF018=-1.0D20 RETURN ENDIF DG1=-1.0D0*COEFF(47)+COEFF(9)+COEFF(10)*TEMPR+COEFF(11) $ *TEMPR**2 DG2=-1.0D0*COEFF(9)*(DLOG(TEMPR/COEFF(49))+1.0D0) $ -COEFF(10)*(2.0D0*TEMPR-COEFF(49)) $ -0.5D0*COEFF(11)*(3.0D0*TEMPR**2-COEFF(49)**2) DG3=(COEFF(3)+2.0D0*COEFF(4)*TEMPR)*(PRESSR-COEFF(50)) FMF018=DG1+DG2+DG3 RETURN END c -------------------------------------------------------------- c Temp_R,Pressure_R ===> (dg^g_{R,H2O}/T_R) c -------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF019(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,DG1,DG2,DG3,D3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF019' J=0 IF (TEMPR.LT.ERRLIM(1).OR.PRESSR.LT.ERRLIM(1)) THEN J=-3 CALL FMF049(J,PRNAME) FMF019=-1.0D20 RETURN ENDIF D3=3.3333333333333D-1 DG1=-1.0D0*COEFF(54)+COEFF(26)+COEFF(27)*TEMPR $ +COEFF(28)*TEMPR**2 DG2=-1.0D0*COEFF(26)*(DLOG(TEMPR/COEFF(55))+1.0D0) $ -COEFF(27)*(2.0D0*TEMPR-COEFF(55))-0.5D0*COEFF(28) $ *(3.0D0*TEMPR**2-COEFF(55)**2)+DLOG(PRESSR/COEFF(56)) DG3=COEFF(20)*(-3.0D0*PRESSR/TEMPR**4+3.0D0*COEFF(56) $ /COEFF(55)**4)+COEFF(21)*(-1.1D1*PRESSR/TEMPR**12 $ +1.1D1*COEFF(56)/COEFF(55)**12)+COEFF(22) $ *(-1.1D1*PRESSR**3/TEMPR**12+1.1D1*COEFF(56)**3 $ /COEFF(55)**12)*D3 FMF019=DG1+DG2+DG3 RETURN END c -------------------------------------------------------------- c Temp_R,Pressure_R ===> (dg^g_{R,NH3}/T_R) c -------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF020(COEFF,TEMPR,PRESSR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,DG1,DG2,DG3,D3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF020' J=0 IF (TEMPR.LT.ERRLIM(1).OR.PRESSR.LT.ERRLIM(1)) THEN J=-3 CALL FMF049(J,PRNAME) FMF020=-1.0D20 RETURN ELSEIF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF020=-1.0D20 RETURN ENDIF D3=3.3333333333333D-1 DG1=-1.0D0*COEFF(48)+COEFF(23)+COEFF(24)*TEMPR $ +COEFF(25)*TEMPR**2 DG2=-1.0D0*COEFF(23)*(DLOG(TEMPR/COEFF(49))+1.0D0) $ -COEFF(24)*(2.0D0*TEMPR-COEFF(49))-0.5D0*COEFF(25) $ *(3.0D0*TEMPR**2-COEFF(49)**2)+DLOG(PRESSR/COEFF(50)) DG3=COEFF(16)*(-3.0D0*PRESSR/TEMPR**4+3.0D0*COEFF(50) $ /COEFF(49)**4)+COEFF(17)*(-1.1D1*PRESSR/TEMPR**12 $ +1.1D1*COEFF(50)/COEFF(49)**12)+COEFF(18) $ *(-1.1D1*PRESSR**3/TEMPR**12+1.1D1*COEFF(50)**3 $ /COEFF(49)**12)*D3 FMF020=DG1+DG2+DG3 RETURN END c ------------------------------------------------------------- c Temp_R,Pressure_R,x ===> (dg_{ER}/dT_R) c ------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF021(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,F1,F2,F3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF021' J=0 IF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF021=-1.0D20 RETURN ENDIF F1=COEFF(31)+COEFF(32)*PRESSR-COEFF(33)/TEMPR**2 $ -2.0D0*COEFF(34)/TEMPR**3 F2=COEFF(37)+COEFF(38)*PRESSR-COEFF(39)/TEMPR**2 $ -2.0D0*COEFF(40)/TEMPR**3 F3=-1.0D0*COEFF(43)/TEMPR**2-2.0D0*COEFF(44)/TEMPR**3 FMF021=MOLFR*(1.0D0-MOLFR)*(F1+F2*(2.0D0*MOLFR-1.0D0) $ +F3*(2.0D0*MOLFR-1.0D0)**2) RETURN END c -------------------------------------------------------------- c Temp_R,Pressure_R,x ===> (dg_{ER}/dx) c -------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF022(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,F1,F2,F3 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF022' IF (DABS(TEMPR).LT.ERRLIM(2)) THEN J=-4 CALL FMF049(J,PRNAME) FMF022=-1.0D20 RETURN ENDIF F1=COEFF(29)+COEFF(30)*PRESSR+(COEFF(31)+COEFF(32)*PRESSR) $ *TEMPR+COEFF(33)/TEMPR+COEFF(34)/TEMPR**2 F2=COEFF(35)+COEFF(36)*PRESSR+(COEFF(37)+COEFF(38)*PRESSR) $ *TEMPR+COEFF(39)/TEMPR+COEFF(40)/TEMPR**2 F3=COEFF(41)+COEFF(42)*PRESSR+COEFF(43)/TEMPR+COEFF(44)/TEMPR**2 FMF022=F1*(1.0D0-2.0D0*MOLFR)+F2*(-6.0D0*MOLFR**2 $ +6.0D0*MOLFR-1.0D0)+F3*(-1.6D1*MOLFR**3+2.4D1*MOLFR**2 $ -1.0D1*MOLFR+1.0D0) RETURN END c -------------------------------------------------------------- c Temp_R,Pressure_R,x ===> v^l_R c -------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF023(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR $ ,FMF007,FMF008,FMF009 c FMF023=(1.0D0-MOLFR)*FMF007(COEFF,TEMPR,PRESSR)+MOLFR $ *FMF008(COEFF,TEMPR,PRESSR)+FMF009(COEFF,TEMPR,MOLFR) RETURN END c ------------------------------------------------------------------- c Temp_R,Pressure_R,y ===> v^g_R c ------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF024(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR $ ,FMF010,FMF011 c FMF024=(1.0D0-MOLFR)*FMF010(COEFF,TEMPR,PRESSR)+MOLFR $ *FMF011(COEFF,TEMPR,PRESSR) RETURN END c ----------------------------------------------------------------- c Temp_R,Pressure_R,x ===> h^l_R c ----------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF025(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR $ ,FMF012,FMF013,FMF016 c FMF025=-TEMPR**2*((1.0D0-MOLFR)*FMF012(COEFF,TEMPR,PRESSR) $ +MOLFR*FMF013(COEFF,TEMPR,PRESSR) $ +FMF016(COEFF,TEMPR,PRESSR,MOLFR)) RETURN END c ---------------------------------------------------------------- c Temp_R,Pressure_R,y ===> h^g_R c ---------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF026(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR $ ,FMF014,FMF015 c FMF026=-TEMPR**2*((1.0D0-MOLFR)*FMF014(COEFF,TEMPR,PRESSR) $ +MOLFR*FMF015(COEFF,TEMPR,PRESSR)) RETURN END c ---------------------------------------------------------------- c Temp_R,Pressure_R,x ===> s^l_R c ---------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF027(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,EL1,EL2 $ ,FMF017,FMF018,FMF021 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c EL1=1.0D0-ERRLIM(1) EL2=ERRLIM(1) IF ((MOLFR.GT.EL1).OR.(MOLFR.LT.EL2)) THEN FMF027=-((1.0D0-MOLFR)*FMF017(COEFF,TEMPR,PRESSR)+MOLFR $ *FMF018(COEFF,TEMPR,PRESSR) $ +FMF021(COEFF,TEMPR,PRESSR,MOLFR)) ELSE FMF027=-((1.0D0-MOLFR)*FMF017(COEFF,TEMPR,PRESSR)+MOLFR $ *FMF018(COEFF,TEMPR,PRESSR)+(1.0D0-MOLFR) $ *DLOG(1.0D0-MOLFR)+MOLFR*DLOG(MOLFR) $ +FMF021(COEFF,TEMPR,PRESSR,MOLFR)) ENDIF RETURN END c --------------------------------------------------------------- c Temp_R,Pressure_R,y ===> s^g_R c --------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF028(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,EL1,EL2 $ ,FMF019,FMF020 REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c EL1=1.0D0-ERRLIM(1) EL2=ERRLIM(1) IF ((MOLFR.GT.EL1).OR.(MOLFR.LT.EL2)) THEN FMF028=-((1.0D0-MOLFR)*FMF019(COEFF,TEMPR,PRESSR)+MOLFR $ *FMF020(COEFF,TEMPR,PRESSR)) ELSE FMF028=-((1.0D0-MOLFR)*FMF019(COEFF,TEMPR,PRESSR)+MOLFR $ *FMF020(COEFF,TEMPR,PRESSR)+(1.0D0-MOLFR) $ *DLOG(1.0D0-MOLFR)+MOLFR*DLOG(MOLFR)) ENDIF RETURN END c ---------------------------------------------------------------- c Temp_R,Pressure_R,x ===> mu^l_{R,H2O} c ---------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF029(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,EL1,EL2 $ ,FMF005,FMF001,FMF002,FMF022 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF029' J=0 EL1=1.0D0-ERRLIM(1) EL2=ERRLIM(1) IF ((MOLFR.GT.EL1).OR.(MOLFR.LT.EL2)) THEN J=-3 CALL FMF049(J,PRNAME) FMF029=-1.0D20 RETURN ENDIF FMF029=FMF005(COEFF,TEMPR,PRESSR,MOLFR)-MOLFR $ *(FMF002(COEFF,TEMPR,PRESSR)-FMF001(COEFF,TEMPR,PRESSR) C $ +TEMPR*(DLOG(MOLFR)-DLOG(1.0D0-MOLFR)) $ +TEMPR*DLOG(MOLFR/(1.0D0-MOLFR)) $ +FMF022(COEFF,TEMPR,PRESSR,MOLFR)) RETURN END c -------------------------------------------------------------- c Temp_R,Pressure_R,x ===> mu^l_{R,NH3} c -------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF030(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,EL1,EL2 $ ,FMF005,FMF001,FMF002,FMF022 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF030' J=0 EL1=1.0D0-ERRLIM(1) EL2=ERRLIM(1) IF ((MOLFR.GT.EL1).OR.(MOLFR.LT.EL2)) THEN J=-3 CALL FMF049(J,PRNAME) FMF030=-1.0D20 RETURN ENDIF FMF030=FMF005(COEFF,TEMPR,PRESSR,MOLFR)+(1.0D0-MOLFR) $ *(FMF002(COEFF,TEMPR,PRESSR)-FMF001(COEFF,TEMPR,PRESSR) C $ +TEMPR*(DLOG(MOLFR)-DLOG(1.0D0-MOLFR)) $ +TEMPR*DLOG(MOLFR/(1.0D0-MOLFR)) $ +FMF022(COEFF,TEMPR,PRESSR,MOLFR)) RETURN END c ------------------------------------------------------------- c Temp_R,Pressure_R,y ===> mu^g_{R,H2O} c ------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF031(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,EL1,EL2 $ ,FMF006,FMF003,FMF004 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF031' J=0 EL1=1.0D0-ERRLIM(1) EL2=ERRLIM(1) IF (MOLFR.GT.EL1.OR.MOLFR.LT.EL2) THEN J=-3 CALL FMF049(J,PRNAME) FMF031=-1.0D20 RETURN ENDIF FMF031=FMF006(COEFF,TEMPR,PRESSR,MOLFR)-MOLFR $ *(FMF004(COEFF,TEMPR,PRESSR)-FMF003(COEFF,TEMPR,PRESSR) C $ +TEMPR*(DLOG(MOLFR)-DLOG(1.0D0-MOLFR))) $ +TEMPR*DLOG(MOLFR/(1.0D0-MOLFR))) RETURN END c -------------------------------------------------------------- c Temp_R,Pressure_R,y ===> mu^g_{R,NH3} c -------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF032(COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,EL1,EL2 $ ,FMF006,FMF004,FMF003 INTEGER J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF032' J=0 EL1=1.0D0-ERRLIM(1) EL2=ERRLIM(1) IF (MOLFR.GT.EL1.OR.MOLFR.LT.EL2) THEN J=-3 CALL FMF049(J,PRNAME) FMF032=-1.0D20 RETURN ENDIF FMF032=FMF006(COEFF,TEMPR,PRESSR,MOLFR)+(1.0D0-MOLFR) $ *(FMF004(COEFF,TEMPR,PRESSR)-FMF003(COEFF,TEMPR,PRESSR) C $ +TEMPR*(DLOG(MOLFR)-DLOG(1.0D0-MOLFR))) $ +TEMPR*DLOG(MOLFR/(1.0D0-MOLFR))) RETURN END c ------------------------------------------------------------------ c m^l_{R,NH3},Temp_R,Pressure_R ===> x c ------------------------------------------------------------------ DOUBLE PRECISION FUNCTION FMF033(COEFF,TEMPR,PRESSR,CHMPOT) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,CHMPOT,EP1,EP2 $ ,MOLFR0,MOLFR1,MOLFR2,CHMPO0,CHMPO1,CHMPO2,FL1CHM,FL2CHM $ ,FL3CHM,EPMOL,EPCHM $ ,FMF030 INTEGER I c EP1=1.0D-14 EP2=1.0D-14 MOLFR0=1.0D0-1.0D-8 MOLFR1=1.0D-8 CHMPO0=FMF030(COEFF,TEMPR,PRESSR,MOLFR0)-CHMPOT CHMPO1=FMF030(COEFF,TEMPR,PRESSR,MOLFR1)-CHMPOT FL1CHM=CHMPO0*CHMPO1 IF (FL1CHM.GT.0.0D0) THEN FMF033=-1.0D20 WRITE(*,*) ' ERROR IN FMF033. NO RESULT' RETURN ENDIF DO 10 I=1,500 MOLFR2=(MOLFR1+MOLFR0)*0.5D0 CHMPO2=FMF030(COEFF,TEMPR,PRESSR,MOLFR2)-CHMPOT FL2CHM=CHMPO0*CHMPO2 IF (FL2CHM.GT.0.0D0) THEN MOLFR0=MOLFR2 CHMPO0=CHMPO2 ELSE MOLFR1=MOLFR2 CHMPO1=CHMPO2 ENDIF MOLFR2=MOLFR1-CHMPO1*(MOLFR0-MOLFR1)/(CHMPO0-CHMPO1) CHMPO2=FMF030(COEFF,TEMPR,PRESSR,MOLFR2)-CHMPOT FL3CHM=CHMPO0*CHMPO2 IF (FL3CHM.GT.0.0D0) THEN MOLFR0=MOLFR2 CHMPO0=CHMPO2 ELSE MOLFR1=MOLFR2 CHMPO1=CHMPO2 ENDIF EPMOL=(MOLFR1-MOLFR0)**2 EPCHM=CHMPO2**2 IF ((EPMOL.LE.EP1).AND.(EPCHM.LE.EP2)) THEN FMF033=MOLFR2 RETURN ENDIF 10 CONTINUE WRITE(*,*) ' ERROR IN FMF033. NO CONVERGENCE' FMF033=-1.0D10 RETURN END c ------------------------------------------------------------------ c m^g_{R,H2O},Temp_R,Pressure_R ===> 1-y c ------------------------------------------------------------------ DOUBLE PRECISION FUNCTION FMF034(COEFF,TEMPR,PRESSR,CHMPOT) DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,CHMPOT,EP1,EP2 $ ,MOLFR0,MOLFR1,MOLFR2,CHMPO0,CHMPO1,CHMPO2,FL1CHM,FL2CHM $ ,FL3CHM,EPMOL,EPCHM $ ,FMF031 INTEGER I c EP1=1.0D-14 EP2=1.0D-14 MOLFR0=1.0D0-1.0D-8 MOLFR1=1.0D-8 CHMPO0=FMF031(COEFF,TEMPR,PRESSR,MOLFR0)-CHMPOT CHMPO1=FMF031(COEFF,TEMPR,PRESSR,MOLFR1)-CHMPOT FL1CHM=CHMPO0*CHMPO1 IF (FL1CHM.GT.0.0D0) THEN WRITE(*,*) ' ERROR IN FMF034. NO RESULT' FMF034=-1.0D20 RETURN ENDIF DO 10 I=1,500 MOLFR2=(MOLFR1+MOLFR0)*0.5D0 CHMPO2=FMF031(COEFF,TEMPR,PRESSR,MOLFR2)-CHMPOT FL2CHM=CHMPO0*CHMPO2 IF (FL2CHM.GT.0.0D0) THEN MOLFR0=MOLFR2 CHMPO0=CHMPO2 ELSE MOLFR1=MOLFR2 CHMPO1=CHMPO2 ENDIF MOLFR2=MOLFR1-CHMPO1*(MOLFR0-MOLFR1)/(CHMPO0-CHMPO1) CHMPO2=FMF031(COEFF,TEMPR,PRESSR,MOLFR2)-CHMPOT FL3CHM=CHMPO0*CHMPO2 IF (FL3CHM.GT.0.0D0) THEN MOLFR0=MOLFR2 CHMPO0=CHMPO2 ELSE MOLFR1=MOLFR2 CHMPO1=CHMPO2 ENDIF EPMOL=(MOLFR1-MOLFR0)**2 EPCHM=CHMPO2**2 IF ((EPMOL.LE.EP1).AND.(EPCHM.LE.EP2)) THEN FMF034=1.0D0-MOLFR2 RETURN ENDIF 10 CONTINUE FMF034=-1.0D10 WRITE(*,*) ' ERROR IN FMF034. NO CONVERGENCE' RETURN END c --------------------------------------------------------------- c y,Temp_R,Pressure_R ===> x,y_{NH3}+y_{H2O}-1 c --------------------------------------------------------------- SUBROUTINE FMF035(J,COEFF,TEMPR,PRESSR,MOLFR,X,DFY) INTEGER J DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLFR,X,DFY $ ,CHNH3V,MLNH3L,CHH2OL,MLH2OV $ ,FMF029,FMF032,FMF033,FMF034 c c ----------------- m^v_{R,NH3} ----------------- CHNH3V=FMF032(COEFF,TEMPR,PRESSR,MOLFR) c ---------------- x --------------------------- MLNH3L=FMF033(COEFF,TEMPR,PRESSR,CHNH3V) IF (MLNH3L.LE.-1.0D15) THEN J=-2 WRITE(*,*) ' ERROR IN FMF035. X_{HN3} J=',J RETURN ELSEIF(MLNH3L.LE.-9.99D9) THEN J=-1 WRITE(*,*) ' ERROR IN FMF035. X_{HN3} J=',J RETURN ENDIF c --------------- m^l_{R,H2O} ------------------ CHH2OL=FMF029(COEFF,TEMPR,PRESSR,MLNH3L) c --------------- 1-y_{H2O} -------------------- MLH2OV=FMF034(COEFF,TEMPR,PRESSR,CHH2OL) IF (MLH2OV.LE.-1.0D15) THEN J=-2 WRITE(*,*) ' ERROR IN FMF035. 1-Y_{H2O} J=',J RETURN ELSEIF(MLH2OV.LE.-9.99D9) THEN J=-1 WRITE(*,*) ' ERROR IN FMF035. 1-Y_{H2O} J=',J RETURN ENDIF J=0 X=MLNH3L DFY=MOLFR+MLH2OV-1.0D0 RETURN END c ------------------------------------------------------------------ c Temp_R,Pressure_R ===> x,y c ------------------------------------------------------------------ SUBROUTINE FMF036(J,COEFF,TEMPR,PRESSR,MOLX,MOLY) INTEGER I,J DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLX,MOLY,EP1,EP2 $ ,MOLY0,MOLY1,MOLY2,MOLX0,MOLX1,MOLX2,FL1DFY,FL2DFY,FL3DFY $ ,DFY0,DFY1,DFY2,EPMOLY,EPDFY,MOLYQ c EP1=1.0D-14 EP2=1.0D-14 MOLY0=1.0D0-1.0D-7 MOLY1=1.0D-7 CALL FMF035(J,COEFF,TEMPR,PRESSR,MOLY0,MOLX0,DFY0) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF036. Z=1 J=',J RETURN ENDIF CALL FMF035(J,COEFF,TEMPR,PRESSR,MOLY1,MOLX1,DFY1) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF036. Z=0 J=',J RETURN ENDIF FL1DFY=DFY0*DFY1 IF (FL1DFY.GT.0.0D0) THEN WRITE(*,*) ' ERROR IN FMF036. NO RESULT' J=-2 RETURN ENDIF DO 10 I=1,500 MOLY2=(MOLY1+MOLY0)*0.5D0 CALL FMF035(J,COEFF,TEMPR,PRESSR,MOLY2,MOLX2,DFY2) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF036. Z=? J=',J RETURN ENDIF FL2DFY=DFY0*DFY2 IF (FL2DFY.GT.0.0D0) THEN MOLY0=MOLY2 MOLX0=MOLX2 DFY0=DFY2 ELSE MOLY1=MOLY2 MOLX1=MOLX2 DFY1=DFY2 ENDIF MOLY2=MOLY1-DFY1*(MOLY0-MOLY1)/(DFY0-DFY1) CALL FMF035(J,COEFF,TEMPR,PRESSR,MOLY2,MOLX2,DFY2) IF (J.NE.0) THEN WRITE(*,*) ' ERROR IN FMF036. Z=?2 J=',J RETURN ENDIF FL3DFY=DFY0*DFY2 IF (FL3DFY.GT.0.0D0) THEN MOLX0=MOLX2 MOLY0=MOLY2 DFY0=DFY2 ELSE MOLX1=MOLX2 MOLY1=MOLY2 DFY1=DFY2 ENDIF MOLYQ=MOLY2**2 IF (MOLYQ.LE.1) THEN EPMOLY=(MOLY1-MOLY0)**2 EPDFY=DFY2**2 IF ((EPMOLY.LE.EP1).AND.(EPDFY.LE.EP2)) THEN J=0 MOLX=MOLX2 MOLY=MOLY2 RETURN ENDIF ELSE EPMOLY=((MOLY1-MOLY0)/MOLY1)**2 EPDFY=DFY2**2 IF ((EPMOLY.LE.EP1).AND.(EPDFY.LE.EP2)) THEN J=0 MOLX=MOLX2 MOLY=MOLY2 RETURN ENDIF ENDIF 10 CONTINUE WRITE(*,*) ' ERROR IN FMF036. NO CONVERGENCE' J=-1 RETURN END c -------------------------------------------------------------------- c Temp_R ===> Pressure_R(saturated) : Water c -------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF037(CONST,COEFF,TEMPR) DOUBLE PRECISION CONST(1:30),COEFF(1:56),TEMPR,EP,TEMP,P0,P1 $ ,P0R,P1R,P2R,DFG0,DFG1,DFG2,EPDFG,H2OLIM $ ,FMF001,FMF003 INTEGER I,J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF037' J=0 EP=1.0D-14 TEMP=TEMPR*CONST(1) H2OLIM=DBLE(APPLIM(7))/CONST(1)+1.0D-7 IF (TEMPR.GT.H2OLIM) THEN J=2 CALL FMF049(J,PRNAME) FMF037=-1.0D20 RETURN ENDIF P0=DEXP(CONST(25)-CONST(26)/(TEMP+CONST(27))) P1=P0+1.0D-3 P0R=P0/CONST(2) P1R=P1/CONST(2) DFG0=FMF001(COEFF,TEMPR,P0R)-FMF003(COEFF,TEMPR,P0R) DFG1=FMF001(COEFF,TEMPR,P1R)-FMF003(COEFF,TEMPR,P1R) DO 10 I=1,100 P2R=P1R-DFG1*(P1R-P0R)/(DFG1-DFG0) DFG2=FMF001(COEFF,TEMPR,P2R)-FMF003(COEFF,TEMPR,P2R) EPDFG=DFG2**2 IF (EPDFG.LE.EP) THEN FMF037=P2R RETURN ENDIF P0R=P1R DFG0=DFG1 P1R=P2R DFG1=DFG2 10 CONTINUE J=-1 CALL FMF049(J,PRNAME) FMF037=-1.0D10 RETURN END c -------------------------------------------------------------------- c Temp_R ===> Pressure_R(saturated) : Ammonia c -------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF038(CONST,COEFF,TEMPR) DOUBLE PRECISION CONST(1:30),COEFF(1:56),TEMPR,EP,TEMP,P0,P1 $ ,P0R,P1R,P2R,DFG0,DFG1,DFG2,EPDFG,NH3LIM $ ,FMF002,FMF004 INTEGER I,J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF038' J=0 EP=1.0D-14 TEMP=TEMPR*CONST(1) NH3LIM=DBLE(APPLIM(5))/CONST(1) IF (TEMPR.GT.NH3LIM) THEN J=2 CALL FMF049(J,PRNAME) FMF038=-1.0D20 RETURN ENDIF P0=DEXP(CONST(15)-CONST(16)/(TEMP+CONST(17))) P1=P0-1.0D-3 P0R=P0/CONST(2) P1R=P1/CONST(2) DFG0=FMF002(COEFF,TEMPR,P0R)-FMF004(COEFF,TEMPR,P0R) DFG1=FMF002(COEFF,TEMPR,P1R)-FMF004(COEFF,TEMPR,P1R) DO 10 I=1,100 P2R=P1R-DFG1*(P1R-P0R)/(DFG1-DFG0) DFG2=FMF002(COEFF,TEMPR,P2R)-FMF004(COEFF,TEMPR,P2R) EPDFG=DFG2**2 IF (EPDFG.LE.EP) THEN FMF038=P2R RETURN ENDIF P0R=P1R DFG0=DFG1 P1R=P2R DFG1=DFG2 10 CONTINUE J=-1 CALL FMF049(J,PRNAME) FMF038=-1.0D10 RETURN END c ------------------------------------------------------------------- c Temp_R,x ===> Pressure_{R,bubble} c ------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF039(CONST,COEFF,TEMPR,MOLFR) DOUBLE PRECISION CONST(1:30),COEFF(1:56),TEMPR,MOLFR $ ,EP1,EP2,PH2OR,PNH3R,PSATR0,PSATR1,PSATR2,MOLX0,MOLX1,MOLX2 $ ,MOLXD2,MOLYD2,MOLE1,MOLE2,EPMOL,EPPSAT,PSATE,TSTMNR,TSTMXR $ ,MOLXD1,MOLYD1,FLG $ ,FMF037,FMF038 INTEGER I,J,IJ CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF039' J=0 IJ=0 IF ((MOLFR.LT.0.0D0).OR.(MOLFR.GT.1.0D0)) THEN J=-2 CALL FMF049(J,PRNAME) FMF039=-1.0D20 RETURN ELSEIF (MOLFR.LT.1.0D-7) THEN FMF039=FMF037(CONST,COEFF,TEMPR) RETURN ELSEIF (MOLFR.GT.9.999999D-1) THEN FMF039=FMF038(CONST,COEFF,TEMPR) RETURN ENDIF EP1=1.0D-14 EP2=1.0D-14 c ----- Saturated Temperature at 50[bar] ----- c ----- Saturated Temperature at 110[bar] ----- c TSTMNR=3.573918D0 TSTMNR=DBLE(APPLIM(5))/CONST(1) c TSTMXR=5.376596D0 TSTMXR=DBLE(APPLIM(7))/CONST(1) c ------------------------ IF (TEMPR.LT.TSTMNR) THEN PH2OR=FMF037(CONST,COEFF,TEMPR) IF (PH2OR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) C WRITE(*,*) ' ERROR IN FMF039. P_{H2O},OUT OF RANGE' FMF039=-1.0D20 RETURN ELSEIF (PH2OR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) C WRITE(*,*) ' ERROR IN FMF039. P_{H2O},NO CONVERGENCE' FMF039=-1.0D10 RETURN ENDIF PNH3R=FMF038(CONST,COEFF,TEMPR) IF (PNH3R.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) C WRITE(*,*) ' ERROR IN FMF039. P_{NH3},OUT OF RANGE' FMF039=-1.0D20 RETURN ELSEIF (PNH3R.LT.0.0D0) THEN J=-1 C WRITE(*,*) ' ERROR IN FMF039. P_{NH3},NO CONVERGENCE' CALL FMF049(J,PRNAME) FMF039=-1.0D10 RETURN ENDIF PSATR0=PH2OR PSATR1=PNH3R MOLX0=0.0D0-MOLFR MOLX1=1.0D0-MOLFR ELSEIF ((TEMPR.GE.TSTMNR).AND.(TEMPR.LE.TSTMXR)) THEN PSATR1=(22.12-11.304)/(6.473-4.054)*(TEMPR-4.054)+11.304 CALL FMF050(J,CONST,COEFF,TEMPR,PSATR1,MOLXD1,MOLYD1) IF (J.EQ.-2) THEN C WRITE(*,*) ' ERROR IN FMF039. X_{NH3},OUT OF RANGE' CALL FMF049(J,PRNAME) FMF039=-1.0D20 RETURN ELSEIF (J.EQ.-1) THEN C WRITE(*,*) ' ERROR IN FMF039. X_{NH3},NO CONVERGENCE' CALL FMF049(J,PRNAME) FMF039=-1.0D10 RETURN ENDIF MOLX1=MOLXD1-MOLFR PH2OR=FMF037(CONST,COEFF,TEMPR) IF (PH2OR.LE.-1.0D15) THEN J=-2 C WRITE(*,*) ' ERROR IN FMF039. P_{H2O},OUT OF RANGE' CALL FMF049(J,PRNAME) FMF039=-1.0D20 RETURN ELSEIF (PH2OR.LT.0.0D0) THEN J=-1 C WRITE(*,*) ' ERROR IN FMF039. P_{H2O},NO CONVERGENCE' CALL FMF049(J,PRNAME) FMF039=-1.0D10 RETURN ENDIF PSATR0=PH2OR MOLX0=0.0D0-MOLFR FLG=MOLX0*MOLX1 IF (FLG.GE.0.0D0) THEN J=-2 C WRITE(*,*) ' ERROR IN FMF039. OUT OF RANGE' CALL FMF049(J,PRNAME) FMF039=-1.0D20 RETURN ENDIF ELSEIF (TEMPR.GT.TSTMXR) THEN J=-2 CALL FMF049(J,PRNAME) FMF039=-1.0D20 RETURN ENDIF DO 10 I=1,100 PSATR2=(PSATR0+PSATR1)*0.5D0 CALL FMF050(J,CONST,COEFF,TEMPR,PSATR2,MOLXD2,MOLYD2) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FMF039=-1.0D10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FMF039=-1.0D20 RETURN ENDIF MOLX2=MOLXD2-MOLFR MOLE1=MOLX0*MOLX2 IF (MOLE1.GT.0.0D0) THEN MOLX0=MOLX2 PSATR0=PSATR2 ELSE MOLX1=MOLX2 PSATR1=PSATR2 ENDIF PSATR2=PSATR1-MOLX1*(PSATR1-PSATR0)/(MOLX1-MOLX0) CALL FMF050(J,CONST,COEFF,TEMPR,PSATR2,MOLXD2,MOLYD2) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FMF039=-1.0D10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FMF039=-1.0D20 RETURN ENDIF MOLX2=MOLXD2-MOLFR MOLE2=MOLX0*MOLX2 IF (MOLE2.GT.0.0D0) THEN MOLX0=MOLX2 PSATR0=PSATR2 ELSE MOLX1=MOLX2 PSATR1=PSATR2 ENDIF PSATE=PSATR1**2 IF (PSATE.LE.1.0D0) THEN EPPSAT=(PSATR1-PSATR0)**2 EPMOL=MOLX2**2 IF ((EPPSAT.LE.EP1).AND.(EPMOL.LE.EP2)) THEN FMF039=PSATR2 RETURN ENDIF ELSE EPPSAT=((PSATR1-PSATR0)/PSATR1)**2 EPMOL=MOLX2**2 IF ((EPPSAT.LE.EP1).AND.(EPMOL.LE.EP2)) THEN FMF039=PSATR2 RETURN ENDIF ENDIF 10 CONTINUE IJ=-1 CALL FMF049(IJ,PRNAME) FMF039=-1.0D10 RETURN END c ------------------------------------------------------------------- c Temp_R,y ===> Pressure_{R,dew} c ------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF040(CONST,COEFF,TEMPR,MOLFR) DOUBLE PRECISION CONST(1:30),COEFF(1:56),TEMPR,MOLFR $ ,EP1,EP2,PH2OR,PNH3R,PSATR0,PSATR1,PSATR2,MOLY0,MOLY1,MOLY2 $ ,MOLXD2,MOLYD2,MOLE1,MOLE2,EPMOL,EPPSAT,PSATE,TSTMNR,TSTMXR $ ,MOLXD1,MOLYD1,FLG $ ,FMF037,FMF038 INTEGER I,J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF040' J=0 IF ((MOLFR.LT.0.0D0).OR.(MOLFR.GT.1.0D0)) THEN J=-2 CALL FMF049(J,PRNAME) FMF040=-1.0D20 RETURN ELSEIF (MOLFR.LT.1.0D-7) THEN FMF040=FMF037(CONST,COEFF,TEMPR) RETURN ELSEIF (MOLFR.GT.9.999999D-1) THEN FMF040=FMF038(CONST,COEFF,TEMPR) RETURN ENDIF EP1=1.0D-14 EP2=1.0D-14 c ----- Saturated Temperature at 50[bar] ----- c ----- Saturated Temperature at 110[bar] ----- c TSTMNR=3.573918D0 TSTMNR=DBLE(APPLIM(5))/CONST(1) c TSTMXR=5.376596D0 TSTMXR=DBLE(APPLIM(7))/CONST(1) c ------------------------ IF (TEMPR.LT.TSTMNR) THEN PH2OR=FMF037(CONST,COEFF,TEMPR) IF (PH2OR.LE.-1.0D15) THEN J=-2 C WRITE(*,*) ' ERROR IN FMF040. P_{H2O},OUT OF RANGE' CALL FMF049(J,PRNAME) FMF040=-1.0D20 RETURN ELSEIF (PH2OR.LT.0.0D0) THEN J=-1 C WRITE(*,*) ' ERROR IN FMF040. P_{H2O},NO CONVERGENCE' CALL FMF049(J,PRNAME) FMF040=-1.0D10 RETURN ENDIF PNH3R=FMF038(CONST,COEFF,TEMPR) IF (PNH3R.LE.-1.0D15) THEN J=-2 C WRITE(*,*) ' ERROR IN FMF040. P_{NH3},OUT OF RANGE' CALL FMF049(J,PRNAME) FMF040=-1.0D20 RETURN ELSEIF (PNH3R.LT.0.0D0) THEN J=-1 C WRITE(*,*) ' ERROR IN FMF040. P_{NH3},NO CONVERGENCE' CALL FMF049(J,PRNAME) FMF040=-1.0D10 RETURN ENDIF PSATR0=PH2OR PSATR1=PNH3R MOLY0=0.0D0-MOLFR MOLY1=1.0D0-MOLFR ELSEIF ((TEMPR.GE.TSTMNR).AND.(TEMPR.LE.TSTMXR)) THEN PSATR1=(22.12-11.304)/(6.473-4.054)*(TEMPR-4.054)+11.304 CALL FMF050(J,CONST,COEFF,TEMPR,PSATR1,MOLXD1,MOLYD1) IF (J.EQ.-2) THEN C WRITE(*,*) ' ERROR IN FMF039. X_{NH3},OUT OF RANGE' CALL FMF049(J,PRNAME) FMF040=-1.0D20 RETURN ELSEIF (J.EQ.-1) THEN WRITE(*,*) ' ERROR IN FMF039. X_{NH3},NO CONVERGENCE' CALL FMF049(J,PRNAME) FMF040=-1.0D10 RETURN ENDIF MOLY1=MOLYD1-MOLFR PH2OR=FMF037(CONST,COEFF,TEMPR) IF (PH2OR.LE.-1.0D15) THEN J=-2 C WRITE(*,*) ' ERROR IN FMF040. P_{H2O},OUT OF RANGE' CALL FMF049(J,PRNAME) FMF040=-1.0D20 RETURN ELSEIF (PH2OR.LT.0.0D0) THEN J=-1 C WRITE(*,*) ' ERROR IN FMF040. P_{H2O},NO CONVERGENCE' CALL FMF049(J,PRNAME) FMF040=-1.0D10 RETURN ENDIF PSATR0=PH2OR MOLY0=0.0D0-MOLFR FLG=MOLY0*MOLY1 IF (FLG.GE.0.0D0) THEN J=-2 CALL FMF049(J,PRNAME) C WRITE(*,*) ' ERROR IN FMF040. OUT OF RANGE' FMF040=-1.0D20 RETURN ENDIF ELSEIF (TEMPR.GT.TSTMXR) THEN J=-2 CALL FMF049(J,PRNAME) FMF040=-1.0D20 RETURN ENDIF DO 10 I=1,100 PSATR2=(PSATR0+PSATR1)*0.5D0 CALL FMF050(J,CONST,COEFF,TEMPR,PSATR2,MOLXD2,MOLYD2) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FMF040=-1.0D10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FMF040=-1.0D20 RETURN ENDIF MOLY2=MOLYD2-MOLFR MOLE1=MOLY0*MOLY2 IF (MOLE1.GT.0.0D0) THEN MOLY0=MOLY2 PSATR0=PSATR2 ELSE MOLY1=MOLY2 PSATR1=PSATR2 ENDIF PSATR2=PSATR1-MOLY1*(PSATR1-PSATR0)/(MOLY1-MOLY0) CALL FMF050(J,CONST,COEFF,TEMPR,PSATR2,MOLXD2,MOLYD2) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FMF040=-1.0D10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FMF040=-1.0D20 RETURN ENDIF MOLY2=MOLYD2-MOLFR MOLE2=MOLY0*MOLY2 IF (MOLE2.GT.0.0D0) THEN MOLY0=MOLY2 PSATR0=PSATR2 ELSE MOLY1=MOLY2 PSATR1=PSATR2 ENDIF PSATE=PSATR1**2 IF (PSATE.LE.1.0D0) THEN EPPSAT=(PSATR1-PSATR0)**2 EPMOL=MOLY2**2 IF ((EPPSAT.LE.EP1).AND.(EPMOL.LE.EP2)) THEN FMF040=PSATR2 RETURN ENDIF ELSE EPPSAT=((PSATR1-PSATR0)/PSATR1)**2 EPMOL=MOLY2**2 IF ((EPPSAT.LE.EP1).AND.(EPMOL.LE.EP2)) THEN FMF040=PSATR2 RETURN ENDIF ENDIF 10 CONTINUE J=-1 C WRITE(*,*) ' ERROR IN FMF040. NO CONVERGENCE' CALL FMF049(J,PRNAME) FMF040=-1.0D10 RETURN END c -------------------------------------------------------------------- c Pressure_R ===> Temp_R(saturated) : Water c -------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF041(CONST,COEFF,PRESSR) DOUBLE PRECISION CONST(1:30),COEFF(1:56),PRESSR,EP,PRESS,T0,T1 $ ,T0R,T1R,T2R,DFG0,DFG1,DFG2,EPDFG $ ,FMF001,FMF003 INTEGER I,J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF041' J=0 IF (PRESSR.LT.0.0D0) THEN J=-2 C WRITE(*,*) ' ERROR IN FMF041. PRESSR LESS THAN ZERO' CALL FMF049(J,PRNAME) FMF041=-1.0D20 RETURN ENDIF EP=1.0D-14 PRESS=PRESSR*CONST(2) T0=CONST(26)/(CONST(25)-DLOG(PRESS))-CONST(27) T1=T0+1.0D-7 T0R=T0/CONST(1) T1R=T1/CONST(1) DFG0=FMF001(COEFF,T0R,PRESSR)-FMF003(COEFF,T0R,PRESSR) DFG1=FMF001(COEFF,T1R,PRESSR)-FMF003(COEFF,T1R,PRESSR) DO 10 I=1,100 T2R=T1R-DFG1*(T1R-T0R)/(DFG1-DFG0) DFG2=FMF001(COEFF,T2R,PRESSR)-FMF003(COEFF,T2R,PRESSR) EPDFG=DFG2**2 IF (EPDFG.LE.EP) THEN FMF041=T2R RETURN ENDIF T0R=T1R DFG0=DFG1 T1R=T2R DFG1=DFG2 10 CONTINUE J=-1 CALL FMF049(J,PRNAME) FMF041=-1.0D10 C WRITE(*,*) ' ERROR IN FMF041. NO CONVERGENCE' RETURN END c -------------------------------------------------------------------- c Pressure_R ===> Temp_R(saturated) : Ammonia c -------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF042(CONST,COEFF,PRESSR) DOUBLE PRECISION CONST(1:30),COEFF(1:56),PRESSR,EP,PRESS,T0,T1 $ ,T0R,T1R,T2R,DFG0,DFG1,DFG2,EPDFG $ ,FMF002,FMF004 INTEGER I,J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF042' J=0 EP=1.0D-14 IF (PRESSR.LT.0.0D0) THEN J=-2 c WRITE(*,*) ' ERROR IN FMF042. PRESSR LESS THAN ZERO' CALL FMF049(J,PRNAME) FMF042=-1.0D20 RETURN ENDIF PRESS=PRESSR*CONST(2) T0=CONST(16)/(CONST(15)-DLOG(PRESS))-CONST(17) T1=T0+1.0D-7 T0R=T0/CONST(1) T1R=T1/CONST(1) DFG0=FMF002(COEFF,T0R,PRESSR)-FMF004(COEFF,T0R,PRESSR) DFG1=FMF002(COEFF,T1R,PRESSR)-FMF004(COEFF,T1R,PRESSR) DO 10 I=1,100 T2R=T1R-DFG1*(T1R-T0R)/(DFG1-DFG0) DFG2=FMF002(COEFF,T2R,PRESSR)-FMF004(COEFF,T2R,PRESSR) EPDFG=DFG2**2 IF (EPDFG.LE.EP) THEN FMF042=T2R RETURN ENDIF T0R=T1R DFG0=DFG1 T1R=T2R DFG1=DFG2 10 CONTINUE J=-1 CALL FMF049(J,PRNAME) FMF042=-1.0D10 RETURN END c ----------------------------------------------------------- c Temp_R,Pressure_R,Mole fraction => Region c FMF043 = 1 : liquid c 2 : wet c 3 : gas c *4 : super critical c -1 : no convergence c -2 : out of range c ----------------------------------------------------------- INTEGER FUNCTION FMF043(CONST,COEFF,TEMPR,PRESSR,MOLFR) DOUBLE PRECISION CONST(1:30),COEFF(1:56),TEMPR,PRESSR $ ,MOLFR,PSATHR,PSATNR,MOLX,MOLY,LIMNH3,LIMH2O,PLMH2O $ ,FMF037,FMF038,FMF041 CHARACTER*6 PRNAME INTEGER J REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM C PRNAME='FMF043' J=0 C LIMNH3=DBLE(APPLIM(5))/CONST(1) PLMH2O=110.0D0/CONST(2) LIMH2O=FMF041(CONST,COEFF,PLMH2O) IF (TEMPR.LE.LIMNH3) THEN PSATHR=FMF037(CONST,COEFF,TEMPR) IF (PSATHR.LT.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) FMF043=-2 RETURN ELSEIF (PSATHR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) FMF043=-1 RETURN ENDIF PSATNR=FMF038(CONST,COEFF,TEMPR) IF (PSATNR.LT.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) FMF043=-2 RETURN ELSEIF (PSATNR.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) FMF043=-1 RETURN ENDIF IF (PRESSR.GE.PSATNR) THEN FMF043=1 RETURN ELSEIF (PRESSR.LE.PSATHR) THEN FMF043=3 RETURN ELSE CALL FMF050(J,CONST,COEFF,TEMPR,PRESSR,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) FMF043=J RETURN ENDIF IF (MOLFR.LE.MOLX) THEN FMF043=1 RETURN ELSEIF ((MOLFR.GT.MOLX).AND.(MOLFR.LT.MOLY)) THEN FMF043=2 RETURN ELSE FMF043=3 RETURN ENDIF ENDIF ELSEIF (TEMPR.LE.LIMH2O) THEN PSATHR=FMF037(CONST,COEFF,TEMPR) IF (PRESSR.LE.PSATHR) THEN FMF043=3 RETURN ELSE CALL FMF050(J,CONST,COEFF,TEMPR,PRESSR,MOLX,MOLY) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) FMF043=J RETURN ENDIF IF (MOLFR.LE.MOLX) THEN FMF043=1 RETURN ELSEIF ((MOLFR.GT.MOLX).AND.(MOLFR.LT.MOLY)) THEN FMF043=2 RETURN ELSE FMF043=3 RETURN ENDIF ENDIF ELSE FMF043=3 RETURN ENDIF END c ------------------------------------------------------------------- c Pressure_R,x ===> Temp_{R,bubble} c ------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF044(CONST,COEFF,PRESSR,MOLFR) DOUBLE PRECISION CONST(1:30),COEFF(1:56),PRESSR,MOLFR $ ,EP1,EP2,PS1R,PS2R,PS3R,TSAT0R,TSAT1R,TSAT2R,MOLX0,MOLX1 $ ,MOLX2,MOLX2P,MOLY2P,FLMOL,EPMOL,EPTSAT,TSATE $ ,MOLX1P,MOLY1P,TMXNH3,TMXH2O,PS4R $ ,FMF037,FMF038,FMF041,FMF042 INTEGER I,J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF044' J=0 IF ((MOLFR.LT.0.0D0).OR.(MOLFR.GT.1.0D0)) THEN J=-2 CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ELSEIF (MOLFR.LT.1.0D-7) THEN FMF044=FMF041(CONST,COEFF,PRESSR) RETURN ELSEIF (MOLFR.GT.9.999999D-1) THEN FMF044=FMF042(CONST,COEFF,PRESSR) RETURN ENDIF EP1=1.0D-14 EP2=1.0D-14 c ----- Pst_{H2O}(T=230K) ----- PS1R=FMF037(CONST,COEFF,2.3D0) c ----- Pst_{NH3}(T=230K) ----- PS2R=FMF038(CONST,COEFF,2.3D0) c ----- Pst_{NH3}(T=TmaxNH3) ----- TMXNH3=DBLE(APPLIM(5))/CONST(1) PS3R=FMF038(CONST,COEFF,TMXNH3) c ----- Pst_{H2O}(T=TmaxH2O) ----- TMXH2O=DBLE(APPLIM(7))/CONST(1) PS4R=FMF037(CONST,COEFF,TMXH2O) c ----------------------------- IF (PRESSR.LT.PS1R) THEN J=-2 CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ELSEIF ((PRESSR.GE.PS1R).AND.(PRESSR.LT.PS2R)) THEN TSAT1R=2.3D0 CALL FMF050(J,CONST,COEFF,TSAT1R,PRESSR,MOLX1P,MOLY1P) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FMF044=-1.0D10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ENDIF MOLX1=MOLX1P-MOLFR MOLX0=0.0D0-MOLFR FLMOL=MOLX0*MOLX1 IF (FLMOL.GT.0.0D0) THEN J=-2 CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ELSEIF (FLMOL.EQ.0.0D0) THEN IF (MOLX0.EQ.0.0D0) THEN FMF044=TSAT0R RETURN ELSEIF (MOLX1.EQ.0.0D0) THEN FMF044=TSAT1R RETURN ENDIF ENDIF TSAT0R=FMF041(CONST,COEFF,PRESSR) IF (TSAT0R.LT.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ELSEIF (TSAT0R.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) FMF044=-1.0D10 RETURN ENDIF ELSEIF ((PRESSR.GE.PS2R).AND.(PRESSR.LE.PS3R)) THEN MOLX0=0.0D0-MOLFR MOLX1=1.0D0-MOLFR TSAT0R=FMF041(CONST,COEFF,PRESSR) IF (TSAT0R.LT.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ELSEIF (TSAT0R.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) FMF044=-1.0D10 RETURN ENDIF TSAT1R=FMF042(CONST,COEFF,PRESSR) IF (TSAT1R.LT.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ELSEIF (TSAT1R.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) FMF044=-1.0D10 RETURN ENDIF ELSEIF (PRESSR.GT.PS3R.AND.PRESSR.LE.PS4R) THEN TSAT1R=(TMXH2O-TMXNH3)/(PS4R-PS3R)*(PRESSR-PS3R)+TMXNH3 CALL FMF050(J,CONST,COEFF,TSAT1R,PRESSR,MOLX1P,MOLY1P) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FMF044=-1.0D10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ENDIF MOLX1=MOLX1P-MOLFR MOLX0=0.0D0-MOLFR TSAT0R=FMF041(CONST,COEFF,PRESSR) IF (TSAT0R.LT.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ELSEIF (TSAT0R.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) FMF044=-1.0D10 RETURN ENDIF FLMOL=MOLX0*MOLX1 IF (FLMOL.GT.0.0D0) THEN J=-2 CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ELSEIF (FLMOL.EQ.0.0D0) THEN IF (MOLX0.EQ.0.0D0) THEN FMF044=TSAT0R RETURN ELSEIF (MOLX1.EQ.0.0D0) THEN FMF044=TSAT1R RETURN ENDIF ENDIF ELSEIF (PRESSR.GT.PS4R) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ENDIF DO 10 I=1,100 TSAT2R=(TSAT0R+TSAT1R)*0.5D0 CALL FMF050(J,CONST,COEFF,TSAT2R,PRESSR,MOLX2P,MOLY2P) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FMF044=-1.0D10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ENDIF MOLX2=MOLX2P-MOLFR FLMOL=MOLX0*MOLX2 IF (FLMOL.GT.0.0D0) THEN MOLX0=MOLX2 TSAT0R=TSAT2R ELSE MOLX1=MOLX2 TSAT1R=TSAT2R ENDIF TSAT2R=TSAT1R-MOLX1*(TSAT1R-TSAT0R)/(MOLX1-MOLX0) CALL FMF050(J,CONST,COEFF,TSAT2R,PRESSR,MOLX2P,MOLY2P) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FMF044=-1.0D10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FMF044=-1.0D20 RETURN ENDIF MOLX2=MOLX2P-MOLFR FLMOL=MOLX0*MOLX2 IF (FLMOL.GT.0.0D0) THEN MOLX0=MOLX2 TSAT0R=TSAT2R ELSE MOLX1=MOLX2 TSAT1R=TSAT2R ENDIF TSATE=TSAT1R**2 IF (TSATE.LE.1.0D0) THEN EPTSAT=(TSAT1R-TSAT0R)**2 EPMOL=MOLX2**2 IF ((EPTSAT.LE.EP1).AND.(EPMOL.LE.EP2)) THEN FMF044=TSAT2R RETURN ENDIF ELSE EPTSAT=((TSAT1R-TSAT0R)/TSAT1R)**2 EPMOL=MOLX2**2 IF ((EPTSAT.LE.EP1).AND.(EPMOL.LE.EP2)) THEN FMF044=TSAT2R RETURN ENDIF ENDIF 10 CONTINUE J=-1 CALL FMF049(J,PRNAME) FMF044=-1.0D10 RETURN END c ------------------------------------------------------------------- c Pressure_R,y ===> Temp_{R,dew} c ------------------------------------------------------------------- DOUBLE PRECISION FUNCTION FMF045(CONST,COEFF,PRESSR,MOLFR) DOUBLE PRECISION CONST(1:30),COEFF(1:56),PRESSR,MOLFR $ ,EP1,EP2,PS1R,PS2R,PS3R,TSAT0R,TSAT1R,TSAT2R,MOLY0,MOLY1 $ ,MOLY2,MOLX2P,MOLY2P,FLMOL,EPMOL,EPTSAT,TSATE $ ,MOLX1P,MOLY1P,TMXNH3,TMXH2O,PS4R $ ,FMF041,FMF042,FMF037,FMF038 INTEGER I,J CHARACTER*6 PRNAME REAL APPLIM(1:10) DOUBLE PRECISION ERRLIM(1:10) COMMON /FMFC/ ERRLIM,APPLIM c PRNAME='FMF045' J=0 IF ((MOLFR.LT.0.0D0).OR.(MOLFR.GT.1.0D0)) THEN J=-2 CALL FMF049(J,PRNAME) FMF045=-1.0D20 RETURN ELSEIF (MOLFR.LT.1.0D-7) THEN FMF045=FMF041(CONST,COEFF,PRESSR) RETURN ELSEIF (MOLFR.GT.9.999999D-1) THEN FMF045=FMF042(CONST,COEFF,PRESSR) RETURN ENDIF EP1=1.0D-14 EP2=1.0D-14 c ----- Pst_{H2O}(T=230K) ----- PS1R=FMF037(CONST,COEFF,2.3D0) c ----- Pst_{NH3}(T=230K) ----- PS2R=FMF038(CONST,COEFF,2.3D0) c ----- Pst_{NH3}(T=TmaxNH3) ----- TMXNH3=DBLE(APPLIM(5))/CONST(1) PS3R=FMF038(CONST,COEFF,TMXNH3) c ----- Pst_{H2O}(T=TmaxH2O) ----- TMXH2O=DBLE(APPLIM(7))/CONST(1) PS4R=FMF037(CONST,COEFF,TMXH2O) c write(*,*) ' ps4r=',ps4r,pressr c -------------------------------- IF (PRESSR.LT.PS1R) THEN J=-2 CALL FMF049(J,PRNAME) c write(*,*) ' 1' FMF045=-1.0D20 RETURN ELSEIF ((PRESSR.GE.PS1R).AND.(PRESSR.LT.PS2R)) THEN TSAT1R=2.3D0 CALL FMF050(J,CONST,COEFF,TSAT1R,PRESSR,MOLX1P,MOLY1P) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) c write(*,*) ' 2' FMF045=-1.0D10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) c write(*,*) ' 3' FMF045=-1.0D20 RETURN ENDIF MOLY1=MOLY1P-MOLFR MOLY0=0.0D0-MOLFR FLMOL=MOLY0*MOLY1 IF (FLMOL.GT.0.0D0) THEN J=-2 c write(*,*) ' 4' CALL FMF049(J,PRNAME) FMF045=-1.0D20 RETURN ELSEIF (FLMOL.EQ.0.0D0) THEN IF (MOLY0.EQ.0.0D0) THEN FMF045=TSAT0R RETURN ELSEIF (MOLY1.EQ.0.0D0) THEN FMF045=TSAT1R RETURN ENDIF ENDIF TSAT0R=FMF041(CONST,COEFF,PRESSR) IF (TSAT0R.LT.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) c write(*,*) ' 5' FMF045=-1.0D20 RETURN ELSEIF (TSAT0R.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) c write(*,*) ' 6' FMF045=-1.0D10 RETURN ENDIF ELSEIF ((PRESSR.GE.PS2R).AND.(PRESSR.LE.PS3R)) THEN MOLY0=0.0D0-MOLFR MOLY1=1.0D0-MOLFR TSAT0R=FMF041(CONST,COEFF,PRESSR) IF (TSAT0R.LT.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) c write(*,*) ' 7' FMF045=-1.0D20 RETURN ELSEIF (TSAT0R.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) c write(*,*) ' 8' FMF045=-1.0D10 RETURN ENDIF TSAT1R=FMF042(CONST,COEFF,PRESSR) IF (TSAT1R.LT.-1.0D15) THEN J=-2 c write(*,*) ' 9' CALL FMF049(J,PRNAME) FMF045=-1.0D20 RETURN ELSEIF (TSAT1R.LT.0.0D0) THEN J=-1 CALL FMF049(J,PRNAME) c write(*,*) ' 10' FMF045=-1.0D10 RETURN ENDIF ELSEIF (PRESSR.GT.PS3R.AND.PRESSR.LE.PS4R) THEN TSAT1R=(TMXH2O-TMXNH3)/(PS4R-PS3R)*(PRESSR-PS3R)+TMXNH3 CALL FMF050(J,CONST,COEFF,TSAT1R,PRESSR,MOLX1P,MOLY1P) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FMF045=-1.0D10 c write(*,*) ' 11' RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FMF045=-1.0D20 c write(*,*) ' 12' RETURN ENDIF MOLY1=MOLY1P-MOLFR MOLY0=0.0D0-MOLFR TSAT0R=FMF041(CONST,COEFF,PRESSR) IF (TSAT0R.LT.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) c write(*,*) ' 13' FMF045=-1.0D20 RETURN ELSEIF (TSAT0R.LT.0.0D0) THEN J=-1 c write(*,*) ' 14' CALL FMF049(J,PRNAME) FMF045=-1.0D10 RETURN ENDIF FLMOL=MOLY0*MOLY1 IF (FLMOL.GT.0.0D0) THEN J=-2 CALL FMF049(J,PRNAME) c write(*,*) ' 15' FMF045=-1.0D20 RETURN ELSEIF (FLMOL.EQ.0.0D0) THEN IF (MOLY0.EQ.0.0D0) THEN FMF045=TSAT0R RETURN ELSEIF (MOLY1.EQ.0.0D0) THEN FMF045=TSAT1R RETURN ENDIF ENDIF ELSEIF (PRESSR.GT.PS4R) THEN J=-2 c write(*,*) ' 16' CALL FMF049(J,PRNAME) RETURN ENDIF DO 10 I=1,100 TSAT2R=(TSAT0R+TSAT1R)*0.5D0 CALL FMF050(J,CONST,COEFF,TSAT2R,PRESSR,MOLX2P,MOLY2P) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) FMF045=-1.0D10 c write(*,*) ' 17' RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) FMF045=-1.0D20 c write(*,*) ' 18' RETURN ENDIF MOLY2=MOLY2P-MOLFR FLMOL=MOLY0*MOLY2 IF (FLMOL.GT.0.0D0) THEN MOLY0=MOLY2 TSAT0R=TSAT2R ELSE MOLY1=MOLY2 TSAT1R=TSAT2R ENDIF TSAT2R=TSAT1R-MOLY1*(TSAT1R-TSAT0R)/(MOLY1-MOLY0) CALL FMF050(J,CONST,COEFF,TSAT2R,PRESSR,MOLX2P,MOLY2P) IF (J.EQ.-1) THEN CALL FMF049(J,PRNAME) c write(*,*) ' 19' FMF045=-1.0D10 RETURN ELSEIF (J.EQ.-2) THEN CALL FMF049(J,PRNAME) c write(*,*) ' 20' FMF045=-1.0D20 RETURN ENDIF MOLY2=MOLY2P-MOLFR FLMOL=MOLY0*MOLY2 IF (FLMOL.GT.0.0D0) THEN MOLY0=MOLY2 TSAT0R=TSAT2R ELSE MOLY1=MOLY2 TSAT1R=TSAT2R ENDIF TSATE=TSAT1R**2 IF (TSATE.LE.1.0D0) THEN EPTSAT=(TSAT1R-TSAT0R)**2 EPMOL=MOLY2**2 IF ((EPTSAT.LE.EP1).AND.(EPMOL.LE.EP2)) THEN FMF045=TSAT2R RETURN ENDIF ELSE EPTSAT=((TSAT1R-TSAT0R)/TSAT1R)**2 EPMOL=MOLY2**2 IF ((EPTSAT.LE.EP1).AND.(EPMOL.LE.EP2)) THEN FMF045=TSAT2R RETURN ENDIF ENDIF 10 CONTINUE J=-1 c write(*,*) ' 21' CALL FMF049(J,PRNAME) FMF045=-1.0D10 RETURN END c ----------------------------------------------------------------- c Mole fraction,Pressure_R,Enthalpy_R => Temperature_R c ----------------------------------------------------------------- SUBROUTINE FMF046(J,CONST,COEFF,MOLFR,PRESSR,ENTHLR,TEMPR,PHASE) INTEGER I,J,PHASE DOUBLE PRECISION CONST(1:30),COEFF(1:56) $ ,MOLFR,PRESSR,ENTHLR,TEMPR $ ,EP1,EP2,TBUBR,TDEWR,HBUBR,HDEWR,TEMPR0,TEMPR1,TEMPR2 $ ,ENTHR0,ENTHR1,ENTHR2,ENTFL,TEMFL,QUAL2,DBX2,DBY2 $ ,FLTMPR,FLENTR $ ,FMF025,FMF026,FMF044,FMF045 CHARACTER*6 PRNAME C PRNAME='FMF046' J=0 EP1=1.0D-14 EP2=1.0D-14 TBUBR=FMF044(CONST,COEFF,PRESSR,MOLFR) IF (TBUBR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (TBUBR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF TDEWR=FMF045(CONST,COEFF,PRESSR,MOLFR) IF (TDEWR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (TDEWR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF HBUBR=FMF025(COEFF,TBUBR,PRESSR,MOLFR) HDEWR=FMF026(COEFF,TDEWR,PRESSR,MOLFR) IF (ENTHLR.LT.HBUBR) THEN PHASE=1 TEMPR0=TBUBR TEMPR1=TBUBR-1.0D-7 ENTHR0=HBUBR-ENTHLR ENTHR1=FMF025(COEFF,TEMPR1,PRESSR,MOLFR)-ENTHLR DO 10 I=1,100 TEMPR2=TEMPR1-ENTHR1*(TEMPR1-TEMPR0)/(ENTHR1-ENTHR0) ENTHR2=FMF025(COEFF,TEMPR2,PRESSR,MOLFR)-ENTHLR TEMFL=((TEMPR2-TEMPR1)/TEMPR1)**2 IF (TEMFL.LT.EP1) THEN ENTFL=ENTHR2**2 IF (ENTFL.LT.EP2) THEN TEMPR=TEMPR2 RETURN ENDIF ENDIF TEMPR0=TEMPR1 ENTHR0=ENTHR1 TEMPR1=TEMPR2 ENTHR1=ENTHR2 10 CONTINUE J=-1 CALL FMF049(J,PRNAME) RETURN ELSEIF (ENTHLR.EQ.HBUBR) THEN PHASE=1 TEMPR=TBUBR RETURN ELSEIF ((ENTHLR.GT.HBUBR).AND.(ENTHLR.LT.HDEWR)) THEN PHASE=2 TEMPR0=TBUBR TEMPR1=TDEWR ENTHR0=FMF025(COEFF,TEMPR0,PRESSR,MOLFR)-ENTHLR ENTHR1=FMF026(COEFF,TEMPR1,PRESSR,MOLFR)-ENTHLR ENTFL=ENTHR0*ENTHR1 IF (ENTFL.GT.0.0D0) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (ENTFL.EQ.0.0D0) THEN IF (ENTHR0.EQ.0.0D0) THEN TEMPR=TEMPR0 RETURN ELSEIF (ENTHR1.EQ.0.0D0) THEN TEMPR=TEMPR1 RETURN ENDIF ENDIF DO 20 I=1,100 TEMPR2=(TEMPR1+TEMPR0)*0.5D0 CALL FMF050(J,CONST,COEFF,TEMPR2,PRESSR,DBX2,DBY2) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF QUAL2=(MOLFR-DBX2)/(DBY2-DBX2) ENTHR2=(1.0D0-QUAL2)*FMF025(COEFF,TEMPR2,PRESSR,DBX2) $ +QUAL2*FMF026(COEFF,TEMPR2,PRESSR,DBY2)-ENTHLR ENTFL=ENTHR0*ENTHR2 IF (ENTFL.GT.0.0D0) THEN TEMPR0=TEMPR2 ENTHR0=ENTHR2 ELSEIF (ENTFL.EQ.0.0D0) THEN IF (ENTHR2.EQ.0.0D0) THEN TEMPR=TEMPR2 RETURN ENDIF ELSE TEMPR1=TEMPR2 ENTHR1=ENTHR2 ENDIF TEMPR2=TEMPR1-ENTHR1*(TEMPR1-TEMPR0)/(ENTHR1-ENTHR0) CALL FMF050(J,CONST,COEFF,TEMPR2,PRESSR,DBX2,DBY2) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF QUAL2=(MOLFR-DBX2)/(DBY2-DBX2) ENTHR2=(1.0D0-QUAL2)*FMF025(COEFF,TEMPR2,PRESSR,DBX2) $ +QUAL2*FMF026(COEFF,TEMPR2,PRESSR,DBY2)-ENTHLR ENTFL=ENTHR0*ENTHR2 IF (ENTFL.GT.0.0D0) THEN TEMPR0=TEMPR2 ENTHR0=ENTHR2 ELSEIF (ENTFL.EQ.0.0D0) THEN IF (ENTHR2.EQ.0.0D0) THEN TEMPR=TEMPR2 RETURN ENDIF ELSE TEMPR1=TEMPR2 ENTHR1=ENTHR2 ENDIF FLENTR=ENTHR2**2 FLTMPR=((TEMPR0-TEMPR1)/TEMPR1)**2 IF (FLENTR.LT.EP1.AND.FLTMPR.LT.EP2) THEN J=0 TEMPR=TEMPR2 RETURN ENDIF 20 CONTINUE J=-1 CALL FMF049(J,PRNAME) RETURN ELSEIF (ENTHLR.EQ.HDEWR) THEN PHASE=3 TEMPR=TDEWR RETURN ELSEIF (ENTHLR.GT.HDEWR) THEN PHASE=3 TEMPR0=TDEWR TEMPR1=TDEWR+1.0D-7 ENTHR0=HDEWR-ENTHLR ENTHR1=FMF026(COEFF,TEMPR1,PRESSR,MOLFR)-ENTHLR DO 30 I=1,100 TEMPR2=TEMPR1-ENTHR1*(TEMPR1-TEMPR0)/(ENTHR1-ENTHR0) ENTHR2=FMF026(COEFF,TEMPR2,PRESSR,MOLFR)-ENTHLR TEMFL=((TEMPR2-TEMPR1)/TEMPR1)**2 IF (TEMFL.LT.EP1) THEN ENTFL=ENTHR2**2 IF (ENTFL.LT.EP2) THEN TEMPR=TEMPR2 RETURN ENDIF ENDIF TEMPR0=TEMPR1 ENTHR0=ENTHR1 TEMPR1=TEMPR2 ENTHR1=ENTHR2 30 CONTINUE J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF J=-2 CALL FMF049(J,PRNAME) RETURN END c ----------------------------------------------------------------- c Mole fraction,Pressure_R,Entropy_R => Temperature_R c ----------------------------------------------------------------- SUBROUTINE FMF047(J,CONST,COEFF,MOLFR,PRESSR,ENTROR,TEMPR,PHASE) INTEGER I,J,PHASE DOUBLE PRECISION CONST(1:30),COEFF(1:56) $ ,MOLFR,PRESSR,ENTROR,TEMPR $ ,EP1,EP2,TBUBR,TDEWR,SBUBR,SDEWR,TEMPR0,TEMPR1,TEMPR2 $ ,ENTRR0,ENTRR1,ENTRR2,ENTFL,TEMFL,QUAL2,DBX2,DBY2,FLENTR $ ,FLTMPR $ ,FMF027,FMF028,FMF044,FMF045 CHARACTER*6 PRNAME C PRNAME='FMF047' J=0 EP1=1.0D-14 EP2=1.0D-14 TBUBR=FMF044(CONST,COEFF,PRESSR,MOLFR) IF (TBUBR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (TBUBR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF TDEWR=FMF045(CONST,COEFF,PRESSR,MOLFR) IF (TDEWR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (TDEWR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF SBUBR=FMF027(COEFF,TBUBR,PRESSR,MOLFR) SDEWR=FMF028(COEFF,TDEWR,PRESSR,MOLFR) IF (ENTROR.LT.SBUBR) THEN PHASE=1 TEMPR0=TBUBR TEMPR1=TBUBR-1.0D-7 ENTRR0=SBUBR-ENTROR ENTRR1=FMF027(COEFF,TEMPR1,PRESSR,MOLFR)-ENTROR DO 10 I=1,100 TEMPR2=TEMPR1-ENTRR1*(TEMPR1-TEMPR0)/(ENTRR1-ENTRR0) ENTRR2=FMF027(COEFF,TEMPR2,PRESSR,MOLFR)-ENTROR TEMFL=((TEMPR2-TEMPR1)/TEMPR1)**2 IF (TEMFL.LT.EP1) THEN ENTFL=ENTRR2**2 IF (ENTFL.LT.EP2) THEN TEMPR=TEMPR2 RETURN ENDIF ENDIF TEMPR0=TEMPR1 ENTRR0=ENTRR1 TEMPR1=TEMPR2 ENTRR1=ENTRR2 10 CONTINUE ELSEIF (ENTROR.EQ.SBUBR) THEN PHASE=1 TEMPR=TBUBR RETURN ELSEIF ((ENTROR.GT.SBUBR).AND.(ENTROR.LT.SDEWR)) THEN PHASE=2 TEMPR0=TBUBR TEMPR1=TDEWR ENTRR0=FMF027(COEFF,TEMPR0,PRESSR,MOLFR)-ENTROR ENTRR1=FMF028(COEFF,TEMPR1,PRESSR,MOLFR)-ENTROR ENTFL=ENTRR0*ENTRR1 IF (ENTFL.GT.0.0D0) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (ENTFL.EQ.0.0D0) THEN IF (ENTRR0.EQ.0.0D0) THEN TEMPR=TEMPR0 RETURN ELSEIF (ENTRR1.EQ.0.0D0) THEN TEMPR=TEMPR1 RETURN ENDIF ENDIF DO 20 I=1,100 TEMPR2=(TEMPR1+TEMPR0)*0.5D0 CALL FMF050(J,CONST,COEFF,TEMPR2,PRESSR,DBX2,DBY2) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF QUAL2=(MOLFR-DBX2)/(DBY2-DBX2) ENTRR2=(1.0D0-QUAL2)*FMF027(COEFF,TEMPR2,PRESSR,DBX2) $ +QUAL2*FMF028(COEFF,TEMPR2,PRESSR,DBY2)-ENTROR ENTFL=ENTRR0*ENTRR2 IF (ENTFL.GT.0.0D0) THEN TEMPR0=TEMPR2 ENTRR0=ENTRR2 ELSEIF (ENTFL.EQ.0.0D0) THEN IF (ENTRR2.EQ.0.0D0) THEN TEMPR=TEMPR2 RETURN ENDIF ELSE TEMPR1=TEMPR2 ENTRR1=ENTRR2 ENDIF TEMPR2=TEMPR1-ENTRR1*(TEMPR1-TEMPR0)/(ENTRR1-ENTRR0) CALL FMF050(J,CONST,COEFF,TEMPR2,PRESSR,DBX2,DBY2) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF QUAL2=(MOLFR-DBX2)/(DBY2-DBX2) ENTRR2=(1.0D0-QUAL2)*FMF027(COEFF,TEMPR2,PRESSR,DBX2) $ +QUAL2*FMF028(COEFF,TEMPR2,PRESSR,DBY2)-ENTROR ENTFL=ENTRR0*ENTRR2 IF (ENTFL.GT.0.0D0) THEN TEMPR0=TEMPR2 ENTRR0=ENTRR2 ELSEIF (ENTFL.EQ.0.0D0) THEN IF (ENTRR2.EQ.0.0D0) THEN TEMPR=TEMPR2 RETURN ENDIF ELSE TEMPR1=TEMPR2 ENTRR1=ENTRR2 ENDIF FLENTR=ENTRR2**2 FLTMPR=((TEMPR0-TEMPR1)/TEMPR1)**2 IF (FLENTR.LT.EP1.AND.FLTMPR.LT.EP2) THEN J=0 TEMPR=TEMPR2 RETURN ENDIF 20 CONTINUE ELSEIF (ENTROR.EQ.SDEWR) THEN PHASE=3 TEMPR=TDEWR RETURN ELSEIF (ENTROR.GT.SDEWR) THEN PHASE=3 TEMPR0=TDEWR TEMPR1=TDEWR+1.0D-7 ENTRR0=SDEWR-ENTROR ENTRR1=FMF028(COEFF,TEMPR1,PRESSR,MOLFR)-ENTROR DO 30 I=1,100 TEMPR2=TEMPR1-ENTRR1*(TEMPR1-TEMPR0)/(ENTRR1-ENTRR0) ENTRR2=FMF028(COEFF,TEMPR2,PRESSR,MOLFR)-ENTROR TEMFL=((TEMPR2-TEMPR1)/TEMPR1)**2 IF (TEMFL.LT.EP1) THEN ENTFL=ENTRR2**2 IF (ENTFL.LT.EP2) THEN TEMPR=TEMPR2 RETURN ENDIF ENDIF TEMPR0=TEMPR1 ENTRR0=ENTRR1 TEMPR1=TEMPR2 ENTRR1=ENTRR2 30 CONTINUE ENDIF J=-1 CALL FMF049(J,PRNAME) RETURN END c ----------------------------------------------------------------- c Mole fraction,Pressure_R,Volume_R => Temperature_R c ----------------------------------------------------------------- SUBROUTINE FMF048(J,CONST,COEFF,MOLFR,PRESSR,VOLUMR,TEMPR,PHASE) INTEGER I,J,PHASE DOUBLE PRECISION CONST(1:30),COEFF(1:56) $ ,MOLFR,PRESSR,VOLUMR,TEMPR $ ,EP1,EP2,TBUBR,TDEWR,VBUBR,VDEWR,TEMPR0,TEMPR1,TEMPR2 $ ,VOLMR0,VOLMR1,VOLMR2,VOLFL,TEMFL,QUAL2,DBX2,DBY2 $ ,FLTMPR,FLVOLR $ ,FMF023,FMF024,FMF044,FMF045 CHARACTER*6 PRNAME C PRNAME='FMF048' J=0 EP1=1.0D-14 EP2=1.0D-14 TBUBR=FMF044(CONST,COEFF,PRESSR,MOLFR) IF (TBUBR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (TBUBR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF TDEWR=FMF045(CONST,COEFF,PRESSR,MOLFR) IF (TDEWR.LE.-1.0D15) THEN J=-2 CALL FMF049(J,PRNAME) RETURN ELSEIF (TDEWR.LE.-9.99D9) THEN J=-1 CALL FMF049(J,PRNAME) RETURN ENDIF VBUBR=FMF023(COEFF,TBUBR,PRESSR,MOLFR) VDEWR=FMF024(COEFF,TDEWR,PRESSR,MOLFR) IF (VOLUMR.LT.VBUBR) THEN PHASE=1 TEMPR0=TBUBR TEMPR1=TBUBR-1.0D-7 VOLMR0=VBUBR-VOLUMR VOLMR1=FMF023(COEFF,TEMPR1,PRESSR,MOLFR)-VOLUMR DO 10 I=1,100 TEMPR2=TEMPR1-VOLMR1*(TEMPR1-TEMPR0)/(VOLMR1-VOLMR0) VOLMR2=FMF023(COEFF,TEMPR2,PRESSR,MOLFR)-VOLUMR TEMFL=((TEMPR2-TEMPR1)/TEMPR1)**2 IF (TEMFL.LT.EP1) THEN VOLFL=VOLMR2**2 IF (VOLFL.LT.EP2) THEN TEMPR=TEMPR2 RETURN ENDIF ENDIF TEMPR0=TEMPR1 VOLMR0=VOLMR1 TEMPR1=TEMPR2 VOLMR1=VOLMR2 10 CONTINUE ELSEIF (VOLUMR.EQ.VBUBR) THEN PHASE=1 TEMPR=TBUBR RETURN ELSEIF ((VOLUMR.GT.VBUBR).AND.(VOLUMR.LT.VDEWR)) THEN PHASE=2 TEMPR0=TBUBR TEMPR1=TDEWR VOLMR0=FMF023(COEFF,TEMPR0,PRESSR,MOLFR)-VOLUMR VOLMR1=FMF024(COEFF,TEMPR1,PRESSR,MOLFR)-VOLUMR VOLFL=VOLMR0*VOLMR1 IF (VOLFL.GT.0.0D0) THEN J=-2 RETURN ELSEIF (VOLFL.EQ.0.0D0) THEN IF (VOLMR0.EQ.0.0D0) THEN TEMPR=TEMPR0 RETURN ELSEIF (VOLMR1.EQ.0.0D0) THEN TEMPR=TEMPR1 RETURN ENDIF ENDIF DO 20 I=1,100 TEMPR2=(TEMPR1+TEMPR0)*0.5D0 CALL FMF050(J,CONST,COEFF,TEMPR2,PRESSR,DBX2,DBY2) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF QUAL2=(MOLFR-DBX2)/(DBY2-DBX2) VOLMR2=(1.0D0-QUAL2)*FMF023(COEFF,TEMPR2,PRESSR,DBX2) $ +QUAL2*FMF024(COEFF,TEMPR2,PRESSR,DBY2)-VOLUMR VOLFL=VOLMR0*VOLMR2 IF (VOLFL.GT.0.0D0) THEN TEMPR0=TEMPR2 VOLMR0=VOLMR2 ELSEIF (VOLFL.EQ.0.0D0) THEN IF (VOLMR2.EQ.0.0D0) THEN TEMPR=TEMPR2 RETURN ENDIF ELSE TEMPR1=TEMPR2 VOLMR1=VOLMR2 ENDIF TEMPR2=TEMPR1-VOLMR1*(TEMPR1-TEMPR0)/(VOLMR1-VOLMR0) CALL FMF050(J,CONST,COEFF,TEMPR2,PRESSR,DBX2,DBY2) IF (J.NE.0) THEN CALL FMF049(J,PRNAME) RETURN ENDIF QUAL2=(MOLFR-DBX2)/(DBY2-DBX2) VOLMR2=(1.0D0-QUAL2)*FMF023(COEFF,TEMPR2,PRESSR,DBX2) $ +QUAL2*FMF024(COEFF,TEMPR2,PRESSR,DBY2)-VOLUMR VOLFL=VOLMR0*VOLMR2 IF (VOLFL.GT.0.0D0) THEN TEMPR0=TEMPR2 VOLMR0=VOLMR2 ELSEIF (VOLFL.EQ.0.0D0) THEN IF (VOLMR2.EQ.0.0D0) THEN TEMPR=TEMPR2 RETURN ENDIF ELSE TEMPR1=TEMPR2 VOLMR1=VOLMR2 ENDIF FLTMPR=((TEMPR0-TEMPR1)/TEMPR1)**2 FLVOLR=VOLMR2**2 IF (FLVOLR.LT.EP1.AND.FLTMPR.LT.EP2) THEN J=0 TEMPR=TEMPR2 RETURN ENDIF 20 CONTINUE ELSEIF (VOLUMR.EQ.VDEWR) THEN PHASE=3 TEMPR=TDEWR RETURN ELSEIF (VOLUMR.GT.VDEWR) THEN PHASE=3 TEMPR0=TDEWR TEMPR1=TDEWR+1.0D-7 VOLMR0=VDEWR-VOLUMR VOLMR1=FMF024(COEFF,TEMPR1,PRESSR,MOLFR)-VOLUMR DO 30 I=1,100 TEMPR2=TEMPR1-VOLMR1*(TEMPR1-TEMPR0)/(VOLMR1-VOLMR0) VOLMR2=FMF024(COEFF,TEMPR2,PRESSR,MOLFR)-VOLUMR TEMFL=((TEMPR2-TEMPR1)/TEMPR1)**2 IF (TEMFL.LT.EP1) THEN VOLFL=VOLMR2**2 IF (VOLFL.LT.EP2) THEN TEMPR=TEMPR2 RETURN ENDIF ENDIF TEMPR0=TEMPR1 VOLMR0=VOLMR1 TEMPR1=TEMPR2 VOLMR1=VOLMR2 30 CONTINUE ENDIF J=-1 CALL FMF049(J,PRNAME) RETURN END C ---------------------------------------- C Error message C ---------------------------------------- SUBROUTINE FMF049(J,PRNAME) INTEGER KPA,MESS,KSTAN,KAS,J COMMON /UNIT/ KPA,MESS,KSTAN,KAS CHARACTER*(*) PRNAME C IF (MESS.EQ.1) THEN IF (J.EQ.-1) THEN WRITE(*,1) PRNAME ELSEIF (J.EQ.-2) THEN WRITE(*,2) PRNAME ELSEIF (J.EQ.-3) THEN WRITE(*,3) PRNAME ELSEIF (J.EQ.-4) THEN WRITE(*,4) PRNAME ENDIF ENDIF 1 FORMAT (1X,'***** ERROR IN',1X,A6,1X,', NO CONVERGENCE *****') 2 FORMAT (1X,'***** ERROR IN',1X,A6,1X,', OUT OF RANGE *****') 3 FORMAT (1X,'***** ERROR IN',1X,A6,1X,', log DOMAIN ERROR *****') 4 FORMAT (1X,'***** ERROR IN',1X,A6,1X,', DEVIDE BY 0 *****') RETURN END c ------------------------------------------------------------------ c Temp_R,Pressure_R ===> x,y c ------------------------------------------------------------------ SUBROUTINE FMF050(J,CONST,COEFF,TEMPR,PRESSR,MOLX,MOLY) INTEGER I,J DOUBLE PRECISION COEFF(1:56),TEMPR,PRESSR,MOLX,MOLY,MOLX1,MOLY1 $ ,MOLX2,MOLY2,K10,K20,K1,K2,CHMPY1,CHMPY2,CHMPX1,CHMPX2 $ ,KY1,KY2,KX1,KX2,EP1,EP2,KE1,KE2,CONST(1:30) $ ,FMF029,FMF030,FMF031,FMF032 CHARACTER*6 PRNAME C PRNAME='FMF050' J=0 EP1=1.0D-14 EP2=1.0D-14 MOLX1=1.0D-10 MOLX2=1.0D0-MOLX1 MOLY1=1.0D0-1.0D-10 MOLY2=1.0D0-MOLY1 K10=MOLY1/MOLX1 K20=MOLY2/MOLX2 DO 10 I=1,1000 CHMPY1=DEXP(FMF032(COEFF,TEMPR,PRESSR,MOLY1)/(CONST(3)*TEMPR)) CHMPY2=DEXP(FMF031(COEFF,TEMPR,PRESSR,MOLY1)/(CONST(3)*TEMPR)) KY1=CHMPY1/MOLY1 KY2=CHMPY2/MOLY2 CHMPX1=DEXP(FMF030(COEFF,TEMPR,PRESSR,MOLX1)/(CONST(3)*TEMPR)) CHMPX2=DEXP(FMF029(COEFF,TEMPR,PRESSR,MOLX1)/(CONST(3)*TEMPR)) KX1=CHMPX1/MOLX1 KX2=CHMPX2/MOLX2 K1=KX1/KY1 K2=KX2/KY2 KE1=((K10-K1)/K1)**2 KE2=((K20-K2)/K2)**2 c WRITE(*,*) ' K1,K2=',ke1,ke2 IF (KE1.LT.EP1.AND.KE2.LT.EP2) THEN J=0 MOLX=MOLX1 MOLY=MOLY1 RETURN ENDIF MOLX1=(K2-1.0D0)/(K2-K1) MOLY1=K1*MOLX1 MOLX2=1.0D0-MOLX1 MOLY2=1.0D0-MOLY1 K10=K1 K20=K2 10 CONTINUE J=-1 CALL FMF049(J,PRNAME) RETURN END c ----------------------------------------------- c This program transfers the unit of fraction c from kmol to kg. c ----------------------------------------------- DOUBLE PRECISION FUNCTION FMF051(CONST,DMOL) DOUBLE PRECISION CONST(1:30),DMOL c FMF051=DMOL*CONST(11)/(DMOL*CONST(11) $ +(1.0D0-DMOL)*CONST(21)) RETURN END c ----------------------------------------------- c This program transfers the unit of fraction c from kg to kmol. c ----------------------------------------------- DOUBLE PRECISION FUNCTION FMF052(CONST,DMASS) DOUBLE PRECISION CONST(1:30),DMASS c FMF052=DMASS*CONST(21)/(DMASS*CONST(21) $ +(1.0D0-DMASS)*CONST(11)) RETURN END c c ------ SUBROUTINE SUBMXH and SUBMXT were added by R.Akasaka (MAY 14, 1998) ----- C ---------------------------------------------------------- C ** SUBROUTINE SUBMXH ** C ** Thermal Properties of Mixtures of Mixed by ** C ** Two Different Compositions ** C ** [Input] ** C ** PP : Pressure [Pa],[bar] ** C ** ZA : Total Composition of A ** C ** ZB : Total Composition of B ** C ** HA : Enthalpy of A [J/kmol],[J/kg] ** C ** HB : Enthalpy of B [J/kmol],[J/kg] ** C ** W : Fraction of A in Mixing ** C ** [Output] ** C ** J : Error Detection Code ** C ** TT : Temperature after Mixing ** C ** Z : Total Composition after Mixing ** C ** V : Volume after Mixing ** C ** H : Enthalpy after Mixing ** C ** S : Entropy after Mixing ** C ---------------------------------------------------------- SUBROUTINE SUBMXH(J,T,P,ZA,ZB,Z,V,HA,HB,H,S,W) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,ZA,ZB,Z,V,HA,HB,H,S,W CHARACTER*6 PRNAME COMMON /UNIT/KPA,MESS,KSTAN,KAS PRNAME='SUBMXH' J=0 Z=W*ZA+(1-W)*ZB H=W*HA+(1-W)*HB CALL SUBMIX(2,J,T,P,Z,V,HH,S) RETURN END C ---------------------------------------------------------- C ** SUBROUTINE SUBMXT ** C ** Thermal Properties of Mixtures of Mixed by ** C ** Two Different Compositions ** C ** [Input] ** C ** PP : Pressure [Pa],[bar] ** C ** ZA : Total Composition of A ** C ** ZB : Total Composition of B ** C ** TA : Temperature of A [C],[K] ** C ** TB : Temperature of B [C],[K] ** C ** W : Fraction of A in Mixing ** C ** [Output] + ** C ** J : Error Detection Code ** C ** TT : Temperature after Mixing ** C ** Z : Total Composition after Mixing ** C ** V : Volume after Mixing ** C ** H : Enthalpy after Mixing ** C ** S : Entropy after Mixing ** C ---------------------------------------------------------- SUBROUTINE SUBMXT(J,T,P,ZA,ZB,Z,V,TA,TB,H,S,W) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,P,ZA,ZB,Z,V,TA,TB,H,S,W,HA,HB,SA,SB CHARACTER*6 PRNAME COMMON /UNIT/KPA,MESS,KSTAN,KAS PRNAME='SUBMXT' J=0 Z=W*ZA+(1-W)*ZB CALL SUBMIX(1,J,TA,P,ZA,VA,HA,SA) CALL SUBMIX(1,J,TB,P,ZB,VB,HB,SB) H=W*HA+(1-W)*HB CALL SUBMIX(2,J,T,P,Z,V,H,S) RETURN END C C