C *************************************************************** C *** *** C *** The Thermal Properties of Mixtures of Water and Ammonia *** C *** Helmholtz Function Provided by *** C *** Reiner Tillner-Roth and Daniel G.Friend *** C *** *** C *************************************************************** C ========================================================= c ** Helmholtz Free Energy of Ideal-Water Properties ** C ========================================================= SUBROUTINE H_IWHG(D,T,FA) DOUBLE PRECISION T,D,FA DOUBLE PRECISION No(1:8),Gamao(4:8) DATA (No(I),I=1,8) /-8.32044648201,6.6832105268,3.00632, $ 0.012436,0.97315,1.27950,0.96956,0.24873/ DATA (Gamao(I),I=4,8) /1.28728967,3.53734222,7.74073708, $ 9.24437796,27.5075105/ FA=DLOG(D)+No(1)+No(2)*T+No(3)*DLOG(T) $ +No(4)*DLOG(1.-DEXP(-Gamao(4)*T)) $ +No(5)*DLOG(1.-DEXP(-Gamao(5)*T)) $ +No(6)*DLOG(1.-DEXP(-Gamao(6)*T)) $ +No(7)*DLOG(1.-DEXP(-Gamao(7)*T)) $ +No(8)*DLOG(1.-DEXP(-Gamao(8)*T)) RETURN END C ==================================================== C ** The Helmholtz Free Energy of The residual part ** c ** OF WATER ** C ==================================================== SUBROUTINE H_RWHG(Dx,Tx,DRWHG) DOUBLE PRECISION Tx,Dx,DRWHG REAL C(8:51),D(1:54),t(1:54) REAL A(55:56),B(55:56),BB(55:56) REAL Alfa(52:54),Belt(52:56),Gama(52:54),El(52:54) REAL BC(55:56),BD(55:56),BA(55:56) DOUBLE PRECISION CN(1:56) DOUBLE PRECISION F,Delta,Sita,Fal DATA (C(I),I=8,51) /1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2, $ 2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,4,6,6,6,6/ DATA (D(I),I=1,54) /1,1,1,2,2,3,4,1,1,1,2,2,3,4,4,5,7,9,10,11, $ 13,15,1,2,2,2,3,4,4,4,5,6,6,7,9,9,9,9,9,10,10,12,3,4,4,5, $ 14,3,6,6,6,3,3,3/ DATA (t(I),I=1,54) /-0.5,0.875,1,0.5,0.75,0.375,1,4,6,12,1,5,4, $ 2,13,9,3,4,11,4,13,1,7,1,9,10,10,3,7,10,10,6,10,10,1,2,3,4, $ 8,6,9,8,16,22,23,23,10,50,44,46,50,0,1,4/ DATA (CN(I),I=1,56) /0.12533547935523D-1,0.78957634722828D+1, $ -0.87803203303561D+1,0.31802509345418,-0.26145533859358, $ -0.78199751687981D-2,0.88089493102134D-2,-0.66856572307965, $ 0.20433810950965,-0.66212605039687D-4,-0.19232721156002, $ -0.25709043003438,0.16074868486251,-0.40092828925807D-1, $ 0.39343422603254D-6,-0.75941377088144D-5,0.56250979351888D-3, $ -0.15608652257135D-4,0.11537996422951D-8,0.36582165144204D-6, $ -0.13251180074668D-11,-0.62639586912454D-9,-0.10793600908932, $ 0.17611491008752D-1,0.22132295167546,-0.40247669763528, $ 0.58083399985759,0.49969146990806D-2,-0.31358700712549D-1, $ -0.74315929710341,0.47807329915480,0.20527940895948D-1, $ -0.13636435110343,0.14180634400617D-1,0.83326504880713D-2, $ -0.29052336009585D-1,0.38615085574206D-1,-0.20393486513704D-1, $ -0.16554050063734D-2,0.19955571979541D-2,0.15870308324157D-3, $ -0.16388568342530D-4,0.43613615723811D-1,0.34994005463765D-1, $ -0.76788197844621D-1,0.22446277332006D-1,-0.62689710414685D-4, $ -0.55711118565645D-9,-0.19905718354408,0.31777497330738, $ -0.11841182425981,-0.31306260323435D2,0.31546140237781D2, $ -0.25213154341695D4,-0.14874640856724,0.31806110878444 / DATA (A(I),I=55,56) /3.5,3.5/ DATA (B(I),I=55,56) /0.85,0.95/ DATA (BB(I),I=55,56) /0.2,0.2/ DATA (Alfa(I),I=52,54) /20,20,20/ DATA (Belt(I),I=52,56) /150,150,250,0.3,0.3/ DATA (Gama(I),I=52,54) /1.21,1.21,1.25/ DATA (El(I),I=52,54) /1,1,1/ DATA (BC(I),I=55,56) /28,32/ DATA (BD(I),I=55,56) /700,800/ DATA (BA(I),I=55,56) /0.32,0.32/ DRWHG=0.0 DO 100 I=1,7 DRWHG=DRWHG+CN(I)*Dx**D(I)*Tx**t(I) 100 CONTINUE DO 200 I=8,51 DRWHG=DRWHG+CN(I)*Dx**D(I)*Tx**t(I)*DEXP(-Dx**C(I)) 200 CONTINUE DO 300 I=52,54 F=-Alfa(I)*(Dx-El(I))**2-Belt(I)*(Tx-Gama(I))**2 DRWHG=DRWHG+CN(I)*Dx**D(I)*Tx**t(I)*DEXP(F) 300 CONTINUE DO 400 I=55,56 Fal=DEXP(-BC(I)*(Dx-1)**2-BD(I)*(Tx-1)**2) Sita=(1-Tx)+BA(I)*((Dx-1)**2)**(0.5/Belt(I)) Delta=Sita*Sita+BB(I)*((Dx-1)**2)**A(I) DRWHG=DRWHG+CN(I)*Delta**B(I)*Dx*Fal 400 CONTINUE RETURN END C =============================================================== C ** Differential of the Helmholtz Free Energy of ** C ** the Residual Part to Dimensionless Temperature T OF WATER ** C =============================================================== SUBROUTINE D_R_WHT(Dx,Tx,DRWHT) DOUBLE PRECISION Tx,Dx,DRWHT REAL C(8:51),D(1:54),t(1:54) REAL A(55:56),B(55:56),BB(55:56) REAL Alfa(52:54),Belt(52:56),Gama(52:54),El(52:54) REAL BC(55:56),BD(55:56),BA(55:56) DOUBLE PRECISION CN(1:56) DOUBLE PRECISION F,Delta,Sita,Fal,G,Deltat,Falt DATA (C(I),I=8,51) /1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2, $ 2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,4,6,6,6,6/ DATA (D(I),I=1,54) /1,1,1,2,2,3,4,1,1,1,2,2,3,4,4,5,7,9,10,11, $ 13,15,1,2,2,2,3,4,4,4,5,6,6,7,9,9,9,9,9,10,10,12,3,4,4,5, $ 14,3,6,6,6,3,3,3/ DATA (t(I),I=1,54) /-0.5,0.875,1,0.5,0.75,0.375,1,4,6,12,1,5,4, $ 2,13,9,3,4,11,4,13,1,7,1,9,10,10,3,7,10,10,6,10,10,1,2,3,4, $ 8,6,9,8,16,22,23,23,10,50,44,46,50,0,1,4/ DATA (CN(I),I=1,56) /0.12533547935523D-1,0.78957634722828D+1, $ -0.87803203303561D+1,0.31802509345418,-0.26145533859358, $ -0.78199751687981D-2,0.88089493102134D-2,-0.66856572307965, $ 0.20433810950965,-0.66212605039687D-4,-0.19232721156002, $ -0.25709043003438,0.16074868486251,-0.40092828925807D-1, $ 0.39343422603254D-6,-0.75941377088144D-5,0.56250979351888D-3, $ -0.15608652257135D-4,0.11537996422951D-8,0.36582165144204D-6, $ -0.13251180074668D-11,-0.62639586912454D-9,-0.10793600908932, $ 0.17611491008752D-1,0.22132295167546,-0.40247669763528, $ 0.58083399985759,0.49969146990806D-2,-0.31358700712549D-1, $ -0.74315929710341,0.47807329915480,0.20527940895948D-1, $ -0.13636435110343,0.14180634400617D-1,0.83326504880713D-2, $ -0.29052336009585D-1,0.38615085574206D-1,-0.20393486513704D-1, $ -0.16554050063734D-2,0.19955571979541D-2,0.15870308324157D-3, $ -0.16388568342530D-4,0.43613615723811D-1,0.34994005463765D-1, $ -0.76788197844621D-1,0.22446277332006D-1,-0.62689710414685D-4, $ -0.55711118565645D-9,-0.19905718354408,0.31777497330738, $ -0.11841182425981,-0.31306260323435D2,0.31546140237781D2, $ -0.25213154341695D4,-0.14874640856724,0.31806110878444 / DATA (A(I),I=55,56) /3.5,3.5/ DATA (B(I),I=55,56) /0.85,0.95/ DATA (BB(I),I=55,56) /0.2,0.2/ DATA (Alfa(I),I=52,54) /20,20,20/ DATA (Belt(I),I=52,56) /150,150,250,0.3,0.3/ DATA (Gama(I),I=52,54) /1.21,1.21,1.25/ DATA (El(I),I=52,54) /1,1,1/ DATA (BC(I),I=55,56) /28,32/ DATA (BD(I),I=55,56) /700,800/ DATA (BA(I),I=55,56) /0.32,0.32/ DRWHT=0 DO 100 I=1,7 DRWHT=DRWHT+CN(I)*t(I)*Dx**D(I)*Tx**(t(I)-1) 100 CONTINUE DO 200 I=8,51 DRWHT=DRWHT+CN(I)*t(I)*Dx**D(I)*Tx**(t(I)-1) $ *DEXP(-Dx**C(I)) 200 CONTINUE DO 300 I=52,54 F=-Alfa(I)*(Dx-El(I))**2-Belt(I)*(Tx-Gama(I))**2 G=t(I)/Tx-2*Belt(I)*(Tx-Gama(I)) DRWHT=DRWHT+CN(I)*Dx**D(I)*Tx**t(I)*DEXP(F)*G 300 CONTINUE DO 400 I=55,56 Fal=DEXP(-BC(I)*(Dx-1)**2-BD(I)*(Tx-1)**2) Sita=(1-Tx)+BA(I)*((Dx-1)**2)**(0.5/Belt(I)) Delta=Sita*Sita+BB(I)*((Dx-1)**2)**A(I) Deltat=-2.*Sita*B(I)*Delta**(B(I)-1) Falt=-2*BD(I)*(Tx-1)*Fal G=Deltat*Fal+Delta**B(I)*Falt DRWHT=DRWHT+CN(I)*Dx*G 400 CONTINUE RETURN END C ================================================================ C ** Differential of the Helmholtz Free Energy of ** C ** the Residual Part to Dimensionless Temperature TT OF WATER ** C ================================================================ SUBROUTINE D_R_WHTT(Dx,Tx,DRWHTT) DOUBLE PRECISION Tx,Dx,DRWHTT REAL C(8:51),D(1:54),t(1:54) REAL A(55:56),B(55:56),BB(55:56) REAL Alfa(52:54),Belt(52:56),Gama(52:54),El(52:54) REAL BC(55:56),BD(55:56),BA(55:56) DOUBLE PRECISION CN(1:56) DOUBLE PRECISION F,Delta,Sita,Fal,G,Deltat,Falt,Dltatt,Faltt DATA (C(I),I=8,51) /1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2, $ 2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,4,6,6,6,6/ DATA (D(I),I=1,54) /1,1,1,2,2,3,4,1,1,1,2,2,3,4,4,5,7,9,10,11, $ 13,15,1,2,2,2,3,4,4,4,5,6,6,7,9,9,9,9,9,10,10,12,3,4,4,5, $ 14,3,6,6,6,3,3,3/ DATA (t(I),I=1,54) /-0.5,0.875,1,0.5,0.75,0.375,1,4,6,12,1,5,4, $ 2,13,9,3,4,11,4,13,1,7,1,9,10,10,3,7,10,10,6,10,10,1,2,3,4, $ 8,6,9,8,16,22,23,23,10,50,44,46,50,0,1,4/ DATA (CN(I),I=1,56) /0.12533547935523D-1,0.78957634722828D+1, $ -0.87803203303561D+1,0.31802509345418,-0.26145533859358, $ -0.78199751687981D-2,0.88089493102134D-2,-0.66856572307965, $ 0.20433810950965,-0.66212605039687D-4,-0.19232721156002, $ -0.25709043003438,0.16074868486251,-0.40092828925807D-1, $ 0.39343422603254D-6,-0.75941377088144D-5,0.56250979351888D-3, $ -0.15608652257135D-4,0.11537996422951D-8,0.36582165144204D-6, $ -0.13251180074668D-11,-0.62639586912454D-9,-0.10793600908932, $ 0.17611491008752D-1,0.22132295167546,-0.40247669763528, $ 0.58083399985759,0.49969146990806D-2,-0.31358700712549D-1, $ -0.74315929710341,0.47807329915480,0.20527940895948D-1, $ -0.13636435110343,0.14180634400617D-1,0.83326504880713D-2, $ -0.29052336009585D-1,0.38615085574206D-1,-0.20393486513704D-1, $ -0.16554050063734D-2,0.19955571979541D-2,0.15870308324157D-3, $ -0.16388568342530D-4,0.43613615723811D-1,0.34994005463765D-1, $ -0.76788197844621D-1,0.22446277332006D-1,-0.62689710414685D-4, $ -0.55711118565645D-9,-0.19905718354408,0.31777497330738, $ -0.11841182425981,-0.31306260323435D2,0.31546140237781D2, $ -0.25213154341695D4,-0.14874640856724,0.31806110878444 / DATA (A(I),I=55,56) /3.5,3.5/ DATA (B(I),I=55,56) /0.85,0.95/ DATA (BB(I),I=55,56) /0.2,0.2/ DATA (Alfa(I),I=52,54) /20,20,20/ DATA (Belt(I),I=52,56) /150,150,250,0.3,0.3/ DATA (Gama(I),I=52,54) /1.21,1.21,1.25/ DATA (El(I),I=52,54) /1,1,1/ DATA (BC(I),I=55,56) /28,32/ DATA (BD(I),I=55,56) /700,800/ DATA (BA(I),I=55,56) /0.32,0.32/ DRWHTT=0. DO 100 I=1,7 DRWHTT=DRWHTT+CN(I)*t(I)*(t(I)-1)*Dx**D(I)*tx**(t(I)-2) 100 CONTINUE DO 200 I=8,51 DRWHTT=DRWHTT+CN(I)*t(I)*(t(I)-1)* $ Dx**D(I)*Tx**(t(I)-2)*DEXP(-Dx**C(I)) 200 CONTINUE DO 300 I=52,54 F=-Alfa(I)*(Dx-El(I))**2-Belt(I)*(Tx-Gama(I))**2 G=(t(I)/Tx-2*Belt(I)*(Tx-Gama(I)))**2 $ -t(I)/(Tx*Tx)-2*Belt(I) DRWHTT=DRWHTT+CN(I)*Dx**D(I)*Tx**t(I)*DEXP(F)*G 300 CONTINUE DO 400 I=55,56 Fal=DEXP(-BC(I)*(Dx-1)**2-BD(I)*(Tx-1)**2) Sita=(1-Tx)+BA(I)*((Dx-1)**2)**(0.5/Belt(I)) Delta=Sita*Sita+BB(I)*((Dx-1)**2)**A(I) Deltat=-2.*Sita*B(I)*Delta**(B(I)-1) Falt=-2*BD(I)*(Tx-1)*Fal Dltatt=2*B(I)*Delta**(B(I)-1) $ +4*Sita*Sita*B(I)*(B(I)-1)*Delta*(B(I)-2) Faltt=2*(2*BD(I)*(Tx-1)*(Tx-1)-1)*BD(I)*Fal G=Dltatt*Fal+2.*Deltat*Falt+Delta**B(I)*Faltt DRWHTT=DRWHTT+CN(I)*Dx*G 400 CONTINUE RETURN END C =============================================================== C ** Differential of the Helmholtz Free Energy of ** C ** the Residual Part to Dimensionless Volume V OF WATER ** C =============================================================== SUBROUTINE D_R_WHV(Dx,Tx,DRWHV) DOUBLE PRECISION Tx,Dx,DRWHV REAL C(8:51),D(1:54),t(1:54) REAL A(55:56),B(55:56),BB(55:56) REAL Alfa(52:54),Belt(52:56),Gama(52:54),El(52:54) REAL BC(55:56),BD(55:56),BA(55:56) DOUBLE PRECISION CN(1:56) DOUBLE PRECISION F,Delta,Sita,Fal,G,Deltav,Falv DATA (C(I),I=8,51) /1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2, $ 2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,4,6,6,6,6/ DATA (D(I),I=1,54) /1,1,1,2,2,3,4,1,1,1,2,2,3,4,4,5,7,9,10,11, $ 13,15,1,2,2,2,3,4,4,4,5,6,6,7,9,9,9,9,9,10,10,12,3,4,4,5, $ 14,3,6,6,6,3,3,3/ DATA (t(I),I=1,54) /-0.5,0.875,1,0.5,0.75,0.375,1,4,6,12,1,5,4, $ 2,13,9,3,4,11,4,13,1,7,1,9,10,10,3,7,10,10,6,10,10,1,2,3,4, $ 8,6,9,8,16,22,23,23,10,50,44,46,50,0,1,4/ DATA (CN(I),I=1,56) /0.12533547935523D-1,0.78957634722828D+1, $ -0.87803203303561D+1,0.31802509345418,-0.26145533859358, $ -0.78199751687981D-2,0.88089493102134D-2,-0.66856572307965, $ 0.20433810950965,-0.66212605039687D-4,-0.19232721156002, $ -0.25709043003438,0.16074868486251,-0.40092828925807D-1, $ 0.39343422603254D-6,-0.75941377088144D-5,0.56250979351888D-3, $ -0.15608652257135D-4,0.11537996422951D-8,0.36582165144204D-6, $ -0.13251180074668D-11,-0.62639586912454D-9,-0.10793600908932, $ 0.17611491008752D-1,0.22132295167546,-0.40247669763528, $ 0.58083399985759,0.49969146990806D-2,-0.31358700712549D-1, $ -0.74315929710341,0.47807329915480,0.20527940895948D-1, $ -0.13636435110343,0.14180634400617D-1,0.83326504880713D-2, $ -0.29052336009585D-1,0.38615085574206D-1,-0.20393486513704D-1, $ -0.16554050063734D-2,0.19955571979541D-2,0.15870308324157D-3, $ -0.16388568342530D-4,0.43613615723811D-1,0.34994005463765D-1, $ -0.76788197844621D-1,0.22446277332006D-1,-0.62689710414685D-4, $ -0.55711118565645D-9,-0.19905718354408,0.31777497330738, $ -0.11841182425981,-0.31306260323435D2,0.31546140237781D2, $ -0.25213154341695D4,-0.14874640856724,0.31806110878444 / DATA (A(I),I=55,56) /3.5,3.5/ DATA (B(I),I=55,56) /0.85,0.95/ DATA (BB(I),I=55,56) /0.2,0.2/ DATA (Alfa(I),I=52,54) /20,20,20/ DATA (Belt(I),I=52,56) /150,150,250,0.3,0.3/ DATA (Gama(I),I=52,54) /1.21,1.21,1.25/ DATA (El(I),I=52,54) /1,1,1/ DATA (BC(I),I=55,56) /28,32/ DATA (BD(I),I=55,56) /700,800/ DATA (BA(I),I=55,56) /0.32,0.32/ DRWHV=0.0 DO 100 I=1,7 DRWHV=DRWHV+CN(I)*D(I)*Dx**(D(I)-1)*Tx**t(I) 100 CONTINUE DO 200 I=8,51 DRWHV=DRWHV+CN(I)*DEXP(-Dx**C(I)) $ *(Dx**(D(I)-1)*Tx**t(I)*(D(I)-C(I)*Dx**C(I))) 200 CONTINUE DO 300 I=52,54 F=-Alfa(I)*(Dx-El(I))**2-Belt(I)*(Tx-Gama(I))**2 G=D(I)/Dx-2.*Alfa(I)*(Dx-El(I)) DRWHV=DRWHV+CN(I)*Dx**D(I)*Tx**t(I)*DEXP(F)*G 300 CONTINUE DO 400 I=55,56 Fal=DEXP(-BC(I)*(Dx-1)**2-BD(I)*(Tx-1)**2) Sita=(1-Tx)+BA(I)*((Dx-1)**2)**(0.5/Belt(I)) Delta=Sita*Sita+BB(I)*((Dx-1)**2)**A(I) Deltav=(Dx-1)*(BA(I)*Sita*2/Belt(I)* $ ((Dx-1)**2)**(0.5/Belt(I)-1) $ +2*BB(I)*A(I)*((Dx-1)**2)**(A(I)-1)) Deltav=B(I)*Delta**(B(I)-1)*Deltav Falv=-2*BC(I)*(Dx-1)*Fal DRWHV=DRWHV+ $ CN(I)*(Delta**B(I)*(Fal+Dx*Falv)+Deltav*Dx*Fal) 400 CONTINUE RETURN END C =============================================================== C ** Differential of the Helmholtz Free Energy of ** C ** the Residual Part to Dimensionless Volume VV OF WATER ** C =============================================================== SUBROUTINE D_R_WHVV(Dx,Tx,DRWHVV) DOUBLE PRECISION Tx,Dx,DRWHVV REAL C(8:51),D(1:54),t(1:54) REAL A(55:56),B(55:56),BB(55:56) REAL Alfa(52:54),Belt(52:56),Gama(52:54),El(52:54) REAL BC(55:56),BD(55:56),BA(55:56) DOUBLE PRECISION CN(1:56) DOUBLE PRECISION F,Delta,Sita,Fal,G,Deltav,Falv,Dltavv,Falvv DATA (C(I),I=8,51) /1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2, $ 2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,4,6,6,6,6/ DATA (D(I),I=1,54) /1,1,1,2,2,3,4,1,1,1,2,2,3,4,4,5,7,9,10,11, $ 13,15,1,2,2,2,3,4,4,4,5,6,6,7,9,9,9,9,9,10,10,12,3,4,4,5, $ 14,3,6,6,6,3,3,3/ DATA (t(I),I=1,54) /-0.5,0.875,1,0.5,0.75,0.375,1,4,6,12,1,5,4, $ 2,13,9,3,4,11,4,13,1,7,1,9,10,10,3,7,10,10,6,10,10,1,2,3,4, $ 8,6,9,8,16,22,23,23,10,50,44,46,50,0,1,4/ DATA (CN(I),I=1,56) /0.12533547935523D-1,0.78957634722828D+1, $ -0.87803203303561D+1,0.31802509345418,-0.26145533859358, $ -0.78199751687981D-2,0.88089493102134D-2,-0.66856572307965, $ 0.20433810950965,-0.66212605039687D-4,-0.19232721156002, $ -0.25709043003438,0.16074868486251,-0.40092828925807D-1, $ 0.39343422603254D-6,-0.75941377088144D-5,0.56250979351888D-3, $ -0.15608652257135D-4,0.11537996422951D-8,0.36582165144204D-6, $ -0.13251180074668D-11,-0.62639586912454D-9,-0.10793600908932, $ 0.17611491008752D-1,0.22132295167546,-0.40247669763528, $ 0.58083399985759,0.49969146990806D-2,-0.31358700712549D-1, $ -0.74315929710341,0.47807329915480,0.20527940895948D-1, $ -0.13636435110343,0.14180634400617D-1,0.83326504880713D-2, $ -0.29052336009585D-1,0.38615085574206D-1,-0.20393486513704D-1, $ -0.16554050063734D-2,0.19955571979541D-2,0.15870308324157D-3, $ -0.16388568342530D-4,0.43613615723811D-1,0.34994005463765D-1, $ -0.76788197844621D-1,0.22446277332006D-1,-0.62689710414685D-4, $ -0.55711118565645D-9,-0.19905718354408,0.31777497330738, $ -0.11841182425981,-0.31306260323435D2,0.31546140237781D2, $ -0.25213154341695D4,-0.14874640856724,0.31806110878444 / DATA (A(I),I=55,56) /3.5,3.5/ DATA (B(I),I=55,56) /0.85,0.95/ DATA (BB(I),I=55,56) /0.2,0.2/ DATA (Alfa(I),I=52,54) /20,20,20/ DATA (Belt(I),I=52,56) /150,150,250,0.3,0.3/ DATA (Gama(I),I=52,54) /1.21,1.21,1.25/ DATA (El(I),I=52,54) /1,1,1/ DATA (BC(I),I=55,56) /28,32/ DATA (BD(I),I=55,56) /700,800/ DATA (BA(I),I=55,56) /0.32,0.32/ DRWHVV=0.0 DO 100 I=1,7 DRWHVV=DRWHVV+CN(I)*D(I)*(D(I)-1)*Dx**(D(I)-2)*Tx**t(I) 100 CONTINUE DO 200 I=8,51 F=Dx**(D(I)-2)*Tx**t(I)* $ ((D(I)-C(I)*Dx**C(I))*(D(I)-1-C(I)*Dx**C(I)) $ -C(I)*C(I)*Dx**C(I)) DRWHVV=DRWHVV+CN(I)*DEXP(-Dx**C(I))*F 200 CONTINUE DO 300 I=52,54 F=-Alfa(I)*(Dx-El(I))**2-Belt(I)*(Tx-Gama(I))**2 G=-2*Alfa(I)*Dx**D(I) $ +4*Alfa(I)*Alfa(I)*Dx**D(I)*(Dx-El(I))*(Dx-El(I)) $ -4*D(I)*Alfa(I)*Dx**(D(I)-1)*(Dx-El(I)) $ +D(I)*(D(I)-1)*Dx**(D(I)-2) DRWHVV=DRWHVV+CN(I)*Tx**t(I)*DEXP(F)*G 300 CONTINUE DO 400 I=55,56 Fal=DEXP(-BC(I)*(Dx-1)**2-BD(I)*(Tx-1)**2) Sita=(1-Tx)+BA(I)*((Dx-1)**2)**(0.5/Belt(I)) Delta=Sita*Sita+BB(I)*((Dx-1)**2)**A(I) Deltav=(Dx-1)*(BA(I)*Sita*2/Belt(I)* $ ((Dx-1)**2)**(0.5/Belt(I)-1) $ +2*BB(I)*A(I)*((Dx-1)**2)**(A(I)-1)) Dltavv=1/(Dx-1)*Deltav+(Dx-1)*(Dx-1)*( $ 4*BB(I)*A(I)*(A(I)-1)*((Dx-1)**2)**(A(I)-2) $ +2*BA(I)**2*(1/Belt(I))**2* $ (((Dx-1)**2)**(0.5/Belt(I)-1))**2 $ +BA(I)*Sita*4/Belt(I)*(0.5/Belt(I)-1)*((Dx-1)**2)** $ (0.5/Belt(I)-2)) Dltavv=B(I)*(Delta**(B(I)-1)*Dltavv $ +(B(I)-1)*Delta**(B(I)-2)*Deltav*Deltav) Deltav=B(I)*Delta**(B(I)-1)*Deltav Falv=-2*BC(I)*(Dx-1)*Fal Falvv=(2*BC(I)*(Dx-1)**2-1)*2*BD(I)*Fal F=Delta**B(I)*(2*Falv+Dx*Falvv) $ +2*Deltav*(Fal+Dx*Falv)+Dltavv*Dx*Fal DRWHVV=DRWHVV+CN(I)*F 400 CONTINUE RETURN END C =============================================================== C ** Differential of the Helmholtz Free Energy of ** C ** the Residual Part to Dimensionless Temperature and ** C ** Volume VV OF WATER ** C =============================================================== SUBROUTINE D_R_WHTV(Dx,Tx,DRWHTV) DOUBLE PRECISION Tx,Dx,DRWHTV REAL C(8:51),D(1:54),t(1:54) REAL A(55:56),B(55:56),BB(55:56) REAL Alfa(52:54),Belt(52:56),Gama(52:54),El(52:54) REAL BC(55:56),BD(55:56),BA(55:56) DOUBLE PRECISION CN(1:56) DOUBLE PRECISION F,Delta,Sita,Fal,G DOUBLE PRECISION Dltat,Dltav,Dlta,Dltatv,Falv,Falt,Faltv C CALL COE_W(C,D,t,CN,Alfa,Belt,Gama,El,A,B,BB,BC,BD,BA) DATA (C(I),I=8,51) /1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2, $ 2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,4,6,6,6,6/ DATA (D(I),I=1,54) /1,1,1,2,2,3,4,1,1,1,2,2,3,4,4,5,7,9,10,11, $ 13,15,1,2,2,2,3,4,4,4,5,6,6,7,9,9,9,9,9,10,10,12,3,4,4,5, $ 14,3,6,6,6,3,3,3/ DATA (t(I),I=1,54) /-0.5,0.875,1,0.5,0.75,0.375,1,4,6,12,1,5,4, $ 2,13,9,3,4,11,4,13,1,7,1,9,10,10,3,7,10,10,6,10,10,1,2,3,4, $ 8,6,9,8,16,22,23,23,10,50,44,46,50,0,1,4/ DATA (CN(I),I=1,56) /0.12533547935523D-1,0.78957634722828D+1, $ -0.87803203303561D+1,0.31802509345418,-0.26145533859358, $ -0.78199751687981D-2,0.88089493102134D-2,-0.66856572307965, $ 0.20433810950965,-0.66212605039687D-4,-0.19232721156002, $ -0.25709043003438,0.16074868486251,-0.40092828925807D-1, $ 0.39343422603254D-6,-0.75941377088144D-5,0.56250979351888D-3, $ -0.15608652257135D-4,0.11537996422951D-8,0.36582165144204D-6, $ -0.13251180074668D-11,-0.62639586912454D-9,-0.10793600908932, $ 0.17611491008752D-1,0.22132295167546,-0.40247669763528, $ 0.58083399985759,0.49969146990806D-2,-0.31358700712549D-1, $ -0.74315929710341,0.47807329915480,0.20527940895948D-1, $ -0.13636435110343,0.14180634400617D-1,0.83326504880713D-2, $ -0.29052336009585D-1,0.38615085574206D-1,-0.20393486513704D-1, $ -0.16554050063734D-2,0.19955571979541D-2,0.15870308324157D-3, $ -0.16388568342530D-4,0.43613615723811D-1,0.34994005463765D-1, $ -0.76788197844621D-1,0.22446277332006D-1,-0.62689710414685D-4, $ -0.55711118565645D-9,-0.19905718354408,0.31777497330738, $ -0.11841182425981,-0.31306260323435D2,0.31546140237781D2, $ -0.25213154341695D4,-0.14874640856724,0.31806110878444 / DATA (A(I),I=55,56) /3.5,3.5/ DATA (B(I),I=55,56) /0.85,0.95/ DATA (BB(I),I=55,56) /0.2,0.2/ DATA (Alfa(I),I=52,54) /20,20,20/ DATA (Belt(I),I=52,56) /150,150,250,0.3,0.3/ DATA (Gama(I),I=52,54) /1.21,1.21,1.25/ DATA (El(I),I=52,54) /1,1,1/ DATA (BC(I),I=55,56) /28,32/ DATA (BD(I),I=55,56) /700,800/ DATA (BA(I),I=55,56) /0.32,0.32/ DRWHTV=0.0 DO 100 I=1,7 DRWHTV=DRWHTV+CN(I)*D(I)*t(I)*Dx**(D(I)-1)*Tx**(t(I)-1) 100 CONTINUE DO 200 I=8,51 DRWHTV=DRWHTV+CN(I)*t(I)*Dx**(D(I)-1)*Tx**(t(I)-1)* $ (D(I)-C(I)*Dx**C(I))*DEXP(-Dx**C(I)) 200 CONTINUE DO 300 I=52,54 F=-Alfa(I)*(Dx-El(I))**2-Belt(I)*(Tx-Gama(I))**2 G=(D(I)/Dx-2*Alfa(I)*(Dx-El(I)))* $ (t(I)/Tx-2*Belt(I)*(Tx-Gama(I))) DRWHTV=DRWHTV+CN(I)*Dx**D(I)*Tx**t(I)*DEXP(F)*G 300 CONTINUE DO 400 I=55,56 Fal=DEXP(-BC(I)*(Dx-1)**2-BD(I)*(Tx-1)**2) Sita=(1-Tx)+BA(I)*((Dx-1)**2)**(0.5/Belt(I)) Delta=Sita*Sita+BB(I)*((Dx-1)**2)**A(I) Dltat=-2.*Sita*B(I)*Delta**(B(I)-1) Falt=-2*BD(I)*(Tx-1)*Fal Dlta=(Dx-1)*(BA(I)*Sita*2/Belt(I)* $ ((Dx-1)**2)**(0.5/Belt(I)-1) $ +2*BB(I)*A(I)*((Dx-1)**2)**(A(I)-1)) Dltav=B(I)*Delta**(B(I)-1)*Dlta Falv=-2*BC(I)*(Dx-1)*Fal Dltatv=-BA(I)*B(I)*2/Belt(I)*(Dx-1)* $ ((Dx-1)**2)**(0.5/Belt(I)-1)- $ 2*Sita*B(I)*(B(I)-1)*Delta**(B(I)-2)*Dlta Faltv=4*BC(I)*BD(I)*(Dx-1)*(Tx-1)*Fal G=Delta**B(I)*(Falt+Dx*Faltv)+Dx*Dltav*Falt+ $ Dltat*(Fal+Dx*Falv)+Dltatv*Dx*Fal DRWHTV=DRWHTV+CN(I)*G 400 CONTINUE RETURN END C =============================================================== C ** Helmholtz Free Energy of Ideal-Ammonia Properties ** C =============================================================== SUBROUTINE H_IAHG(D,T,FA) DOUBLE PRECISION T,D,FA DOUBLE PRECISION A(5) DATA (A(J),J=1,5) /-15.815020,4.255726,11.474340,-1.296211, $ 0.5706757/ FA=DLOG(D)+A(1)+A(2)*T-DLOG(T)+A(3)*T**(1./3.)+A(4)*T**(-1.5) $ +A(5)*T**(-1.75) RETURN END C ========================================================== C ** The Helmholtz Free Energy of The residual part ** C ** OF AMMONIA ** C ========================================================== SUBROUTINE H_RAHG(Dx,Tx,FA) DOUBLE PRECISION A(1:21) REAL t(1:21),D(1:21) DOUBLE PRECISION Tx,Dx,FA,FA1,FA2,FA3,FA4 DATA (A(I),I=1,21) /0.4554431D-1,0.7238548D0,0.1229470D-1, $ -0.1858814D+1,0.2141882D-10,-0.1430020D-1,0.3441324D0, $ -0.2873571D0,0.2352589D-4,-0.3497111D-1,0.2397852D-1, $ 0.1831117D-2,-0.4085375D-1,0.2379275D0,-0.3548972D-1, $ -0.1823729D0,0.2281556D-1,-0.6663444D-2,-0.8847486D-2, $ 0.2272635D-2,-0.5588655D-3/ DATA(t(I),I=1,21) /-0.5,0.5,1.,1.5,3.,0.,3.,4.,4.,5.,3., $ 5.,6.,8.,8.,10.,10.,5.,7.5,15.,30./ DATA(D(I),I=1,21) /2.,1.,4.,1.,15.,3.,3.,1.,8.,2.,1.,8., $ 1.,2.,3.,2.,4.,3.,1.,2.,4./ FA1=A(1)*Tx**t(1)*Dx**D(1)+A(2)*Tx**t(2)*Dx**D(2)+ $ A(3)*Tx**t(3)*Dx**D(3)+A(4)*Tx**t(4)*Dx**D(4)+ $ A(5)*Tx**t(5)*Dx**D(5) FA2=A(6)*Tx**t(6)*Dx**D(6)+A(7)*Tx**t(7)*Dx**D(7)+ $ A(8)*Tx**t(8)*Dx**D(8)+A(9)*Tx**t(9)*Dx**D(9)+ $ A(10)*Tx**t(10)*Dx**D(10) FA3=A(11)*Tx**t(11)*Dx**D(11)+A(12)*Tx**t(12)*Dx**D(12)+ $ A(13)*Tx**t(13)*Dx**D(13)+A(14)*Tx**t(14)*Dx**D(14)+ $ A(15)*Tx**t(15)*Dx**D(15)+A(16)*Tx**t(16)*Dx**D(16)+ $ A(17)*Tx**t(17)*Dx**D(17) FA4=A(18)*Tx**t(18)*Dx**D(18)+A(19)*Tx**t(19)*Dx**D(19)+ $ A(20)*Tx**t(20)*Dx**D(20)+A(21)*Tx**t(21)*Dx**D(21) FA=FA1+DEXP(-Dx)*FA2+DEXP(-Dx*Dx)*FA3+DEXP(-Dx*Dx*Dx)*FA4 RETURN END C =============================================================== C ** Differential of the Helmholtz Free Energy of ** C **the Residual Part to Dimensionless Temperature T OF AMMONIA** C =============================================================== SUBROUTINE D_R_AHT(Dx,Tx,DRAHT) DOUBLE PRECISION A(1:21) REAL t(1:21),D(1:21) DOUBLE PRECISION Tx,Dx,DRAHT,FA1,FA2,FA3,FA4 DATA (A(I),I=1,21) /0.4554431D-1,0.7238548D0,0.1229470D-1, $ -0.1858814D+1,0.2141882D-10,-0.1430020D-1,0.3441324D0, $ -0.2873571D0,0.2352589D-4,-0.3497111D-1,0.2397852D-1, $ 0.1831117D-2,-0.4085375D-1,0.2379275D0,-0.3548972D-1, $ -0.1823729D0,0.2281556D-1,-0.6663444D-2,-0.8847486D-2, $ 0.2272635D-2,-0.5588655D-3/ DATA(t(I),I=1,21) /-0.5,0.5,1.,1.5,3.,0.,3.,4.,4.,5.,3., $ 5.,6.,8.,8.,10.,10.,5.,7.5,15.,30./ DATA(D(I),I=1,21) /2.,1.,4.,1.,15.,3.,3.,1.,8.,2.,1.,8., $ 1.,2.,3.,2.,4.,3.,1.,2.,4./ FA1=A(1)*t(1)*Tx**(t(1)-1.)*Dx**D(1)+ $ A(2)*t(2)*Tx**(t(2)-1.)*Dx**D(2)+ $ A(3)*t(3)*Tx**(t(3)-1.)*Dx**D(3)+ $ A(4)*t(4)*Tx**(t(4)-1.)*Dx**D(4)+ $ A(5)*t(5)*Tx**(t(5)-1.)*Dx**D(5) FA2=A(6)*t(6)*Tx**(t(6)-1.)*Dx**D(6)+ $ A(7)*t(7)*Tx**(t(7)-1.)*Dx**D(7)+ $ A(8)*t(8)*Tx**(t(8)-1.)*Dx**D(8)+ $ A(9)*t(9)*Tx**(t(9)-1.)*Dx**D(9)+ $ A(10)*t(10)*Tx**(t(10)-1.)*Dx**D(10) FA3=A(11)*t(11)*Tx**(t(11)-1.)*Dx**D(11)+ $ A(12)*t(12)*Tx**(t(12)-1.)*Dx**D(12)+ $ A(13)*t(13)*Tx**(t(13)-1.)*Dx**D(13)+ $ A(14)*t(14)*Tx**(t(14)-1.)*Dx**D(14)+ $ A(15)*t(15)*Tx**(t(15)-1.)*Dx**D(15)+ $ A(16)*t(16)*Tx**(t(16)-1.)*Dx**D(16)+ $ A(17)*t(17)*Tx**(t(17)-1.)*Dx**D(17) FA4=A(18)*t(18)*Tx**(t(18)-1.)*Dx**D(18)+ $ A(19)*t(19)*Tx**(t(19)-1.)*Dx**D(19)+ $ A(20)*t(20)*Tx**(t(20)-1.)*Dx**D(20)+ $ A(21)*t(21)*Tx**(t(21)-1.)*Dx**D(21) DRAHT=FA1+DEXP(-Dx)*FA2+DEXP(-Dx*Dx)*FA3+DEXP(-Dx*Dx*Dx)*FA4 RETURN END C =============================================================== C ** Differential of the Helmholtz Free Energy of the Residual ** C ** Part to Dimensionless Temperature TT OF AMMONIA ** C =============================================================== SUBROUTINE D_R_AHTT(Dx,Tx,DRAHTT) DOUBLE PRECISION A(1:21) REAL t(1:21),D(1:21) DOUBLE PRECISION Tx,Dx,DRAHTT,FA1,FA2,FA3,FA4 DATA (A(I),I=1,21) /0.4554431D-1,0.7238548D0,0.1229470D-1, $ -0.1858814D+1,0.2141882D-10,-0.1430020D-1,0.3441324D0, $ -0.2873571D0,0.2352589D-4,-0.3497111D-1,0.2397852D-1, $ 0.1831117D-2,-0.4085375D-1,0.2379275D0,-0.3548972D-1, $ -0.1823729D0,0.2281556D-1,-0.6663444D-2,-0.8847486D-2, $ 0.2272635D-2,-0.5588655D-3/ DATA(t(I),I=1,21) /-0.5,0.5,1.,1.5,3.,0.,3.,4.,4.,5.,3., $ 5.,6.,8.,8.,10.,10.,5.,7.5,15.,30./ DATA(D(I),I=1,21) /2.,1.,4.,1.,15.,3.,3.,1.,8.,2.,1.,8., $ 1.,2.,3.,2.,4.,3.,1.,2.,4./ FA1=A(1)*t(1)*(t(1)-1.)*Tx**(t(1)-2.)*Dx**D(1)+ $ A(2)*t(2)*(t(2)-1.)*Tx**(t(2)-2.)*Dx**D(2)+ $ A(3)*t(3)*(t(3)-1.)*Tx**(t(3)-2.)*Dx**D(3)+ $ A(4)*t(4)*(t(4)-1.)*Tx**(t(4)-2.)*Dx**D(4)+ $ A(5)*t(5)*(t(5)-1.)*Tx**(t(5)-2.)*Dx**D(5) FA2=A(6)*t(6)*(t(6)-1.)*Tx**(t(6)-2.)*Dx**D(6)+ $ A(7)*t(7)*(t(7)-1.)*Tx**(t(7)-2.)*Dx**D(7)+ $ A(8)*t(8)*(t(8)-1.)*Tx**(t(8)-2.)*Dx**D(8)+ $ A(9)*t(9)*(t(9)-1.)*Tx**(t(9)-2.)*Dx**D(9)+ $ A(10)*t(10)*(t(10)-1.)*Tx**(t(10)-2.)*Dx**D(10) FA3=A(11)*t(11)*(t(11)-1.)*Tx**(t(11)-2.)*Dx**D(11)+ $ A(12)*t(12)*(t(12)-1.)*Tx**(t(12)-2.)*Dx**D(12)+ $ A(13)*t(13)*(t(13)-1.)*Tx**(t(13)-2.)*Dx**D(13)+ $ A(14)*t(14)*(t(14)-1.)*Tx**(t(14)-2.)*Dx**D(14)+ $ A(15)*t(15)*(t(15)-1.)*Tx**(t(15)-2.)*Dx**D(15)+ $ A(16)*t(16)*(t(16)-1.)*Tx**(t(16)-2.)*Dx**D(16)+ $ A(17)*t(17)*(t(17)-1.)*Tx**(t(17)-2.)*Dx**D(17) FA4=A(18)*t(18)*(t(18)-1.)*Tx**(t(18)-2.)*Dx**D(18)+ $ A(19)*t(19)*(t(19)-1.)*Tx**(t(19)-2.)*Dx**D(19)+ $ A(20)*t(20)*(t(20)-1.)*Tx**(t(20)-2.)*Dx**D(20)+ $ A(21)*t(21)*(t(21)-1.)*Tx**(t(21)-2.)*Dx**D(21) DRAHTT=FA1+DEXP(-Dx)*FA2+DEXP(-Dx*Dx)*FA3+DEXP(-Dx*Dx*Dx)*FA4 RETURN END C =============================================================== C ** Differential of the Helmholtz Free Energy of the Residual ** C ** Part to Dimensionless Volume V OF AMMONIA ** C =============================================================== SUBROUTINE D_R_AHV(Dx,Tx,DRAHV) DOUBLE PRECISION A(1:21) REAL t(1:21),D(1:21) DOUBLE PRECISION Tx,Dx,DRAHV,FA1,FA2,FA3,FA4 DATA (A(I),I=1,21) /0.4554431D-1,0.7238548D0,0.1229470D-1, $ -0.1858814D+1,0.2141882D-10,-0.1430020D-1,0.3441324D0, $ -0.2873571D0,0.2352589D-4,-0.3497111D-1,0.2397852D-1, $ 0.1831117D-2,-0.4085375D-1,0.2379275D0,-0.3548972D-1, $ -0.1823729D0,0.2281556D-1,-0.6663444D-2,-0.8847486D-2, $ 0.2272635D-2,-0.5588655D-3/ DATA(t(I),I=1,21) /-0.5,0.5,1.,1.5,3.,0.,3.,4.,4.,5.,3., $ 5.,6.,8.,8.,10.,10.,5.,7.5,15.,30./ DATA(D(I),I=1,21) /2.,1.,4.,1.,15.,3.,3.,1.,8.,2.,1.,8., $ 1.,2.,3.,2.,4.,3.,1.,2.,4./ FA1=A(1)*Tx**t(1)*D(1)*Dx**(D(1)-1.)+ $ A(2)*Tx**t(2)*D(2)*Dx**(D(2)-1.)+ $ A(3)*Tx**t(3)*D(3)*Dx**(D(3)-1.)+ $ A(4)*Tx**t(4)*D(4)*Dx**(D(4)-1.)+ $ A(5)*Tx**t(5)*D(5)*Dx**(D(5)-1.) FA2=A(6)*Tx**t(6)*(D(6)*Dx**(D(6)-1.)-Dx**D(6))+ $ A(7)*Tx**t(7)*(D(7)*Dx**(D(7)-1.)-Dx**D(7))+ $ A(8)*Tx**t(8)*(D(8)*Dx**(D(8)-1.)-Dx**D(8))+ $ A(9)*Tx**t(9)*(D(9)*Dx**(D(9)-1.)-Dx**D(9))+ $ A(10)*Tx**t(10)*(D(10)*Dx**(D(10)-1.)-Dx**D(10)) FA3=A(11)*Tx**t(11)*(D(11)*Dx**(D(11)-1.)-2.*Dx**(D(11)+1.))+ $ A(12)*Tx**t(12)*(D(12)*Dx**(D(12)-1.)-2.*Dx**(D(12)+1.))+ $ A(13)*Tx**t(13)*(D(13)*Dx**(D(13)-1.)-2.*Dx**(D(13)+1.))+ $ A(14)*Tx**t(14)*(D(14)*Dx**(D(14)-1.)-2.*Dx**(D(14)+1.))+ $ A(15)*Tx**t(15)*(D(15)*Dx**(D(15)-1.)-2.*Dx**(D(15)+1.))+ $ A(16)*Tx**t(16)*(D(16)*Dx**(D(16)-1.)-2.*Dx**(D(16)+1.))+ $ A(17)*Tx**t(17)*(D(17)*Dx**(D(17)-1.)-2.*Dx**(D(17)+1.)) FA4=A(18)*Tx**t(18)*(D(18)*Dx**(D(18)-1.)-3.*Dx**(D(18)+2.))+ $ A(19)*Tx**t(19)*(D(19)*Dx**(D(19)-1.)-3.*Dx**(D(19)+2.))+ $ A(20)*Tx**t(20)*(D(20)*Dx**(D(20)-1.)-3.*Dx**(D(20)+2.))+ $ A(21)*Tx**t(21)*(D(21)*Dx**(D(21)-1.)-3.*Dx**(D(21)+2.)) DRAHV=FA1+DEXP(-Dx)*FA2+DEXP(-Dx*Dx)*FA3+DEXP(-Dx*Dx*Dx)*FA4 RETURN END C ================================================================ C ** Differential of the Helmholtz Free Energy of the Residual ** C ** Part to Dimensionless Volume and Temperature TV OF AMMONIA ** C ================================================================ SUBROUTINE D_R_AHTV(Dx,Tx,DRAHTV) DOUBLE PRECISION A(1:21) REAL t(1:21),D(1:21) DOUBLE PRECISION Tx,Dx,DRAHTV,FA1,FA2,FA3,FA4 DATA (A(I),I=1,21) /0.4554431D-1,0.7238548D0,0.1229470D-1, $ -0.1858814D+1,0.2141882D-10,-0.1430020D-1,0.3441324D0, $ -0.2873571D0,0.2352589D-4,-0.3497111D-1,0.2397852D-1, $ 0.1831117D-2,-0.4085375D-1,0.2379275D0,-0.3548972D-1, $ -0.1823729,0.2281556D-1,-0.6663444D-2,-0.8847486D-2, $ 0.2272635D-2,-0.5588655D-3/ DATA(t(I),I=1,21) /-0.5,0.5,1.,1.5,3.,0.,3.,4.,4.,5.,3., $ 5.,6.,8.,8.,10.,10.,5.,7.5,15.,30./ DATA(D(I),I=1,21) /2.,1.,4.,1.,15.,3.,3.,1.,8.,2.,1.,8., $ 1.,2.,3.,2.,4.,3.,1.,2.,4./ FA1=A(1)*t(1)*Tx**(t(1)-1.)*D(1)*Dx**(D(1)-1.)+ $ A(2)*t(2)*Tx**(t(2)-1.)*D(2)*Dx**(D(2)-1.)+ $ A(3)*t(3)*Tx**(t(3)-1.)*D(3)*Dx**(D(3)-1.)+ $ A(4)*t(4)*Tx**(t(4)-1.)*D(4)*Dx**(D(4)-1.)+ $ A(5)*t(5)*Tx**(t(5)-1.)*D(5)*Dx**(D(5)-1.) FA2=A(6)*t(6)*Tx**(t(6)-1.)*(D(6)*Dx**(D(6)-1.)-Dx**D(6))+ $ A(7)*t(7)*Tx**(t(7)-1.)*(D(7)*Dx**(D(7)-1.)-Dx**D(7))+ $ A(8)*t(8)*Tx**(t(8)-1.)*(D(8)*Dx**(D(8)-1.)-Dx**D(8))+ $ A(9)*t(9)*Tx**(t(9)-1.)*(D(9)*Dx**(D(9)-1.)-Dx**D(9))+ $ A(10)*t(10)*Tx**(t(10)-1.)*(D(10)*Dx**(D(10)-1.)-Dx**D(10)) FA3= $ A(11)*t(11)*Tx**(t(11)-1)*(D(11)*Dx**(D(11)-1)-2*Dx**(D(11)+1))+ $ A(12)*t(12)*Tx**(t(12)-1)*(D(12)*Dx**(D(12)-1)-2*Dx**(D(12)+1))+ $ A(13)*t(13)*Tx**(t(13)-1)*(D(13)*Dx**(D(13)-1)-2*Dx**(D(13)+1))+ $ A(14)*t(14)*Tx**(t(14)-1)*(D(14)*Dx**(D(14)-1)-2*Dx**(D(14)+1))+ $ A(15)*t(15)*Tx**(t(15)-1)*(D(15)*Dx**(D(15)-1)-2*Dx**(D(15)+1))+ $ A(16)*t(16)*Tx**(t(16)-1)*(D(16)*Dx**(D(16)-1)-2*Dx**(D(16)+1))+ $ A(17)*t(17)*Tx**(t(17)-1)*(D(17)*Dx**(D(17)-1)-2*Dx**(D(17)+1)) FA4= $ A(18)*t(18)*Tx**(t(18)-1)*(D(18)*Dx**(D(18)-1)-3*Dx**(D(18)+2))+ $ A(19)*t(19)*Tx**(t(19)-1)*(D(19)*Dx**(D(19)-1)-3*Dx**(D(19)+2))+ $ A(20)*t(20)*Tx**(t(20)-1)*(D(20)*Dx**(D(20)-1)-3*Dx**(D(20)+2))+ $ A(21)*t(21)*Tx**(t(21)-1)*(D(21)*Dx**(D(21)-1)-3*Dx**(D(21)+2)) DRAHTV=FA1+DEXP(-Dx)*FA2+DEXP(-Dx*Dx)*FA3+DEXP(-Dx*Dx*Dx)*FA4 RETURN END C =============================================================== C ** Differential of the Helmholtz Free Energy of the Residual ** C ** Part to Dimensionless Volume VX OF AMMONIA ** C =============================================================== SUBROUTINE D_R_AHVV(Dx,Tx,DRAHVV) DOUBLE PRECISION A(1:21) REAL t(1:21),D(1:21) DOUBLE PRECISION Tx,Dx,DRAHVV,FA1,FA2,FA3,FA4 DATA (A(I),I=1,21) /0.4554431D-1,0.7238548,0.1229470D-1, $ -0.1858814D1,0.2141882D-10,-0.1430020D-1,0.3441324D0, $ -0.2873571D0,0.2352589D-4,-0.3497111D-1,0.2397852D-1, $ 0.1831117D-2,-0.4085375D-1,0.2379275D0,-0.3548972D-1, $ -0.1823729,0.2281556D-1,-0.6663444D-2,-0.8847486D-2, $ 0.2272635D-2,-0.5588655D-3/ DATA(t(I),I=1,21) /-0.5,0.5,1.,1.5,3.,0.,3.,4.,4.,5.,3., $ 5.,6.,8.,8.,10.,10.,5.,7.5,15.,30./ DATA(D(I),I=1,21) /2.,1.,4.,1.,15.,3.,3.,1.,8.,2.,1.,8., $ 1.,2.,3.,2.,4.,3.,1.,2.,4./ FA1=A(1)*Tx**t(1)*D(1)*(D(1)-1.)*Dx**(D(1)-2.)+ $ A(2)*Tx**t(2)*D(2)*(D(2)-1.)*Dx**(D(2)-2.)+ $ A(3)*Tx**t(3)*D(3)*(D(3)-1.)*Dx**(D(3)-2.)+ $ A(4)*Tx**t(4)*D(4)*(D(4)-1.)*Dx**(D(4)-2.)+ $ A(5)*Tx**t(5)*D(5)*(D(5)-1.)*Dx**(D(5)-2.) FA2=A(6)*Tx**t(6)*(D(6)*(D(6)-1)*Dx**(D(6)-2) $ -2*D(6)*Dx**(D(6)-1)+DX**D(6)) $ +A(7)*Tx**t(7)*(D(7)*(D(7)-1)*Dx**(D(7)-2) $ -2*D(7)*Dx**(D(7)-1)+DX**D(7)) $ +A(8)*Tx**t(8)*(D(8)*(D(8)-1)*Dx**(D(8)-2) $ -2*D(8)*Dx**(D(8)-1)+DX**D(8)) $ +A(9)*Tx**t(9)*(D(9)*(D(9)-1)*Dx**(D(9)-2) $ -2*D(9)*Dx**(D(9)-1)+DX**D(9)) $ +A(10)*Tx**t(10)*(D(10)*(D(10)-1)*Dx**(D(10)-2) $ -2*D(10)*Dx**(D(10)-1)+DX**D(10)) FA3=A(11)*Tx**T(11)*(D(11)*(D(11)-1)*Dx**(D(11)-2) $ -(4*D(11)+2)*Dx**D(11)+4*Dx**(D(11)+2)) $ +(12)*Tx**T(12)*(D(12)*(D(12)-1)*Dx**(D(12)-2) $ -(4*D(12)+2)*Dx**D(12)+4*Dx**(D(12)+2)) $ +(13)*Tx**T(13)*(D(13)*(D(13)-1)*Dx**(D(13)-2) $ -(4*D(13)+2)*Dx**D(13)+4*Dx**(D(13)+2)) $ +(14)*Tx**T(14)*(D(14)*(D(14)-1)*Dx**(D(14)-2) $ -(4*D(14)+2)*Dx**D(14)+4*Dx**(D(14)+2)) $ +(15)*Tx**T(15)*(D(15)*(D(15)-1)*Dx**(D(15)-2) $ -(4*D(15)+2)*Dx**D(15)+4*Dx**(D(15)+2)) $ +(16)*Tx**T(16)*(D(16)*(D(16)-1)*Dx**(D(16)-2) $ -(4*D(16)+2)*Dx**D(16)+4*Dx**(D(16)+2)) $ +(17)*Tx**T(17)*(D(17)*(D(17)-1)*Dx**(D(17)-2) $ -(4*D(17)+2)*Dx**D(17)+4*Dx**(D(17)+2)) FA4= A(18)*Tx**t(18)*(D(18)*(D(18)-1)*Dx**(D(18)-2) $ -6*(D(18)+1)*Dx**(D(18)+1)+9*Dx**(D(18)+4)) $ +A(19)*Tx**t(19)*(D(19)*(D(19)-1)*Dx**(D(19)-2) $ -6*(D(19)+1)*Dx**(D(19)+1)+9*Dx**(D(19)+4)) $ +A(20)*Tx**t(20)*(D(20)*(D(20)-1)*Dx**(D(20)-2) $ -6*(D(20)+1)*Dx**(D(20)+1)+9*Dx**(D(20)+4)) $ +A(21)*Tx**t(21)*(D(21)*(D(21)-1)*Dx**(D(21)-2) $ -6*(D(21)+1)*Dx**(D(21)+1)+9*Dx**(D(21)+4)) DRAHVV=FA1+DEXP(-Dx)*FA2+DEXP(-Dx*Dx)*FA3+DEXP(-Dx*Dx*Dx)*FA4 RETURN END c ----- The Thermal Properties of Ammonia END, END, END ----- C ========================================================== C ** TNX, VNX ** C ========================================================== SUBROUTINE TN_X(X,TNX,VNX) DOUBLE PRECISION X,TNX,VNX DOUBLE PRECISION Kv,Kt,Alfa,Belt,Tc12,Vc12 DOUBLE PRECISION M1,M2 M1=0.018015268 M2=0.01703026 c ------- Critical Temperature and V --------------- Tc1=647.096 Tc2=405.40 Vc1=1.0/322.0*M1 Vc2=1.0/225.0*M2 c ------ Coefficiences of Reducing Function --------- Kv=1.2395117 Kt=0.9648407 Alfa=1.125455 Belt=0.8978069 c -- Dimensionless T and V For Residual Part of HFG-- Tc12=0.5*Kt*(Tc1+Tc2) Vc12=0.5*Kv*(Vc1+Vc2) TNX=(1.0-X)**2.*Tc1+X**2.*Tc2+2.0*X*(1.0-X**Alfa)*Tc12 VNX=(1.0-x)**2.*Vc1+X**2.*Vc2+2.0*X*(1.0-X**Belt)*Vc12 RETURN END C ====================================================== C ** DIFFERENT OF DTX,DDY ** C ====================================================== SUBROUTINE DTV_X(X,DTNX,DDNX) DOUBLE PRECISION X,DTNX,DDNX DOUBLE PRECISION Kv,Kt,Alfa,Belt,Tc12,Vc12 DOUBLE PRECISION M1,M2 M1=0.018015268 M2=0.01703026 c ------- Critical Temperature and V ---------------- Tc1=647.096 Tc2=405.40 Vc1=1.0/322.0*M1 Vc2=1.0/225.0*M2 c ------ Coefficiences of Reducing Function --------- Kv=1.2395117 Kt=0.9648407 Alfa=1.125455 Belt=0.8978069 c -- Dimensionless T and V For Residual Part of HFG-- Tc12=0.5*Kt*(Tc1+Tc2) Vc12=0.5*Kv*(Vc1+Vc2) DTNX=-2*(1-X)*Tc1+2*X*Tc2+2*(1-(Alfa+1)*X**Alfa)*Tc12 DDNX=-2*(1-X)*Vc1+2*X*Vc2+2*(1-(Belt+1)*X**Belt)*Vc12 RETURN END C ====================================================== C ** T, V ---> Tx,Dx; To,Do ** C ====================================================== SUBROUTINE TRAN_TV(T,V,X,Tx,Dx,To,Do) DOUBLE PRECISION T,V,Tx,Dx,To,Do,X DOUBLE PRECISION Tno,Vno,TNX,VNX c --------------------- Tno and Vno -------------------- Tno=500.0 Vno=1.0/15000.0 c -------Uniform Variables To and Do For Ideal Gas------ To=Tno/T Do=Vno/V CALL TN_X(X,TNX,VNX) Dx=VNX/V Tx=TNX/T RETURN END C ========================================================== C ** Helmholtz Free Energy of Ideal-Gas Properties ** C ========================================================== SUBROUTINE H_IMHG(Do,To,X,FA) DOUBLE PRECISION To,Do,FA1,FA2,FA3,X,FA DOUBLE PRECISION A(1:14) DOUBLE PRECISION Si(4:8),t(12:14) DATA (A(J),J=1,14) /-7.720435,8.649358,3.00632,0.012436,0.97315 $ ,1.27950,0.96956,0.24873,-16.444285,4.036946,-1.0,10.69955 $ ,-1.775436,0.82374034/ DATA (Si(J),J=4,8) /1.666,4.578,10.018,11.964,35.600/ DATA (t(J),J=12,14) /0.3333333333,-1.5,-1.75/ FA1=DLOG(Do) IF (X.GE.1.) THEN FA2=0.0 ELSE FA2=(1.-X)*(A(1)+A(2)*To+A(3)*DLOG(To)+DLOG(1.-X)+ $ A(4)*DLOG(1.0-DEXP(-Si(4)*To))+A(5)*DLOG(1.0-DEXP(-Si(5)*To)) $ +A(6)*DLOG(1.0-DEXP(-Si(6)*To))+A(7)*DLOG(1.0-DEXP(-Si(7)*To)) $ +A(8)*DLOG(1.0-DEXP(-Si(8)*To))) END IF IF (X.LE.0.0) THEN FA3=0.0 ELSE FA3=X*(A(9)+A(10)*To+A(11)*DLOG(To)+DLOG(X)+A(12)*To**t(12) $ +A(13)*To**t(13)+A(14)*To**t(14)) END IF FA=FA1+FA2+FA3 RETURN END C ==================================================== C ** The Helmholtz Free Energy of The residual part ** C ==================================================== SUBROUTINE H_RDHG(Dx,Tx,X,FA) DOUBLE PRECISION A(1:14) REAL tl(1:14),D(1:14),E(2:14) DOUBLE PRECISION Tx,Dx,R,FA1,FA2,FA,X DATA (A(I),I=1,14) /-1.855822D-2,5.258010D-2,3.552874D-10 $ ,5.451379D-6,-5.998546D-13,-3.687808D-6,0.2586192 $ ,-1.368072D-8,1.226146D-2,-7.181443D-2,9.970849D-2 $ ,1.0584086D-3,-0.1963687,-0.7777897/ DATA (tl(I),I=1,14) / 1.5,0.5,6.5,1.75,15.0,6.0,-1.0,4.0 $ ,3.5,0.0,-1.0,8.0,7.5,4.0 / DATA (D(I),I=1,14) / 4.0,5.0,15.0,12.0,12.0,15.0,4.0,15.0 $ ,4.0,5.0,6.0,10.0,6.0,2.0 / DATA (E(I),I=2,14) / 1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0 $ ,2.0,2.0,2.0 / R=0.5248379 FA1=A(2)*DEXP(-Dx**E(2))*Tx**Tl(2)*Dx**D(2) $ +A(3)*DEXP(-Dx**E(3))*Tx**Tl(3)*Dx**D(3) $ +A(4)*DEXP(-Dx**E(4))*Tx**Tl(4)*Dx**D(4) $ +A(5)*DEXP(-Dx**E(5))*Tx**Tl(5)*Dx**D(5) $ +A(6)*DEXP(-Dx**E(6))*Tx**Tl(6)*Dx**D(6) FA2=A(7)*DEXP(-Dx**E(7))*Tx**Tl(7)*Dx**D(7) $ +A(8)*DEXP(-Dx**E(8))*Tx**Tl(8)*Dx**D(8) $ +A(9)*DEXP(-Dx**E(9))*Tx**Tl(9)*Dx**D(9) $ +A(10)*DEXP(-Dx**E(10))*Tx**Tl(10)*Dx**D(10) $ +A(11)*DEXP(-Dx**E(11))*Tx**Tl(11)*Dx**D(11) $ +A(12)*DEXP(-Dx**E(12))*Tx**Tl(12)*Dx**D(12) $ +A(13)*DEXP(-Dx**E(13))*Tx**Tl(13)*Dx**D(13) FA=A(1)*Tx**Tl(1)*Dx**D(1)+FA1+X*FA2+ $ A(14)*X*X*DEXP(-Dx**E(14))*Tx**Tl(14)*Dx**D(14) FA=X*(1-X**R)*FA RETURN END C ============================================================== c ** Differential of Helmholtz Free Energy of Ideal Gas ** c ** to Dimensionless Temperature ** C ============================================================== SUBROUTINE D_I_HT(To,X,DIHT) DOUBLE PRECISION To,DIHT,X DOUBLE PRECISION A(1:14) DOUBLE PRECISION Si(4:8),t(12:14) DATA (A(J),J=1,14) /-7.720435,8.649358,3.00632,0.012436,0.97315 $ ,1.27950,0.96956,0.24873,-16.444285,4.036946,-1.0,10.69955 $ ,-1.775436,0.82374034/ DATA (Si(J),J=4,8) /1.666,4.578,10.018,11.964,35.600/ DATA (t(J),J=12,14) /0.3333333333,-1.5,-1.75/ DIHT=(1.0-X)*(A(2)+A(3)/To $ +A(4)*Si(4)/(1.0-DEXP(-Si(4)*To))*DEXP(-Si(4)*To) $ +A(5)*Si(5)/(1.0-DEXP(-Si(5)*To))*DEXP(-Si(5)*To) $ +A(6)*Si(6)/(1.0-DEXP(-Si(6)*To))*DEXP(-Si(6)*To) $ +A(7)*Si(7)/(1.0-DEXP(-Si(7)*To))*DEXP(-Si(7)*To) $ +A(8)*Si(8)/(1.0-DEXP(-Si(8)*To))*DEXP(-Si(8)*To)) DIHT=DIHT+X*(A(10)+A(11)/To+A(12)*t(12)*To**(t(12)-1.0) $ +A(13)*t(13)*To**(t(13)-1.0) $ +A(14)*t(14)*To**(t(14)-1.0)) RETURN END C ============================================================== c ** Differential of Helmholtz Free Energy of Ideal Gas ** c ** to Dimensionless TT ** C ============================================================== SUBROUTINE D_I_HTT(To,X,DIHTT) DOUBLE PRECISION To,DIHTT,X DOUBLE PRECISION A(1:14) DOUBLE PRECISION Si(4:8),t(12:14) DATA (A(J),J=1,14) /-7.720435,8.649358,3.00632,0.012436,0.97315 $ ,1.27950,0.96956,0.24873,-16.444285,4.036946,-1.0,10.69955 $ ,-1.775436,0.82374034/ DATA (Si(J),J=4,8) /1.666,4.578,10.018,11.964,35.600/ DATA (t(J),J=12,14) /0.3333333333,-1.5,-1.75/ DIHTT=-A(3)/(To*To) $ -A(4)*Si(4)**2*DEXP(-Si(4)*To)/(1.0-DEXP(-Si(4)*To))**2.0 $ -A(5)*Si(5)**2*DEXP(-Si(5)*To)/(1.0-DEXP(-Si(5)*To))**2.0 $ -A(6)*Si(6)**2*DEXP(-Si(6)*To)/(1.0-DEXP(-Si(6)*To))**2.0 $ -A(7)*Si(7)**2*DEXP(-Si(7)*To)/(1.0-DEXP(-Si(7)*To))**2.0 $ -A(8)*Si(8)**2*DEXP(-Si(8)*To)/(1.0-DEXP(-Si(8)*To))**2.0 DIHTT=(1.0-X)*DIHTT+X*(-A(11)/(To*To) $ +A(12)*t(12)*(t(12)-1.0)*To**(t(12)-2.0) $ +A(13)*t(13)*(t(13)-1.0)*To**(t(13)-2.0) $ +A(14)*t(14)*(t(14)-1.0)*To**(t(14)-2.0)) RETURN END C ================================================================ c ** Differential of Departure Function of the Residual Part of ** c ** the Helmholtz Free Energy to Dimensionless Temperature ** C ================================================================ SUBROUTINE D_R_DHT(Dx,Tx,X,DRDHT) DOUBLE PRECISION A(1:14) REAL tl(1:14),D(1:14),E(2:14) DOUBLE PRECISION Tx,Dx,R,DRDHT,X DATA (A(I),I=1,14) /-1.855822D-2,5.258010D-2,3.552874D-10 $ ,5.451379D-6,-5.998546D-13,-3.687808D-6,0.2586192 $ ,-1.368072D-8,1.226146D-2,-7.181443D-2,9.970849D-2 $ ,1.05840860D-3,-0.1963687,-0.7777897/ DATA (tl(I),I=1,14) / 1.5,0.5,6.5,1.75,15.0,6.0,-1.0,4.0 $ ,3.5,0.0,-1.0,8.0,7.5,4.0 / DATA (D(I),I=1,14) / 4.0,5.0,15.0,12.0,12.0,15.0,4.0,15.0 $ ,4.0,5.0,6.0,10.0,6.0,2.0 / DATA (E(I),I=2,14) / 1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0 $ ,2.0,2.0,2.0 / c ------ Coefficiences of Reducing Function --------- R=0.5248379 DRDHT=A(1)*Dx**D(1)*tl(1)*Tx**(tl(1)-1.0) $ +A(2)*DEXP(-Dx**E(2))*Dx**D(2)*tl(2)*Tx**(tl(2)-1.) $ +A(3)*DEXP(-Dx**E(3))*Dx**D(3)*tl(3)*Tx**(tl(3)-1.) $ +A(4)*DEXP(-Dx**E(4))*Dx**D(4)*tl(4)*Tx**(tl(4)-1.) $ +A(5)*DEXP(-Dx**E(5))*Dx**D(5)*tl(5)*Tx**(tl(5)-1.) $ +A(6)*DEXP(-Dx**E(6))*Dx**D(6)*tl(6)*Tx**(tl(6)-1.) DRDHT=DRDHT $ +X*(A(7)*DEXP(-Dx**E(7))*Dx**D(7)*tl(7)*Tx**(tl(7)-1.) $ +A(8)*DEXP(-Dx**E(8))*Dx**D(8)*tl(8)*Tx**(tl(8)-1.) $ +A(9)*DEXP(-Dx**E(9))*Dx**D(9)*tl(9)*Tx**(tl(9)-1.) $ +A(10)*DEXP(-Dx**E(10))*Dx**D(10)*tl(10)*Tx**(tl(10)-1.) $ +A(11)*DEXP(-Dx**E(11))*Dx**D(11)*tl(11)*Tx**(tl(11)-1.) $ +A(12)*DEXP(-Dx**E(12))*Dx**D(12)*tl(12)*Tx**(tl(12)-1.) $ +A(13)*DEXP(-Dx**E(13))*Dx**D(13)*tl(13)*Tx**(tl(13)-1.)) DRDHT=DRDHT+ $ X*X*A(14)*DEXP(-Dx**E(14))*Dx**D(14)*tl(14)*Tx**(tl(14)-1.) DRDHT=X*(1.0-X**R)*DRDHT RETURN END C ================================================================ c ** Differential of Departure Function of the Residual Part of ** c ** the Helmholtz Free Energy to Dimensionless T and X ** C ================================================================ SUBROUTINE D_R_DHTX(Dx,Tx,X,DRDHTX) DOUBLE PRECISION A(1:14) REAL tl(1:14),D(1:14),E(2:14) DOUBLE PRECISION Tx,Dx,R,DRDHTX,X,Xr DOUBLE PRECISION C1,C2,C3 DATA (A(I),I=1,14) /-1.855822D-2,5.258010D-2,3.552874D-10 $ ,5.451379D-6,-5.998546D-13,-3.687808D-6,0.2586192 $ ,-1.368072D-8,1.226146D-2,-7.181443D-2,9.970849D-2 $ ,1.05840860D-3,-0.1963687,-0.7777897/ DATA (tl(I),I=1,14) / 1.5,0.5,6.5,1.75,15.0,6.0,-1.0,4.0 $ ,3.5,0.0,-1.0,8.0,7.5,4.0 / DATA (D(I),I=1,14) / 4.0,5.0,15.0,12.0,12.0,15.0,4.0,15.0 $ ,4.0,5.0,6.0,10.0,6.0,2.0 / DATA (E(I),I=2,14) / 1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0 $ ,2.0,2.0,2.0 / c ------ Coefficiences of Reducing Function --------- R=0.5248379 C1=A(1)*Dx**D(1)*tl(1)*Tx**(tl(1)-1.0) $ +A(2)*DEXP(-Dx**E(2))*Dx**D(2)*tl(2)*Tx**(tl(2)-1.) $ +A(3)*DEXP(-Dx**E(3))*Dx**D(3)*tl(3)*Tx**(tl(3)-1.) $ +A(4)*DEXP(-Dx**E(4))*Dx**D(4)*tl(4)*Tx**(tl(4)-1.) $ +A(5)*DEXP(-Dx**E(5))*Dx**D(5)*tl(5)*Tx**(tl(5)-1.) $ +A(6)*DEXP(-Dx**E(6))*Dx**D(6)*tl(6)*Tx**(tl(6)-1.) C2=A(7)*DEXP(-Dx**E(7))*Dx**D(7)*tl(7)*Tx**(tl(7)-1.) $ +A(8)*DEXP(-Dx**E(8))*Dx**D(8)*tl(8)*Tx**(tl(8)-1.) $ +A(9)*DEXP(-Dx**E(9))*Dx**D(9)*tl(9)*Tx**(tl(9)-1.) $ +A(10)*DEXP(-Dx**E(10))*Dx**D(10)*tl(10)*Tx**(tl(10)-1.) $ +A(11)*DEXP(-Dx**E(11))*Dx**D(11)*tl(11)*Tx**(tl(11)-1.) $ +A(12)*DEXP(-Dx**E(12))*Dx**D(12)*tl(12)*Tx**(tl(12)-1.) $ +A(13)*DEXP(-Dx**E(13))*Dx**D(13)*tl(13)*Tx**(tl(13)-1.) C3=A(14)*DEXP(-Dx**E(14))*Dx**D(14)*tl(14)*Tx**(tl(14)-1.) Xr=X**R DRDHTX=C1*(1.-(R+1.)*Xr)+C2*X*(2.-(R+2.)*Xr) $ +C3*X*X*(3.-(R+3.)*Xr) RETURN END C ================================================================ c ** Differential of Departure Function of the Residual Part of ** c ** the Helmholtz Free Energy to Dimensionless V ** C ================================================================ SUBROUTINE D_R_DHV(Dx,Tx,X,DRDHV) DOUBLE PRECISION A(1:14) REAL tl(1:14),D(1:14),E(2:14) DOUBLE PRECISION Tx,Dx,R,DRDHV,X DATA (A(I),I=1,14) /-1.855822D-2,5.258010D-2,3.552874D-10 $ ,5.451379D-6,-5.998546D-13,-3.687808D-6,0.2586192 $ ,-1.368072D-8,1.226146D-2,-7.181443D-2,9.970849D-2 $ ,1.0584086D-3,-0.1963687,-0.7777897/ DATA (tl(I),I=1,14) / 1.5,0.5,6.5,1.75,15.0,6.0,-1.0,4.0 $ ,3.5,0.0,-1.0,8.0,7.5,4.0 / DATA (D(I),I=1,14) / 4.0,5.0,15.0,12.0,12.0,15.0,4.0,15.0 $ ,4.0,5.0,6.0,10.0,6.0,2.0 / DATA (E(I),I=2,14) / 1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0 $ ,2.0,2.0,2.0 / R=0.5248379 DRDHV=A(1)*Tx**tl(1)*D(1)*Dx**(D(1)-1.0) DO 100 I=2,6 DRDHV=DRDHV+A(I)*Tx**tl(I)*DEXP(-Dx**E(I))*( $ D(I)*Dx**(D(I)-1)-E(I)*Dx**(D(I)+E(I)-1.)) 100 CONTINUE DO 200 I=7,13 DRDHV=DRDHV+X*A(I)*Tx**tl(I)*DEXP(-Dx**E(I))*( $ D(I)*Dx**(D(I)-1)-E(I)*Dx**(D(I)+E(I)-1.)) 200 CONTINUE DRDHV=DRDHV+X*X*A(14)*Tx**tl(14)*DEXP(-Dx**E(14))*( $ D(14)*Dx**(D(14)-1)-E(14)*Dx**(D(14)+E(14)-1.)) DRDHV=X*(1.0-X**R)*DRDHV RETURN END C ================================================================ c ** Differential of Departure Function of the Residual Part of ** c ** the Helmholtz Free Energy to Dimensionless V and X ** C ================================================================ SUBROUTINE D_R_DHVX(Dx,Tx,X,DRDHVX) DOUBLE PRECISION A(1:14) REAL tl(1:14),D(1:14),E(2:14),C1,C2,C3 DOUBLE PRECISION Tx,Dx,R,DRDHVX,X,XR DATA (A(I),I=1,14) /-1.855822D-2,5.258010D-2,3.552874D-10 $ ,5.451379D-6,-5.998546D-13,-3.687808D-6,0.2586192 $ ,-1.368072D-8,1.226146D-2,-7.181443D-2,9.970849D-2 $ ,1.0584086D-3,-0.1963687,-0.7777897/ DATA (tl(I),I=1,14) / 1.5,0.5,6.5,1.75,15.0,6.0,-1.0,4.0 $ ,3.5,0.0,-1.0,8.0,7.5,4.0 / DATA (D(I),I=1,14) / 4.0,5.0,15.0,12.0,12.0,15.0,4.0,15.0 $ ,4.0,5.0,6.0,10.0,6.0,2.0 / DATA (E(I),I=2,14) / 1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0 $ ,2.0,2.0,2.0 / R=0.5248379 XR=X**R C1=A(1)*Tx**tl(1)*D(1)*Dx**(D(1)-1.0) DO 100 I=2,6 C1=C1+A(I)*Tx**tl(I)*DEXP(-Dx**E(I))*( $ D(I)*Dx**(D(I)-1)-E(I)*Dx**(D(I)+E(I)-1.)) 100 CONTINUE DO 200 I=7,13 C2=A(I)*Tx**tl(I)*DEXP(-Dx**E(I))*( $ D(I)*Dx**(D(I)-1)-E(I)*Dx**(D(I)+E(I)-1.)) 200 CONTINUE C3=A(14)*Tx**tl(14)*DEXP(-Dx**E(14))*( $ D(14)*Dx**(D(14)-1)-E(14)*Dx**(D(14)+E(14)-1.)) DRDHVX=(1.-(R+1)*XR)*C1+X*(2.-(R+1)*XR)*C2+X*X*(3.-(R+3)*XR)*C3 RETURN END C ================================================================ c ** Differential of Departure Function of the Residual Part of ** c ** the Helmholtz Free Energy to Mole Fraction of NH3 X ** C ================================================================ SUBROUTINE D_R_DHX(Dx,Tx,X,DRDHX) DOUBLE PRECISION A(1:14) REAL tl(1:14),D(1:14),E(2:14) DOUBLE PRECISION Tx,Dx,R,DRDHX DOUBLE PRECISION C1,C2,C3,Xr,X DATA (A(I),I=1,14) /-1.855822D-2,5.258010D-2,3.552874D-10 $ ,5.451379D-6,-5.998546D-13,-3.687808D-6,0.2586192 $ ,-1.368072D-8,1.226146D-2,-7.181443D-2,9.970849D-2 $ ,1.0584086D-3,-0.1963687,-0.7777897/ DATA (tl(I),I=1,14) / 1.5,0.5,6.5,1.75,15.0,6.0,-1.0,4.0 $ ,3.5,0.0,-1.0,8.0,7.5,4.0 / DATA (D(I),I=1,14) / 4.0,5.0,15.0,12.0,12.0,15.0,4.0,15.0 $ ,4.0,5.0,6.0,10.0,6.0,2.0 / DATA (E(I),I=2,14) / 1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0 $ ,2.0,2.0,2.0 / R=0.5248379 C1=A(1)*Tx**tl(1)*Dx**D(1) $ +A(2)*DEXP(-Dx**E(2))*Tx**Tl(2)*Dx**D(2) $ +A(3)*DEXP(-Dx**E(3))*Tx**Tl(3)*Dx**D(3) $ +A(4)*DEXP(-Dx**E(4))*Tx**Tl(4)*Dx**D(4) $ +A(5)*DEXP(-Dx**E(5))*Tx**Tl(5)*Dx**D(5) $ +A(6)*DEXP(-Dx**E(6))*Tx**Tl(6)*Dx**D(6) C2=A(7)*DEXP(-Dx**E(7))*Tx**Tl(7)*Dx**D(7) $ +A(8)*DEXP(-Dx**E(8))*Tx**Tl(8)*Dx**D(8) $ +A(9)*DEXP(-Dx**E(9))*Tx**Tl(9)*Dx**D(9) $ +A(10)*DEXP(-Dx**E(10))*Tx**Tl(10)*Dx**D(10) $ +A(11)*DEXP(-Dx**E(11))*Tx**Tl(11)*Dx**D(11) $ +A(12)*DEXP(-Dx**E(12))*Tx**Tl(12)*Dx**D(12) $ +A(13)*DEXP(-Dx**E(13))*Tx**Tl(13)*Dx**D(13) C3=A(14)*DEXP(-Dx**E(14))*Tx**Tl(14)*Dx**D(14) Xr=X**R DRDHX=C1*(1.-(R+1.)*Xr)+C2*X*(2.-(R+2.)*Xr) $ +C3*X*X*(3.-(R+3.)*Xr) RETURN END C ================================================================ c ** Differential of Departure Function of the Residual Part of ** c ** the Helmholtz Free Energy to Dimensionless TT ** C ================================================================ SUBROUTINE D_R_DHTT(Dx,Tx,X,DRDHTT) DOUBLE PRECISION A(1:14) REAL tl(1:14),D(1:14),E(2:14) DOUBLE PRECISION Tx,Dx,R,DRDHTT,X DATA (A(I),I=1,14) /-1.855822D-2,5.258010D-2,3.552874D-10 $ ,5.451379D-6,-5.998546D-13,-3.687808D-6,0.2586192 $ ,-1.368072D-8,1.226146D-2,-7.181443D-2,9.970849D-2 $ ,1.0584086D-3,-0.1963687,-0.7777897/ DATA (tl(I),I=1,14) / 1.5,0.5,6.5,1.75,15.0,6.0,-1.0,4.0 $ ,3.5,0.0,-1.0,8.0,7.5,4.0 / DATA (D(I),I=1,14) / 4.0,5.0,15.0,12.0,12.0,15.0,4.0,15.0 $ ,4.0,5.0,6.0,10.0,6.0,2.0 / DATA (E(I),I=2,14) / 1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0 $ ,2.0,2.0,2.0 / R=0.5248379 DRDHTT=A(1)*Dx**D(1)*tl(1)*(tl(1)-1.)*Tx**(tl(1)-2.0) $ +A(2)*DEXP(-Dx**E(2))*Dx**D(2)*tl(2)*(tl(2)-1.)*Tx**(tl(2)-2.) $ +A(3)*DEXP(-Dx**E(3))*Dx**D(3)*tl(3)*(tl(3)-1.)*Tx**(tl(3)-2.) $ +A(4)*DEXP(-Dx**E(4))*Dx**D(4)*tl(4)*(tl(4)-1.)*Tx**(tl(4)-2.) $ +A(5)*DEXP(-Dx**E(5))*Dx**D(5)*tl(5)*(tl(5)-1.)*Tx**(tl(5)-2.) $ +A(6)*DEXP(-Dx**E(6))*Dx**D(6)*tl(6)*(tl(6)-1.)*Tx**(tl(6)-2.) DRDHTT=DRDHTT+X* $ (A(7)*DEXP(-Dx**E(7))*Dx**D(7)*tl(7)*(tl(7)-1.)*Tx**(tl(7)-2.) $ +A(8)*DEXP(-Dx**E(8))*Dx**D(8)*tl(8)*(tl(8)-1.)*Tx**(tl(8)-2.) $ +A(9)*DEXP(-Dx**E(9))*Dx**D(9)*tl(9)*(tl(9)-1.)*Tx**(tl(9)-2.) $ +A(10)*DEXP(-Dx**E(10))*Dx**D(10)*tl(10)*(tl(10)-1.)* $ Tx**(tl(10)-2.) $ +A(11)*DEXP(-Dx**E(11))*Dx**D(11)*tl(11)*(tl(11)-1.)* $ Tx**(tl(11)-2.) $ +A(12)*DEXP(-Dx**E(12))*Dx**D(12)*tl(12)*(tl(12)-1.)* $ Tx**(tl(12)-2.) $ +A(13)*DEXP(-Dx**E(13))*Dx**D(13)*tl(13)*(tl(13)-1.)* $ Tx**(tl(13)-2.)) DRDHTT=DRDHTT+X*X*A(14)* $ DEXP(-Dx**E(14))*Dx**D(14)*tl(14)*(tl(14)-1.)*Tx**(tl(14)-2.) DRDHTT=X*(1.0-X**R)*DRDHTT RETURN END C ================================================================ c ** Differential of Departure Function of the Residual Part of ** c ** the Helmholtz Free Energy to Dimensionless VT ** C ================================================================ SUBROUTINE D_R_DHTV(Dx,Tx,X,DRDHVT) DOUBLE PRECISION A(1:14),C(1:14) REAL tl(1:14),D(1:14),E(2:14) DOUBLE PRECISION Tx,Dx,R,DRDHVT,X DATA (A(I),I=1,14) /-1.855822D-2,5.258010D-2,3.552874D-10 $ ,5.451379D-6,-5.998546D-13,-3.687808D-6,0.2586192 $ ,-1.368072D-8,1.226146D-2,-7.181443D-2,9.970849D-2 $ ,1.0584086D-3,-0.1963687,-0.7777897/ DATA (tl(I),I=1,14) / 1.5,0.5,6.5,1.75,15.0,6.0,-1.0,4.0 $ ,3.5,0.0,-1.0,8.0,7.5,4.0 / DATA (D(I),I=1,14) / 4.0,5.0,15.0,12.0,12.0,15.0,4.0,15.0 $ ,4.0,5.0,6.0,10.0,6.0,2.0 / DATA (E(I),I=2,14) / 1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0 $ ,2.0,2.0,2.0 / R=0.5248379 C(1)=A(1)*D(1)*Dx**(D(1)-1.)*tl(1) I=1 DO 10 WHILE(I.LE.13) I=I+1 C(I)=A(I)*tl(I)*DEXP(-Dx**E(I))* $ (D(I)*Dx**(D(I)-1.)-Dx**D(I)*E(I)*Dx**(E(I)-1.)) 10 CONTINUE DRDHVT=C(1)*Tx**(Tl(1)-1.)+ $ C(2)*Tx**(Tl(2)-1.)+C(3)*Tx**(Tl(3)-1.)+C(4)*Tx**(Tl(4)-1.) $ +C(5)*Tx**(Tl(5)-1.)+C(6)*Tx**(Tl(6)-1.)+ $ X*(C(7)*Tx**(Tl(7)-1.)+C(8)*Tx**(Tl(8)-1.)+C(9)*Tx**(Tl(9)-1.) $ +C(10)*Tx**(Tl(10)-1.)+C(11)*Tx**(Tl(11)-1.) $ +C(12)*Tx**(Tl(12)-1.)+C(13)*Tx**(Tl(13)-1.))+ $ X*X*C(14)*Tx**(Tl(14)-1.) DRDHVT=X*(1.0-X**R)*DRDHVT RETURN END C ================================================================ c ** Differential of Departure Function of the Residual Part of ** c ** the Helmholtz Free Energy to Dimensionless VT ** C ================================================================ SUBROUTINE D_R_DHVV(Dx,Tx,X,DRDHVV) DOUBLE PRECISION A(1:14),C(1:14) REAL tl(1:14),D(1:14),E(2:14) DOUBLE PRECISION Tx,Dx,R,DRDHVV,X DATA (A(I),I=1,14) /-1.855822D-2,5.258010D-2,3.552874D-10 $ ,5.451379D-6,-5.998546D-13,-3.687808D-6,0.2586192 $ ,-1.368072D-8,1.226146D-2,-7.181443D-2,9.970849D-2 $ ,1.0584086D-3,-0.1963687,-0.7777897/ DATA (tl(I),I=1,14) / 1.5,0.5,6.5,1.75,15.0,6.0,-1.0,4.0 $ ,3.5,0.0,-1.0,8.0,7.5,4.0 / DATA (D(I),I=1,14) / 4.0,5.0,15.0,12.0,12.0,15.0,4.0,15.0 $ ,4.0,5.0,6.0,10.0,6.0,2.0 / DATA (E(I),I=2,14) / 1.0,1.0,1.0,1.0,2.0,1.0,1.0,1.0,1.0,2.0 $ ,2.0,2.0,2.0 / R=0.5248379 C(1)=A(1)*Tx**tl(1)*D(1)*(D(1)-1.)*Dx**(D(1)-2.) I=1 DO 10 WHILE(I.LE.13) I=I+1 C(I)=DEXP(-Dx**E(I))*(D(I)*(D(I)-1.)*Dx**(D(I)-2.) $ -E(I)*(2.*D(I)+E(I)-1)*Dx**(D(I)+E(I)-2.) $ +E(I)*E(I)*Dx**(D(I)+2.*E(I)-2.)) C(I)=A(I)*Tx**tl(I)*C(I) 10 CONTINUE DRDHVV=C(1)+C(2)+C(3)+C(4)+C(5)+C(6)+ $ X*(C(7)+C(8)+C(9)+C(10)+C(11)+C(12)+C(13))+X*X*C(14) DRDHVV=X*(1.0-X**R)*DRDHVV RETURN END C ======================================================== c ** Fugacity Coefficients 's Coefficient Ffa ** C ======================================================== SUBROUTINE Fu_C(Dx,Tx,X,F1,F2) DOUBLE PRECISION Tx,Dx,F1,F2,X DOUBLE PRECISION Kv,Kt,Alfa,Belt,Tc12,Vc12,Tnx,Vnx DOUBLE PRECISION M1,M2 M1=0.018015268 M2=0.01703026 Tc1=647.096 Tc2=405.40 Vc1=1.0/322.0*M1 Vc2=1.0/225.0*M2 Kv=1.2395117 Kt=0.9648407 Alfa=1.125455 Belt=0.8978069 Tc12=0.5*Kt*(Tc1+Tc2) Vc12=0.5*Kv*(Vc1+Vc2) Tnx=(1.0-X)**2.*Tc1+X**2.*Tc2+2.0*X*(1.0-X**Alfa)*Tc12 Vnx=(1.0-x)**2.*Vc1+X**2.*Vc2+2.0*X*(1.0-X**Belt)*Vc12 F1=-2.*(1.0-X)*Tc1+2*X*Tc2+2.*(1.-(1.+Alfa)*X**Alfa)*Tc12 F1=Tx/Tnx*F1 F2=-2.*(1.-X)*Vc1+2.*X*Vc2+2.*(1.-(1.+Belt)*X**Belt)*Vc12 F2=Dx/Vnx*F2 RETURN END C ============================================================ C ** The Helmholtz Free Energy of the Mixture ** C ** (Ammonia and Water) of The residual part ** C ============================================================ SUBROUTINE H_RMHG(Dx,Tx,X,HRMHG) DOUBLE PRECISION Dx,Tx,HRMHG,X DOUBLE PRECISION HRDHG,HRAHG,HRWHG IF((X.GT.0.).AND.(X.LT.1.)) THEN CALL H_RDHG(Dx,Tx,X,HRDHG) ELSE HRDHG=0 ENDIF IF(X.GT.0.) THEN CALL H_RAHG(Dx,Tx,HRAHG) ELSE HRAHG=0 ENDIF IF(X.LT.1.) THEN CALL H_RWHG(Dx,Tx,HRWHG) ELSE HRWHG=0 ENDIF HRMHG=(1.-X)*HRWHG+X*HRAHG+HRDHG RETURN END C ================================================================ c ** Differential of Departure Function of the Residual Part of ** c ** the Helmholtz Free Energy to Dimensionless Temperature ** C ================================================================ SUBROUTINE D_R_MHT(Dx,Tx,X,DRMHT) DOUBLE PRECISION Dx,Tx,DRMHT,X DOUBLE PRECISION DRDHT,DRAHT,DRWHT IF((X.LT.1.).AND.(X.GT.0.)) THEN CALL D_R_DHT(Dx,Tx,X,DRDHT) ELSE DRDHT=0 ENDIF IF(X.GT.0.) THEN CALL D_R_AHT(Dx,Tx,DRAHT) ELSE DRAHT=0 ENDIF IF(X.LT.1.) THEN CALL D_R_WHT(Dx,Tx,DRWHT) ELSE DRWHT=0. ENDIF DRMHT=(1.-X)*DRWHT+X*DRAHT+DRDHT RETURN END C =============================================================== c * Differential of Department Function of the Residual Part of * c * the Helmholtz Free Energy to Dimensionless Temperature TT * C =============================================================== SUBROUTINE D_R_MHTT(Dx,Tx,X,DRMHTT) DOUBLE PRECISION Dx,Tx,DRMHTT,X DOUBLE PRECISION DRDHTT,DRWHTT,DRAHTT IF((X.LT.1.).AND.(X.GT.0.)) THEN CALL D_R_DHTT(Dx,Tx,X,DRDHTT) ELSE DRDHTT=0 ENDIF IF(X.LT.1.) THEN CALL D_R_WHTT(Dx,Tx,DRWHTT) ELSE DRWHTT=0 ENDIF IF(X.GT.0.) THEN CALL D_R_AHTT(Dx,Tx,DRAHTT) ELSE DRAHTT=0 ENDIF DRMHTT=(1.-X)*DRWHTT+X*DRAHTT+DRDHTT RETURN END C =============================================================== c * Differential of Department Function of the Residual Part of * c * the Helmholtz Free Energy to Dimensionless Volume * C =============================================================== SUBROUTINE D_R_MHV(Dx,Tx,X,DRMHV) DOUBLE PRECISION Dx,Tx,DRMHV,X DOUBLE PRECISION DRDHV,DRWHV,DRAHV IF((X.LT.1.).AND.(X.GT.0.)) THEN CALL D_R_DHV(Dx,Tx,X,DRDHV) ELSE DRDHV=0.0 ENDIF IF(X.LT.1.) THEN CALL D_R_WHV(Dx,Tx,DRWHV) ELSE DRWHV=0 ENDIF IF(X.GT.0.) THEN CALL D_R_AHV(Dx,Tx,DRAHV) ELSE DRAHV=0 ENDIF DRMHV=(1.-X)*DRWHV+X*DRAHV+DRDHV RETURN END C ================================================================ c * Differential of Department Function of the Residual Part of * c * the Helmholtz Free Energy to Dimension VV * C ================================================================ SUBROUTINE D_R_MHVV(Dx,Tx,X,DRMHVV) DOUBLE PRECISION Dx,Tx,DRMHVV,X DOUBLE PRECISION DRWHVV,DRAHVV,DRDHVV IF((X.LT.1.).AND.(X.GT.0.)) THEN CALL D_R_DHVV(Dx,Tx,X,DRDHVV) ELSE DRDHVV=0 ENDIF IF(X.LT.1.) THEN CALL D_R_WHVV(Dx,Tx,DRWHVV) ELSE DRWHVV=0 ENDIF IF(X.GT.0) THEN CALL D_R_AHVV(Dx,Tx,DRAHVV) ELSE DRAHVV=0 ENDIF DRMHVV=(1.-X)*DRWHVV+X*DRAHVV+DRDHVV RETURN END C ================================================================ c * Differential of Departmeent Function of the Residual Part of * c * the Helmholtz Free Energy to Dimension TV * C ================================================================ SUBROUTINE D_R_MHTV(Dx,Tx,X,DRMHTV) DOUBLE PRECISION Dx,Tx,DRMHTV,X DOUBLE PRECISION DRWHTV,DRAHTV,DRDHTV IF((X.LT.1.).AND.(X.GT.0)) THEN CALL D_R_DHTV(Dx,Tx,X,DRDHTV) ELSE DRDHTV=0 ENDIF IF(X.LT.1.) THEN CALL D_R_WHTV(Dx,Tx,DRWHTV) ELSE DRWHTV=0 ENDIF IF(X.GT.0.) THEN CALL D_R_AHTV(Dx,Tx,DRAHTV) ELSE DRAHTV=0 ENDIF DRMHTV=(1.-X)*DRWHTV+X*DRAHTV+DRDHTV RETURN END C =============================================================== c * Differential of Department Function of the Residual Part of * c * the Helmholtz Free Energy to Dimension X * c =============================================================== SUBROUTINE D_R_MHX(Dx,Tx,X,DRMHX) DOUBLE PRECISION Dx,Tx,DRMHX,X DOUBLE PRECISION HRWHG,HRAHG,DRDHX CALL D_R_DHX(Dx,Tx,X,DRDHX) CALL H_RWHG(Dx,Tx,HRWHG) CALL H_RAHG(Dx,Tx,HRAHG) DRMHX=-HRWHG+HRAHG+DRDHX RETURN END C =============================================================== C * Differential of Department Function of the Residual Part of * C * the Helmholtz Free Energy to X and dimensionless T * C =============================================================== SUBROUTINE D_R_MHTX(Dx,Tx,X,DRMHTX) DOUBLE PRECISION Dx,Tx,X,DRMHTX DOUBLE PRECISION DRDHTX,DRAHT,DRWHT CALL D_R_DHTX(Dx,Tx,X,DRDHTX) CALL D_R_AHT(Dx,Tx,DRAHT) CALL D_R_WHT(Dx,Tx,DRWHT) DRMHTX=-DRWHT+DRAHT+DRMHTX RETURN END C =============================================================== C * Differential of Department Function of the Residual Part of * C * the Helmholtz Free Energy to X and Dimensionless V * C =============================================================== SUBROUTINE D_R_MHVX(Dx,Tx,X,DRMHVX) DOUBLE PRECISION Dx,Tx,X,DRMHVX DOUBLE PRECISION DRDHVX,DRAHV,DRWHV CALL D_R_DHVX(Dx,Tx,X,DRDHVX) CALL D_R_AHV(Dx,Tx,DRAHV) CALL D_R_WHV(Dx,Tx,DRWHV) DRMHVX=-DRWHV+DRAHV+DRMHVX RETURN END C ======================================================== C ** The coefficients of Fugacity Coefficients *** C ======================================================== SUBROUTINE FU_FA(Dx,Tx,X,FUFA) DOUBLE PRECISION Dx,Tx,FUFA,X DOUBLE PRECISION F1,F2 DOUBLE PRECISION DRMHX,DRMHT,DRMHV CALL D_R_MHX(Dx,Tx,X,DRMHX) CALL FU_C(Dx,Tx,X,F1,F2) CALL D_R_MHT(Dx,Tx,X,DRMHT) CALL D_R_MHV(Dx,Tx,X,DRMHV) FUFA=DRMHX+F2*DRMHV+F1*DRMHT RETURN END C ======================================================= C ** Compressibility Factor Z ** C ======================================================= SUBROUTINE COMP_Z(Dx,Tx,X,Z) DOUBLE PRECISION Dx,Tx,Z,X DOUBLE PRECISION DRMHV CALL D_R_MHV(Dx,Tx,X,DRMHV) Z=1+Dx*DRMHV RETURN END C ====================================================== C ** DIMENSIONLESS ENTHALPY H ** C ====================================================== SUBROUTINE DIM_H(Dx,Tx,Do,To,X,H) DOUBLE PRECISION Dx,Tx,Do,To,H,X DOUBLE PRECISION DRMHV,DIHT,DRMHT CALL D_R_MHV(Dx,Tx,X,DRMHV) CALL D_I_HT(To,X,DIHT) CALL D_R_MHT(Dx,Tx,X,DRMHT) H=1.+Dx*DRMHV+To*DIHT+Tx*DRMHT RETURN END C ====================================================== c ** DIMENSIONLESS ENTROPY S ** C ====================================================== SUBROUTINE DIM_S(Dx,Tx,Do,To,X,S) DOUBLE PRECISION Dx,Tx,Do,To,X,S DOUBLE PRECISION DIHT,DRMHT,HIMHG,HRMHG CALL D_I_HT(To,X,DIHT) CALL D_R_MHT(Dx,Tx,X,DRMHT) CALL H_IMHG(Do,To,X,HIMHG) CALL H_RMHG(Dx,Tx,X,HRMHG) S=To*DIHT+Tx*DRMHT-HIMHG-HRMHG RETURN END C ====================================================== C ** DIMENSIONLESS ISOCHORIC HEAT CAPACITY ** C ====================================================== SUBROUTINE DIM_CV(Dx,Tx,Do,To,X,CV) DOUBLE PRECISION Dx,Tx,Do,To,X,CV DOUBLE PRECISION DIHTT,DRMHTT CALL D_I_HTT(To,X,DIHTT) CALL D_R_MHTT(Dx,Tx,X,DRMHTT) CV=-To*To*DIHTT-Tx*Tx*DRMHTT RETURN END C ======================================================= C ** DIMENSIONLESS ISOBARIC HEAT CAPACITY ** C ======================================================= SUBROUTINE DIM_CP(Dx,Tx,Do,To,X,CP) DOUBLE PRECISION Dx,Tx,Do,To,X,CP DOUBLE PRECISION CV,DRMHV,DRMHTV,DRMHVV CALL DIM_CV(Dx,Tx,Do,To,X,CV) CALL D_R_MHV(Dx,Tx,X,DRMHV) CALL D_R_MHTV(Dx,Tx,X,DRMHTV) CALL D_R_MHVV(Dx,Tx,X,DRMHVV) CP=(1+Dx*DRMHV-Dx*Tx*DRMHTV)**2/(1.+2*Dx*DRMHV+Dx*Dx*DRMHVV) CP=CV+CP RETURN END C ======================================================= C ** FUGACITY COEFFICIENTS ** C ======================================================= SUBROUTINE FUGA_C(Dx,Tx,X,FUGAC1,FUGAC2) DOUBLE PRECISION Dx,Tx,X,FUGAC1,FUGAC2 DOUBLE PRECISION FUFA,HRMHG,DRMHV,Z CALL FU_FA(Dx,Tx,X,FUFA) CALL H_RMHG(Dx,Tx,X,HRMHG) CALL D_R_MHV(Dx,Tx,X,DRMHV) CALL COMP_Z(Dx,Tx,X,Z) FUGAC1=HRMHG+Dx*DRMHV-X*FUFA FUGAC2=HRMHG+Dx*DRMHV+(1.-X)*FUFA FUGAC1=DEXP(FUGAC1)/Z FUGAC2=DEXP(FUGAC2)/Z RETURN END C =========================================================== C ** Single Phase Properties: P,V,X ====>T ** C =========================================================== SUBROUTINE PVX_T(P,V,X,T) DOUBLE PRECISION P,V,X,T DOUBLE PRECISION TN,VN,Tx,Dx,To,Do,DRMHV,DRMHTV,Z DOUBLE PRECISION R,F,P1,F1 EPS=1.0 R=8.314510 CALL TN_X(X,TN,VN) T=4*40D6*V/R CALL TRAN_TV(T,V,X,Tx,Dx,To,Do) F=R*TN/V DO 100 WHILE(EPS.GT.1.0E-10) CALL D_R_MHV(Dx,Tx,X,DRMHV) CALL D_R_MHTV(Dx,Tx,X,DRMHTV) CALL COMP_Z(Dx,Tx,X,Z) P1=Z*R*TN/(V*Tx) F1=F*(-1./(Tx*Tx)-Dx/(Tx*Tx)*DRMHV+Dx/Tx*DRMHTV) Tx=Tx+((P-P1)/F1) EPS=DABS((P-P1)/P) 100 CONTINUE T=TN/Tx RETURN END C ===================================================== C ** SINGLE PHASE P,T,X ==> V ** C ** KL=0: LIQUID; KL=OTHERS: GAS ** C ===================================================== SUBROUTINE PTX_V(P,T,X,V,KL) DOUBLE PRECISION P,T,X,V DOUBLE PRECISION R,Tx,Dx,To,Do,Z DOUBLE PRECISION TN,VN,DRMHV,DRMHVV DOUBLE PRECISION EPS,F,P1,F1,M1,M2,M R=8.314510 EPS=1.E6 CALL TN_X(X,TN,VN) F=R*T/VN I=0 IF(X.GT.-1E-6.AND.X.LT.1E-6) THEN S=1 ELSE IF(KL.EQ.1) THEN S=1 ELSE S=18 ENDIF ENDIF IF(KL.EQ.1) THEN V=R*T/P ELSE V=1.5E-5 M1=0.018015268 M2=0.01703026 M=X*M2+(1.-X)*M1 V=1008-330*(T-273)/(350)-340*X**1.1 V=M*1.0/V ENDIF CALL TRAN_TV(T,V,X,Tx,Dx,To,Do) DO 200 WHILE(EPS.GT.1.0E-5) CALL D_R_MHV(Dx,Tx,X,DRMHV) CALL D_R_MHVV(Dx,Tx,X,DRMHVV) CALL COMP_Z(Dx,Tx,X,Z) P1=Z*Dx*R*T/VN F1=F*(1+2*Dx*DRMHV+Dx*Dx*DRMHVV) F1=DABS(F1) Dx=Dx+S*(P-P1)/F1 EPS0=DABS((P-P1)/P) EPS=EPS0 200 CONTINUE V=VN/Dx RETURN END C ======================================================= C ** SINGLE PHASE : P,H,X ===> V,T ** C ** KJ=0, LIQUID; KJ=OTHERS, GAS ** C ======================================================= SUBROUTINE PHX_VT(P,H,X,V,T,KJ) DOUBLE PRECISION P,H,X,V,T DOUBLE PRECISION TN,VN,R,Z,DIMH,Dx,Tx,Do,To DOUBLE PRECISION DRMHV,DRMHVV,DRMHTV,DRMHT,DRMHTT,DIHT,DIHTT DOUBLE PRECISION P1,H1,FP,FH DOUBLE PRECISION M1,M2,M REAL EPS1,EPS2,EPS,EPS0,REL1 M1=0.018015268 M2=0.01703026 M=X*M2+(1.-X)*M1 EPS=1.0 EPS0=1.0E8 R=8.314510 IF(P.LE.10000000) THEN IF(KJ.EQ.0) THEN V=1008-340*1.1*X V=M/V T=300.0 J=0 ELSE V=M*1.0E5/P T=500.0 J=1 ENDIF ELSE IF(KJ.EQ.0) THEN V=1008-340*1.1*X V=M/V T=350.0 J=0 ELSE V=M*1.0E5/P T=550 J=1 ENDIF ENDIF I=1 CALL TN_X(X,TN,VN) DO 100 WHILE(EPS.GT.1.0E-5) IF(V.LE.0.0) V=1.0E-5 IF(T.LE.0.0) T=0.01 Dx=VN/V Tx=TN/T To=500/T Do=1/(15000*V) CALL COMP_Z(Dx,Tx,X,Z) P1=R*TN*Dx/(Tx*VN)*Z CALL D_R_MHV(Dx,Tx,X,DRMHV) CALL D_R_MHVV(Dx,Tx,X,DRMHVV) FP=R*TN/(Tx*VN)*(1+2*Dx*DRMHV+Dx*Dx*DRMHVV) FP=DABS(FP) EPS1=DABS((P-P1)/P) EPS=MAX(EPS1,EPS2) IF((EPS.GE.EPS0).OR.(X.LT.1E-5.AND.X.GT.-1E-5)) THEN I=I+30 ELSE I=I-1 ENDIF IF(I.GT.300) I=300 IF(I.LT.-30) I=-30 IF(J.EQ.0) THEN IF(I.GE.1) THEN REL1=0.8 ELSE REL1=20 ENDIF ELSE REL1=1 ENDIF IF(X.LT.1E-5.AND.X.GT.-1E-5) THEN IF(J.EQ.0) REL1=1 IF(J.EQ.1) REL1=0.3 ENDIF Dx=Dx+REL1*(P-P1)/FP CALL DIM_H(Dx,Tx,Do,To,X,DIMH) H1=R*TN/Tx*DIMH CALL D_R_MHTV(Dx,Tx,X,DRMHTV) CALL D_R_MHT(Dx,Tx,X,DRMHT) CALL D_R_MHTT(Dx,Tx,X,DRMHTT) CALL D_I_HT(To,X,DIHT) CALL D_I_HTT(To,X,DIHTT) FH=R*TN*(-1.0/(Tx*Tx)*DIMH $ +1.0/Tx*(Dx*DRMHTV+DRMHT+Tx*DRMHTT $ +500./TN*DIHT+500*To/TN*DIHTT)) EPS2=DABS((H-H1)/H) Tx=Tx+0.1*(H-H1)/FH T=TN/Tx V=Vn/Dx EPS0=EPS 100 CONTINUE RETURN END C ========================================================== C ** SINGLE PHASE : P,S,X ==> V,T ** C ========================================================== SUBROUTINE PSX_VT(P,S,X,V,T) DOUBLE PRECISION P,S,X,V,T DOUBLE PRECISION TN,VN,R,Z,Dx,Tx,To,Do DOUBLE PRECISION FP,FS DOUBLE PRECISION S1,P1,DRMHV,DRMHVV,DIMS,DIHTT,DRMHTT REAL EPS1,EPS2,EPS,EPS0,REL1,REL2 DOUBLE PRECISION M1,M2,M,SC,SC1 M1=0.018015268 M2=0.01703026 M=X*M2+(1.-X)*M1 I=1 J=0 EPS0=1.0E8 EPS=1.0 R=8.314510 T=273.15 V=((1-X)*8E-4+X*1.15E-3)*R*T*0.000009 CALL TN_X(X,TN,VN) DO 100 WHILE(EPS.GT.1.0D-5) IF(V.LE.0.0) V=1.0E-5 IF(T.LE.0.0) T=0.01 IF(T.GT.2000) T=700 Dx=VN/V Tx=TN/T To=500/T Do=1/(15000*V) CALL DIM_S(Dx,Tx,Do,To,X,DIMS) S1=R*DIMS CALL D_I_HTT(To,X,DIHTT) CALL D_R_MHTT(Dx,Tx,X,DRMHTT) FS=R*(To*500./TN*DIHTT+Tx*DRMHTT) EPS2=DABS((S-S1)/S) REL2=0.1 Tx=Tx+REL2*(S-S1)/FS SC=S/M/1000.0 SC1=S1/M/1000.0 CALL COMP_Z(Dx,Tx,X,Z) P1=R*TN*Dx/(Tx*VN)*Z CALL D_R_MHV(Dx,Tx,X,DRMHV) CALL D_R_MHVV(Dx,Tx,X,DRMHVV) FP=R*TN/(Tx*VN)*(1+2*Dx*DRMHV+Dx*Dx*DRMHVV) FP=DABS(FP) EPS1=DABS((P-P1)/P) EPS=MAX(EPS1,EPS2) IF((EPS.GE.EPS0).OR.(X.LT.1E-5.AND.X.GT.-1E-5)) THEN I=I+20 ELSE I=I-1 J=0 ENDIF IF(I.GT.300) THEN I=300 J=J+1 ENDIF IF(X.LT.1E-5.AND.X.GT.-1E-5) J=0 IF(I.LT.-30) I=-30 IF(I.GT.1) THEN REL1=1 ELSE REL1=20 ENDIF IF(X.LT.1E-5.AND.X.GT.-1E-5) REL1=1 Dx=Dx+REL1*(P-P1)/FP EPS0=EPS T=TN/Tx V=Vn/Dx 100 CONTINUE RETURN END C ======================================================= C ** P,X ===> T,XL,XG,VL,VG ** C ** BUBBLE POINT : XL=X ** C ======================================================= SUBROUTINE BPX_TXV(P,X,T,XL,XG,VL,VG) DOUBLE PRECISION P,X,T,XL,XG,VL,VG DOUBLE PRECISION R,TNG,VNG,TNL,VNL,DRMHV,DRMHVV DOUBLE PRECISION DGx,DLx,TGx,TLx DOUBLE PRECISION FCL1,FCL2,FCG1,FCG2 DOUBLE PRECISION A,B,C,A1,T1,S1,S2 DOUBLE PRECISION HRMHGG,HRMHGL,DRMHVL,DRMHVG,DRMHXL,DRMHXG DOUBLE PRECISION DRMHTL,DRMHTG DOUBLE PRECISION DDNXL,DDNXG,DTNXL,DTNXG,ZL,ZG DOUBLE PRECISION F,F1,COEX1,COEX2,COEX4,M1,M2,M DOUBLE PRECISION DRMHVTL,DRMHVTG,DRMHTXL,DRMHTXG,DRMHTTL,DRMHTTG DOUBLE PRECISION DLx0,DGx0,XG0 REAL A0(0:10),B0(0:10),XI(0:10) M1=0.018015268 M2=0.01703026 M=X*M2+(1.-X)*M1 R=8.314510 EPS=1.0 EPS01=1.0e6 EPS02=1.0e6 EPS4=1.E-3 EPS2=10 C ************ INITIAL TEPERATURE ***************** DATA(A0(I),I=0,10)/24.61, 23.63, 23.11, 22.71, 22.35, $ 22.25, 22.08, 22.11, 22.29, 22.47, 23.06/ DATA(B0(I),I=0,10)/4932.7, 4185.7, 3744.8, 3387.5, 3072.8, $ 2869.2, 2694.0, 2611.0, 2607.1, 2623.6, 2768.9/ DATA(XI(I),I=0,10)/0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1./ I=10.*X IF(I.EQ.10) THEN A=A0(10) B=B0(10) ELSE A=A0(I)+10.*(X-XI(I))*(A0(I+1)-A0(I)) B=B0(I)+10.*(X-XI(I))*(B0(I+1)-B0(I)) ENDIF T=B/(A-LOG(P)) C ********************* END ************************* IF(T.LT.200) T=200 VL=2E-5 VG=R*T/P XL=X IF(XL.LT.1.0E-5.AND.XL.GT.-1.0E-5) THEN XG=0 ELSE IF(X.LT.(1.+1.0E-5).AND.X.GT.(1.-1.0E-5)) THEN XG=1 ELSE XG=0.99 ENDIF J=0 C ******************** TOTAL CYCLE *********************** DO 100 WHILE(EPS.GT.1.0E-5) J=J+1 IF(X.GT.-1E-5.AND.X.LT.1E-5) THEN COEX1=1 COEX2=0.8 ELSE IF(DABS(X-1).LT.1E-5) THEN COEX1=500*(1-(T-200)/250) COEX2=1 ELSE IF(T.GT.413.15) THEN IF(T.GT.500) THEN COEX1=10 ELSE COEX1=100*(1-0.1*X*X) ENDIF IF(T.GE.500) THEN COEX4=1 COEX2=1 ELSE IF(T.GE.470) THEN COEX4=5/J**0.3 COEX2=1 ELSE COEX4=10/((X*X+0.1)*J**0.3) COEX2=2 ENDIF IF(T.GE.433.15.AND.X.GE.0.9) COEX4=0.8 IF(T.GE.453.15.AND.X.GE.0.8) COEX4=0.8 IF(T.GE.473.15.AND.X.GE.0.7) COEX4=0.8 IF(T.GE.503.15.AND.X.GE.0.6) COEX4=0.5 IF(T.GE.523.15.AND.X.GE.0.5) COEX4=0.5 IF(T.GE.553.15.AND.X.GE.0.4) COEX4=0.5 IF(T.GE.573.15.AND.X.GE.0.3) COEX4=0.5 IF(T.GE.593.15.AND.X.GE.0.2) COEX4=0.5 IF(T.GE.613.15.AND.X.GE.0.8) COEX4=1 ELSE IF(T.GT.343.15) THEN COEX1=100*(1-0.1*X*X) COEX2=1*(1+X) COEX4=10/(X*X*X+0.1)/J**0.3 ELSE IF(T.GT.313.15) THEN COEX1=100 COEX2=1 COEX4=50/(X*X*X+0.1)/J**0.3 ELSE IF(T.GT.265) THEN COEX1=300 COEX2=1 COEX4=50/(X*X*X+0.01)/J**0.3 ELSE IF(T.GT.220) THEN COEX1=550 COEX2=1 COEX4=200/(X**3*J**0.3) ELSE COEX1=550 COEX2=1 COEX4=200/(X**5*(1-X)*J**0.3) ENDIF COEX3=1 ENDIF I=0 C --------- THE VOLUME OF SATURATION LIQUID ----------- DO 99 WHILE(EPS01.GT.0.8*EPS4) I=I+1 CALL TN_X(XL,TNL,VNL) DLx=VNL/VL TLx=TNL/T CALL D_R_MHV(DLx,TLx,XL,DRMHV) CALL D_R_MHVV(DLx,TLx,XL,DRMHVV) CALL COMP_Z(DLx,TLx,XL,ZL) P1=ZL*DLx*R*T/VNL F=R*T/VNL F1=F*(1.+2*DLx*DRMHV+DLx*DLx*DRMHVV) F1=DABS(F1) DLx0=DLx+COEX1*((P-P1)/F1) IF(DLx0.LT.1.) THEN T=0.99*T DLx0=3 ENDIF IF(DLx0.LE.0.0) DLx0=1.0E-5 EPS1=DABS((P-P1)/P) IF(EPS1.LT.EPS01) COEX1=0.99*COEX1 IF(EPS1.LT.(EPS01+1.E-6).AND.EPS1.GT.(EPS01-1.E-6)) $ EPS1=1.1E-5 IF(COEX1.LT.1) THEN IF(XL.GT.1.0E-5) THEN VL=3.*VNL COEX1=10 ENDIF ENDIF DLx=DLx0 VL=VNL/DLx IF(I.LT.500) THEN EPS01=EPS1 ELSE EPS01=0 ENDIF 99 CONTINUE EPS01=1.0e6 C ------------ END OF SATURATION LIQUID VOLUM ------------- I=0 C ---------- THE VOLUME OF SATURATION GAS --------- DO 999 WHILE(EPS02.GT.0.8*EPS4) I=I+1 CALL TN_X(XG,TNG,VNG) DGx=VNG/VG TGx=TNG/T CALL D_R_MHV(DGx,TGx,XG,DRMHV) CALL D_R_MHVV(DGx,TGx,XG,DRMHVV) CALL COMP_Z(DGx,TGx,XG,ZG) P1=ZG*DGx*R*T/VNG F=R*T/VNG F1=F*(1.+2*DGx*DRMHV+DGx*DGx*DRMHVV) DGx0=DGx+COEX2*(P-P1)/F1 IF(DGx0.GT.1.) DGx0=0.8 IF(DGx0.LE.0.0) DGx0=1.0E-5 EPS2=DABS((P-P1)/P) IF(EPS2.GE.EPS02) COEX2=0.95*COEX2 IF(EPS2.LT.(EPS02+1.E-6).AND.EPS2.GT.(EPS02-1.E-6)) $ EPS2=1.1E-5 DGx=DGx0 VG=VNG/DGx IF(I.LT.500) THEN EPS02=EPS2 ELSE EPS02=0 ENDIF 999 CONTINUE C -------- END OF THE SATURATION GAS VOLUM -------- EPS02=1.e6 C ----------- THE MOLE FRACTION OF GAS ------------ IF(XL.LT.1.0E-5.AND.XL.GT.-1.0E-5) THEN XG0=0 XG=XG0 ELSE IF(X.LT.(1.+1.0E-5).AND.XL.GT.(1.-1.0E-5)) THEN XG0=1.0 XG=XG0 ELSE CALL FUGA_C(DLx,TLx,XL,FCL1,FCL2) CALL FUGA_C(DGx,TGx,XG,FCG1,FCG2) IF(FCG1.EQ.0) FCG1=1.0E-9 XG0=XG+COEX3*(1.-(1.-XL)*FCL1/FCG1-XG) IF(XG0.LT.0) XG0=0.1 IF(XG0.GT.1) XG0=0.9 ENDIF EPS3=DABS(XG0-XG) XG=XG0 C --------- END OF MOLE FRACTION OF GAS ----------- C ----------- CALCULATION TEMPERATURE ------------- CALL D_R_MHT(DLx,TLx,XL,DRMHTL) CALL D_R_MHT(DGx,TGx,XG,DRMHTG) CALL H_RMHG(DLx,TLx,XL,HRMHGL) CALL H_RMHG(DGx,TGx,XG,HRMHGG) IF(DABS(X).LT.1.0E-5.OR.DABS(X-1).LT.1.0E-5) THEN CALL TN_X(XL,TNL,VNL) CALL TN_X(XG,TNG,VNG) S1=(HRMHGL-HRMHGG)+DLOG(VG/VL)-P*(VG-VL)/(R*T) S2=((TNL*DRMHTL-TNG*DRMHTG)-P*(VG-VL)/R)/(T*T) T1=T+S1/(S2*J**0.3) EPS4=100*ABS(T1-T)/T1 T=T1 ELSE CALL DTV_X(XL,DTNXL,DDNXL) CALL COMP_Z(DLx,TLx,XL,ZL) CALL D_R_MHV(DLx,TLx,XL,DRMHVL) CALL D_R_MHX(DLx,TLx,XL,DRMHXL) CALL D_R_MHTV(DLx,TLx,XL,DRMHVTL) CALL D_R_MHTX(DLx,TLx,XL,DRMHTXL) CALL D_R_MHTT(DLx,TLx,XL,DRMHTTL) CALL DTV_X(XG,DTNXG,DDNXG) CALL COMP_Z(DGx,TGx,XG,ZG) CALL D_R_MHV(DGx,TGx,XG,DRMHVG) CALL D_R_MHX(DGx,TGx,XG,DRMHXG) CALL D_R_MHTV(DGx,TGx,XG,DRMHVTG) CALL D_R_MHTX(DGx,TGx,XG,DRMHTXG) CALL D_R_MHTT(DGx,TGx,XG,DRMHTTG) A1=-DLOG(DABS(ZG/ZL))+ $ (HRMHGG-HRMHGL)+(DGx*DRMHVG-DLx*DRMHVL)+ $ (1-XG)*(DRMHXG+DDNXG/VG*DRMHVG+DTNXG/T*DRMHTG)- $ (1-XL)*(DRMHXL+DDNXL/VL*DRMHVL+DTNXL/T*DRMHTL) A=DLOG(XL/XG) B=TNG/(T*T)*(1/ZG*DGx-DRMHTG-DGx*DRMHVTG-(1-XG)* $ (DRMHTXG+DRMHTG+TGx*DRMHTTG+DGx*DRMHVTG)) C=TNL/(T*T)*(1/ZL*DLx-DRMHTL-DLx*DRMHVTL-(1-XL)* $ (DRMHTXL+DRMHTL+TLx*DRMHTTL+DLx*DRMHVTL)) T1=T+COEX4*(A-A1)/(B-C) EPS4=2*ABS(T1-T)/T1 T=T1 ENDIF IF(T.GT.2000) T=1000.0 IF(T.LT.150) T=150.0 EPS=MAX(EPS1,EPS2,EPS3,EPS4) EPS4=EPS IF(EPS4.LT.1.E-5) EPS4=1.E-5 IF(EPS4.GT.1.E-3) EPS4=1.E-3 IF(J.GT.30) EPS=0 C ------------ END OF THE TEMPERATURE ------------- 100 CONTINUE RETURN END C ======================================================= C ** P,X ===> T,XL,XG,VL,VG ** C ** DEW POINT: (XG=X) ** C ======================================================= SUBROUTINE DPX_TXV(P,X,T,XL,XG,VL,VG) DOUBLE PRECISION P,X,T,XL,XG,VL,VG DOUBLE PRECISION R,TNG,VNG,TNL,VNL,DRMHV,DRMHVV DOUBLE PRECISION DGx,DLx,TGx,TLx DOUBLE PRECISION FCL1,FCL2,FCG1,FCG2 DOUBLE PRECISION A,B,C,A1,T1,S1,S2 DOUBLE PRECISION HRMHGG,HRMHGL,DRMHVL,DRMHVG,DRMHXL,DRMHXG DOUBLE PRECISION DRMHTL,DRMHTG DOUBLE PRECISION DDNXL,DDNXG,DTNXL,DTNXG,ZL,ZG DOUBLE PRECISION F,F1,COEX1,COEX2,COEX4,M1,M2,M DOUBLE PRECISION DRMHVTL,DRMHVTG,DRMHTXL,DRMHTXG,DRMHTTL,DRMHTTG DOUBLE PRECISION DLx0,DGx0,XL0 REAL A0(0:10),B0(0:10),XI(0:10) M1=0.018015268 M2=0.01703026 M=X*M2+(1.-X)*M1 R=8.314510 EPS=1.0 EPS01=1.0e6 EPS02=1.0e6 EPS4=1.E-3 EPS2=10 C ************ INITIAL TEPERATURE ***************** DATA(A0(I),I=0,10)/24.61, 24.79, 25.00, 25.22, 25.43, $ 25.62, 25.88, 26.19, 26.55, 27.16, 23.06/ DATA(B0(I),I=0,10)/4932.7, 4996.2, 5020.6, 5045.8, $ 5059.5, 5064.3, 5077.3, 5087.4, 5083.1, 5080.7, 2768.9/ DATA(XI(I),I=0,10)/0., .1, .2, .3, .4, .5, .6, .7, .8, .9, 1./ I=10.*X IF(I.EQ.10) THEN A=A0(10) B=B0(10) ELSE IF(I.EQ.9) THEN CA1=(A0(I+1)-2*A0(I)+A0(I-1))/2 CA2=A0(I)-A0(I-1)-CA1*(XI(I)+XI(I-1)) CA3=A0(I-1)-CA1*XI(I-1)*XI(I-1)-CA2*XI(I-1) A=CA1*X*X+CA2*X+CA3 CB1=(B0(I+1)-2*B0(I)+B0(I-1))/2 CB2=B0(I)-B0(I-1)-CB1*(XI(I)+XI(I-1)) CB3=B0(I-1)-CB1*XI(I-1)*XI(I-1)-CB2*XI(I-1) B=CB1*X*X+CB2*X+CB3 ELSE A=A0(I)+10.*(X-XI(I))*(A0(I+1)-A0(I)) B=B0(I)+10.*(X-XI(I))*(B0(I+1)-B0(I)) ENDIF T=B/(A-LOG(P)) C ********************* END ************************* IF(T.LT.200) T=200 VL=2E-5 VG=R*T/P XG=X IF(X.LT.1.0E-5.AND.X.GT.-1.0E-5) THEN XL=0 ELSE IF(X.LT.(1.+1.0E-5).AND.X.GT.(1.-1.0E-5)) THEN XL=1 ELSE IF(T.GT.500) THEN XL=0.25*X ELSE XL=0.5*XG*XG ENDIF ENDIF J=0 C ******************** TOTAL CYCLE *********************** DO 100 WHILE(EPS.GT.1.0E-5) J=J+1 IF(X.GT.-1E-5.AND.X.LT.1E-5) THEN COEX1=1 COEX2=0.8 ELSE IF(DABS(X-1).LT.1E-5) THEN COEX1=500*(1-(T-200)/250) COEX2=1 ELSE IF(T.GT.585.15) THEN COEX1=10 COEX2=1 COEX3=1 COEX4=40*(0.35-X)/J**0.3 IF(X.GT.0.4) COEX1=0.05/J**0.2 ELSE IF(T.GT.553.15) THEN COEX1=100*X*X COEX2=1 COEX3=1 COEX4=40*(0.5-X)/J**0.3 IF(X.GT.0.5) COEX4=0.1/J**0.2 ELSE IF(T.GT.523.15) THEN COEX1=20*X COEX2=1 COEX3=1 COEX4=40*(0.6-X)/J**0.3 IF(X.GT.0.6) COEX4=1/J**0.2 ELSE IF(T.GT.503.15) THEN COEX1=20*X COEX2=1 COEX3=1 COEX4=40*(0.7-X)/J**0.3 IF(X.GT.0.6) COEX4=1/J**0.2 ELSE IF(T.GT.473.15) THEN COEX1=100*X*X+50*X COEX2=1 COEX3=1 COEX4=20/J**0.3 IF(X.GT.0.7) COEX4=1/J**0.2 ELSE IF(T.GT.453.15) THEN COEX1=100*X COEX2=1 COEX3=1 COEX4=100*(0.9-X)/J**0.3 IF(X.GT.0.8) COEX4=1/J**0.2 ELSE IF (T.GE.403.) THEN COEX1=100*X**0.8*(273/T) COEX2=1. COEX3=1 COEX4=(201-200*X)/J**0.3 ELSE IF(T.GT.363.15) THEN COEX1=100*X**0.8*(273/T) COEX2=1. COEX3=1 COEX4=100/(X+0.1)*(273/T)/J**0.3 ELSE IF(T.GT.323.15) THEN COEX1=100*X**0.5*(273/T) COEX2=1. COEX3=1 COEX4=200/(X*X*X+0.1)*(273/T)/J**0.3 ELSE COEX1=200*X**0.5*(273/T)**2.5 COEX2=1. COEX3=1 COEX4=1000/(X**3+0.1)*(273/T)**2.5/J**0.3 ENDIF ENDIF IF(P.GE.7.9E6.AND.X.GT.0.85.AND.X.LT.1) THEN COEX4=20*(1-X)*11.36E6/(P*J**0.1) ENDIF I=0 C --------- THE VOLUME OF SATURATION LIQUID ----------- DO 99 WHILE(EPS01.GT.0.8*EPS4) I=I+1 CALL TN_X(XL,TNL,VNL) DLx=VNL/VL TLx=TNL/T CALL D_R_MHV(DLx,TLx,XL,DRMHV) CALL D_R_MHVV(DLx,TLx,XL,DRMHVV) CALL COMP_Z(DLx,TLx,XL,ZL) P1=ZL*DLx*R*T/VNL F=R*T/VNL F1=F*(1.+2*DLx*DRMHV+DLx*DLx*DRMHVV) F1=DABS(F1) DLx0=DLx+COEX1*((P-P1)/F1) IF(DLx0.LT.1.) THEN T=0.99*T DLx0=3 ENDIF IF(DLx0.LE.0.0) DLx0=1.0E-5 EPS1=DABS((P-P1)/P) IF(EPS1.LT.EPS01) COEX1=0.99*COEX1 IF(EPS1.LT.(EPS01+1.E-6).AND.EPS1.GT.(EPS01-1.E-6)) $ EPS1=1.1E-5 IF(COEX1.LT.1) THEN IF(XL.GT.1.0E-5) THEN VL=3.*VNL COEX1=10 ENDIF ENDIF DLx=DLx0 VL=VNL/DLx IF(I.LT.500) THEN EPS01=EPS1 ELSE EPS01=0 ENDIF 99 CONTINUE EPS01=1.0e6 C ------------ END OF SATURATION LIQUID VOLUM ------------- I=0 C ---------- THE VOLUME OF SATURATION GAS --------- DO 999 WHILE(EPS02.GT.0.8*EPS4) I=I+1 CALL TN_X(XG,TNG,VNG) DGx=VNG/VG TGx=TNG/T CALL D_R_MHV(DGx,TGx,XG,DRMHV) CALL D_R_MHVV(DGx,TGx,XG,DRMHVV) CALL COMP_Z(DGx,TGx,XG,ZG) P1=ZG*DGx*R*T/VNG F=R*T/VNG F1=F*(1.+2*DGx*DRMHV+DGx*DGx*DRMHVV) DGx0=DGx+COEX2*(P-P1)/F1 IF(DGx0.GT.1.) DGx0=0.8 IF(DGx0.LE.0.0) DGx0=1.0E-5 EPS2=DABS((P-P1)/P) IF(EPS2.GE.EPS02) COEX2=0.95*COEX2 IF(EPS2.LT.(EPS02+1.E-6).AND.EPS2.GT.(EPS02-1.E-6)) $ EPS2=1.1E-5 DGx=DGx0 VG=VNG/DGx IF(I.LT.500) THEN EPS02=EPS2 ELSE EPS02=0 ENDIF 999 CONTINUE C -------- END OF THE SATURATION GAS VOLUM -------- EPS02=1.e6 C ----------- CALCULATION TEMPERATURE ------------- CALL D_R_MHT(DLx,TLx,XL,DRMHTL) CALL D_R_MHT(DGx,TGx,XG,DRMHTG) CALL H_RMHG(DLx,TLx,XL,HRMHGL) CALL H_RMHG(DGx,TGx,XG,HRMHGG) IF((XL.LT.1.0E-5.AND.XL.GT.-1.0E-5).OR. $ (XL.LT.(1.+1.E-5).AND.XL.GT.(1.-1.E-5))) THEN CALL TN_X(XL,TNL,VNL) CALL TN_X(XG,TNG,VNG) S1=(HRMHGL-HRMHGG)+DLOG(VG/VL)-P*(VG-VL)/(R*T) S2=((TNL*DRMHTL-TNG*DRMHTG)-P*(VG-VL)/R)/(T*T) T1=T+S1/S2 EPS4=ABS(T1-T)/T1 T=T+1*(T1-T) ELSE CALL DTV_X(XL,DTNXL,DDNXL) CALL COMP_Z(DLx,TLx,XL,ZL) CALL D_R_MHV(DLx,TLx,XL,DRMHVL) CALL D_R_MHX(DLx,TLx,XL,DRMHXL) CALL D_R_MHTV(DLx,TLx,XL,DRMHVTL) CALL D_R_MHTX(DLx,TLx,XL,DRMHTXL) CALL D_R_MHTT(DLx,TLx,XL,DRMHTTL) CALL DTV_X(XG,DTNXG,DDNXG) CALL COMP_Z(DGx,TGx,XG,ZG) CALL D_R_MHV(DGx,TGx,XG,DRMHVG) CALL D_R_MHX(DGx,TGx,XG,DRMHXG) CALL D_R_MHTV(DGx,TGx,XG,DRMHVTG) CALL D_R_MHTX(DGx,TGx,XG,DRMHTXG) CALL D_R_MHTT(DGx,TGx,XG,DRMHTTG) A1=DLOG((1-XL)/(1-XG))+DLOG(DABS(ZG/ZL))- $ (HRMHGG-HRMHGL)-(DGx*DRMHVG-DLx*DRMHVL)+ $ XG*(DRMHXG+DDNXG/VG*DRMHVG+DTNXG/T*DRMHTG)- $ XL*(DRMHXL+DDNXL/VL*DRMHVL+DTNXL/T*DRMHTL) B=TNG/(T*T)*(1/ZG*DGx-DRMHTG-DGx*DRMHVTG+XG* $ (DRMHTXG+DRMHTG+TGx*DRMHTTG+DGx*DRMHVTG)) C=TNL/(T*T)*(1/ZL*DLx-DRMHTL-DLx*DRMHVTL+XL* $ (DRMHTXL+DRMHTL+TLx*DRMHTTL+DLx*DRMHVTL)) A=A1/(B-C) T1=T+COEX4*A EPS4=ABS(T1-T)/T T=T1 ENDIF C ---------- END OF THE TEMPERATURE CALCULATION --------- C -------------- THE MOLE FRACTION OF GAS --------------- IF(X.LT.1.0E-5.AND.X.GT.-1.0E-5) THEN XL0=0 XL=XL0 EPS3=0 ELSE IF(X.LT.(1.+1.0E-5).AND.X.GT.(1.-1.0E-5)) THEN XL0=1.0 XL=XL0 EPS3=0 ELSE CALL FUGA_C(DLx,TLx,XL,FCL1,FCL2) CALL FUGA_C(DGx,TGx,XG,FCG1,FCG2) IF(FCL2.EQ.0) FCL2=1.0E-9 XL0=XL+COEX3*(XG*FCG2/FCL2-XL) IF(XL0.LT.0) XL0=0.01 IF(XL0.GT.1) XL0=0.99 EPS3=DABS(XL0-XL) EPS03=EPS3 XL=XL0 ENDIF EPS03=1E6 C --------- END OF MOLE FRACTION OF GAS ----------- IF(T.GT.2000) T=1000.0 IF(T.LT.150) T=150.0 EPS=MAX(EPS1,EPS2,EPS3,EPS4) EPS4=EPS IF(EPS4.LT.1.E-5) EPS4=1.E-5 IF(EPS4.GT.1.E-3) EPS4=1.E-3 IF(J.GT.35) EPS=0 C ------------ END OF THE TEMPERATURE ------------- 100 CONTINUE RETURN END C ======================================================= C ** SATUATION POINT: P,T ===> T,XL,XG,VL,VG ** C ======================================================= SUBROUTINE SPT_XV(P,T,XL,XG,VL,VG) DOUBLE PRECISION P,T,XL,XG,VL,VG,X0 DOUBLE PRECISION R,TNG,VNG,TNL,VNL,DRMHV,DRMHVV DOUBLE PRECISION DGx,DLx,TGx,TLx DOUBLE PRECISION FCL1,FCL2,FCG1,FCG2 DOUBLE PRECISION ZL,ZG DOUBLE PRECISION F,F1,COEX1,COEX2 DOUBLE PRECISION DLx0,DGx0,XG0,XL0 I=1 R=8.314510 EPS=0.1 EPS01=1.0E6 EPS02=1.0E6 IF(T.GT.415) THEN XL=0.1 XG=0.7 ELSE XL=0.3 XG=0.99 ENDIF VL=1008-330*DABS(T-273)/350 VL=0.017/VL VG=0.5*R*T/P DO 100 WHILE(EPS.GT.5.0E-5) I=I+1 IF(XL.GT.-0.5E-3.AND.XL.LT.0.5E-3) THEN COEX1=1 COEX2=1 COEX3=1 COEX4=1 ELSE IF(T.LT.300) THEN COEX1=100 COEX2=1 COEX3=0.5 COEX4=0.5 ELSE IF(T.LT.450) THEN COEX1=20 COEX2=1 COEX3=0.5 COEX4=0.5 ELSE COEX1=5 COEX2=1 COEX3=0.5 COEX4=0.5 ENDIF IF(COEX1.GT.50) COEX1=50 ENDIF J=0 C --------- THE VOLUME OF SATURATION LIQUID ----------- DO 999 WHILE(EPS01.GT.0.5*EPS) CALL TN_X(XL,TNL,VNL) DLx=VNL/VL TLx=TNL/T CALL D_R_MHV(DLx,TLx,XL,DRMHV) CALL D_R_MHVV(DLx,TLx,XL,DRMHVV) CALL COMP_Z(DLx,TLx,XL,ZL) P1=ZL*DLx*R*T/VNL F=R*T/VNL F1=F*(1.+2*DLx*DRMHV+DLx*DLx*DRMHVV) F1=DABS(F1) DLx0=DLx+COEX1*((P-P1)/F1) IF(DLx0.LT.1.) THEN DLx0=3. ENDIF EPS1=DABS((P-P1)/P) IF(EPS1.GE.EPS01+1E-5) COEX1=0.99*COEX1 IF(ABS(EPS1-EPS01).LT.1.E-6) EPS1=1.1E-5 DLx=DLx0 VL=VNL/DLx EPS01=EPS1 IF(J.GT.400) EPS01=0 999 CONTINUE C ------------ END OF THE LIQUIT VOLUM ------------- EPS01=1E6 J=0 C ---------- THE VOLUME OF SATURATION GAS --------- DO 9999 WHILE(EPS02.GT.0.5*EPS) J=J+1 CALL TN_X(XG,TNG,VNG) DGx=VNG/VG TGx=TNG/T CALL D_R_MHV(DGx,TGx,XG,DRMHV) CALL D_R_MHVV(DGx,TGx,XG,DRMHVV) CALL COMP_Z(DGx,TGx,XG,ZG) P1=ZG*DGx*R*T/VNG F=R*T/VNG F1=F*(1.+2*DGx*DRMHV+DGx*DGx*DRMHVV) F1=DABS(F1) DGx0=DGx+COEX2*(P-P1)/F1 IF(DGx0.GT.1) THEN DGx0=0.8 ENDIF IF(DGx0.LT.0) DGx0=1.0E-5 EPS2=DABS((P-P1)/P) IF(ABS(EPS2-EPS02).LT.1.0E-6) EPS2=1.0E-5 IF(EPS2.GE.EPS02) COEX2=0.99*COEX2 IF(COEX2.LT.0.1) COEX2=0.1 EPS02=EPS2 IF(J.GT.400) EPS02=0 DGx=DGx0 VG=VNG/DGx 9999 CONTINUE C -------- END OF THE SATURATION GAS VOLUM -------- EPS02=1.0E6 C ----------- THE MOLE FRACTION OF GAS AND LIQUID ------------ CALL FUGA_C(DLx,TLx,XL,FCL1,FCL2) CALL FUGA_C(DGx,TGx,XG,FCG1,FCG2) IF(FCG1.EQ.0) FCG1=1.0E-10 IF(FCL2.EQ.0) FCL2=1.0E-10 XG0=XG+COEX3*((1-(1-XL)*FCL1/FCG1)-XG) XL0=XL+COEX4*(XG*FCG2/FCL2-XL) IF(T.LT.410) THEN X0=1 ELSE X0=0.95 ENDIF IF(XG0.GT.1) XG0=0.99*X0 IF(XG0.LT.0) XG0=0.005*XG IF(XL0.GT.1) XL0=0.99*X0 IF(XL0.LT.0) XL0=0.005*XL EPS3=10*DABS(XG0-XG) EPS4=10*DABS(XL0-XL) IF(XL0.GT.XG0) XL0=XG0*XG0 C IF(XL0.LT.1.0E-8) XL0=0.8*XG0 IF(XL0.LT.1.0E-10.AND.XG0.GT.0.1) XG0=0.1 IF(XL0.GT.(1-1.0E-10).AND.XG0.LT.0.9) XG0=1 IF(XG0.LT.1.0E-10.AND.XL0.GT.0.1) XL0=0.1 EPS34=MAX(EPS3,EPS4) C --------- END OF MOLE FRACTION OF GAS ----------- XG=XG0 XL=XL0 VL=VNL/DLx VG=VNG/DGx EPS=MAX(EPS1,EPS2,EPS34) IF(EPS.GT.0.001) EPS=0.001 IF(I.GT.30) EPS=0 100 CONTINUE RETURN END C ======================================================== C ** T,X ===> P,XL,XG,VL,VG ** C ** BUBBLE POINT (XL=X) ** C ======================================================== SUBROUTINE BTX_PXV(T,X,P,XL,XG,VL,VG) DOUBLE PRECISION T,X,P,XL,XG,VL,VG DOUBLE PRECISION R,TNG,VNG,TNL,VNL,DRMHV,DRMHVV DOUBLE PRECISION DGx,DLx,TGx,TLx DOUBLE PRECISION FCL1,FCL2,FCG1,FCG2 DOUBLE PRECISION A,B DOUBLE PRECISION HRMHGG,HRMHGL DOUBLE PRECISION ZL,ZG DOUBLE PRECISION F,F1,COEX1,COEX2,COEX4,M1,M2,M,P0 DOUBLE PRECISION XG0,DTNXL,DTNXG,DDNXL,DDNXG,DRMHXL,DRMHXG, $ DRMHTL,DRMHTG,DRMHVL,DRMHVG REAL A0(0:10),B0(0:10),XI(0:10) M1=0.018015268 M2=0.01703026 M=X*M2+(1.-X)*M1 R=8.314510 EPS=1.0E6 EPS01=1.0e6 EPS02=1.0e6 EPS4=1.E-2 EPS2=10 C ************ INITIAL PRESSURE ***************** DATA(A0(I),I=0,10)/24.61, 23.63, 23.11, 22.71, 22.35, $ 22.25, 22.08, 22.11, 22.29, 22.47, 23.06/ DATA(B0(I),I=0,10)/4932.7, 4185.7, 3744.8, 3387.5, 3072.8, $ 2869.2, 2694.0, 2611.0, 2607.1, 2623.6, 2768.9/ DATA(XI(I),I=0,10)/0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1./ I=10.*X IF(I.EQ.10) THEN A=A0(10) B=B0(10) ELSE A=A0(I)+10.*(X-XI(I))*(A0(I+1)-A0(I)) B=B0(I)+10.*(X-XI(I))*(B0(I+1)-B0(I)) ENDIF P0=DEXP(A-B/T) P=P0 C ******************* END *********************** VL=1008-330*(T-273)/350-300*XL VL=M/VL VG=R*T/P XL=X IF(DABS(XL).LT.1.0E-5) THEN XG=0 ELSE IF(DABS(1.-X).LT.1.0E-5) THEN XG=1 ELSE IF(T.LT.460) XG=0.99 IF(T.GE.460.AND.T.LT.530) XG=0.9 IF(T.GE.530.AND.T.LT.573) XG=0.6 IF(T.GE.573) XG=0.4 ENDIF J=1 C ****************** TOTAL CYCLE ********************** DO 100 WHILE(EPS.GT.1.0E-5) J=J+1 IF(DABS(XL).LT.1E-5) THEN COEX1=1 COEX2=0.8 ELSE IF(DABS(XL-1.).LT.1E-5) THEN IF(T.LT.200) THEN COEX1=100 ELSE COEX1=25*(400/T)**2 ENDIF COEX2=1 ELSE IF(T.LT.273.15) COEX1=30*(400/T)**2 IF(T.GE.273.15.AND.T.LT.420) COEX1=50 IF(T.GE.420.AND.T.LT.460) COEX1=100 IF(T.GE.460.AND.T.LT.503.15) COEX1=40 IF(T.GE.503.15) COEX1=40 COEX2=1 COEX3=1 IF(T.LE.473.15) COEX4=0.8 IF(T.GT.473.15.AND.T.LT.513.15) COEX4=0.4 IF(T.GE.513.15) COEX4=0.4*513/(T*J**0.3) IF(X.LT.0.1.OR.X.GT.0.9) THEN COEX1=400*X*(1-X) COEX4=2*X*(1-X) ENDIF ENDIF I=0 C --------- THE VOLUME OF SATURATION LIQUID ----------- CALL TN_X(XL,TNL,VNL) F=R*T/VNL TLx=TNL/T DO 99 WHILE(EPS01.GT.0.8*EPS4) I=I+1 DLx=VNL/VL CALL D_R_MHV(DLx,TLx,XL,DRMHV) CALL D_R_MHVV(DLx,TLx,XL,DRMHVV) CALL COMP_Z(DLx,TLx,XL,ZL) P1=ZL*DLx*F F1=F*(1.+2*DLx*DRMHV+DLx*DLx*DRMHVV) F1=DABS(F1) DLx=DLx+COEX1*((P-P1)/F1) IF(T.GT.400) THEN IF(DLx.LT.1) DLx=3 ENDIF IF(DLx.LE.0.0) DLx=1.0E-5 EPS1=DABS((P-P1)/P) IF(EPS1.LT.EPS01) COEX1=0.99*COEX1 IF(EPS1.LT.(EPS01+1.E-10).AND.EPS1.GT.(EPS01-1.E-10)) $ EPS1=1.0E-6 VL=VNL/DLx IF(I.GT.500) THEN EPS01=0 ELSE EPS01=EPS1 ENDIF 99 CONTINUE EPS01=1.0e6 I=0 C ------------ END OF SATURATION LIQUID VOLUM ------------- C -------------- THE VOLUME OF SATURATION GAS ------------- CALL TN_X(XG,TNG,VNG) F=R*T/VNG TGx=TNG/T DO 999 WHILE(EPS02.GT.0.8*EPS4) I=I+1 DGx=VNG/VG CALL D_R_MHV(DGx,TGx,XG,DRMHV) CALL D_R_MHVV(DGx,TGx,XG,DRMHVV) CALL COMP_Z(DGx,TGx,XG,ZG) P1=ZG*DGx*F F1=F*(1.+2*DGx*DRMHV+DGx*DGx*DRMHVV) F1=DABS(F1) DGx=DGx+COEX2*(P-P1)/F1 IF(DGx.GT.1.) THEN DGx=0.8 ENDIF IF(DGx.LE.0.0) DGx=1.0E-5 EPS2=DABS((P-P1)/P) IF(EPS2.GE.EPS02) COEX2=0.99*COEX2 IF(EPS2.LT.(EPS02+1.E-10).AND.EPS2.GT.(EPS02-1.E-10)) $ EPS2=1.0E-6 VG=VNG/DGx IF(I.GT.500) THEN EPS02=0 ELSE EPS02=EPS2 ENDIF 999 CONTINUE C -------- END OF THE SATURATION GAS VOLUM -------- EPS02=1.e6 C ----------- THE MOLE FRACTION OF GAS ------------ IF(DABS(X).LT.1.0E-5) THEN XG0=0 XG=0 EPS3=0 ELSE IF(DABS(1.-X).LT.1.0E-5) THEN XG0=1 XG=1. EPS3=0 ELSE CALL FUGA_C(DLx,TLx,XL,FCL1,FCL2) CALL FUGA_C(DGx,TGx,XG,FCG1,FCG2) IF(FCG1.EQ.0) FCG1=1.0E-9 XG0=XG+COEX3*(1.-(1.-XL)*FCL1/FCG1-XG) IF(XG0.LT.0) XG0=0.01 IF(XG0.GT.1) XG0=0.999 EPS3=DABS(XG0-XG) XG=XG0 ENDIF C --------- END OF MOLE FRACTION OF GAS ----------- C ----------- PRESSURE CALCULATION ------------- CALL H_RMHG(DLx,TLx,XL,HRMHGL) CALL H_RMHG(DGx,TGx,XG,HRMHGG) IF(DABS(X).LT.1.0E-5.OR.DABS(1-X).LT.1.0E-5) THEN CALL TN_X(XL,TNL,VNL) CALL TN_X(XG,TNG,VNG) P1=R*T*(HRMHGL-HRMHGG+DLOG(VG/VL))/(VG-VL) EPS4=ABS(P1-P)/P1 P=P1 ELSE CALL DTV_X(XG,DTNXG,DDNXG) CALL DTV_X(XL,DTNXL,DDNXL) CALL TN_X(XG,TNG,VNG) TGx=TNG/T DGx=VNG/VG CALL COMP_Z(DLx,TLx,XL,ZL) CALL COMP_Z(DGx,TGx,XG,ZG) CALL D_R_MHX(DLx,TLx,XL,DRMHXL) CALL D_R_MHX(DGx,TGx,XG,DRMHXG) CALL D_R_MHT(DLx,TLx,XL,DRMHTL) CALL D_R_MHT(DGx,TGx,XG,DRMHTG) CALL D_R_MHV(DLx,TLx,XL,DRMHVL) CALL D_R_MHV(DGx,TGx,XG,DRMHVG) F=DLOG(XL/XG)+DLOG(DABS(ZG/ZL))+HRMHGL-HRMHGG $ +(1-XL)*(DRMHXL+DTNXL/T*DRMHTL+DDNXL/VL*DRMHVL) $ -(1-XG)*(DRMHXG+DTNXG/T*DRMHTG+DDNXG/VG*DRMHVG) P1=R*T*F/(VG-VL) IF(P1.LT.0) P1=P0 EPS4=ABS(P1-P)/P P=P+COEX4*(P1-P) ENDIF C ------------ END OF THE PRESSURE ------------- IF(J.GT.30) THEN EPS=0 ELSE EPS=MAX(EPS1,EPS2,EPS3,EPS4) ENDIF EPS4=EPS IF(EPS4.LT.1.E-5) EPS4=1.E-5 IF(EPS4.GT.1.E-3) EPS4=1.E-3 100 CONTINUE RETURN END C ========================================================== C ** T,X ===> P,XL,XG,VL,VG ** C ** DEW POINT : (XL=X) ** C ========================================================== SUBROUTINE DTX_PXV(T,X,P,XL,XG,VL,VG) DOUBLE PRECISION T,X,P,XL,XG,VL,VG DOUBLE PRECISION R,TNG,VNG,TNL,VNL,DRMHV,DRMHVV DOUBLE PRECISION DGx,DLx,TGx,TLx DOUBLE PRECISION FCL1,FCL2,FCG1,FCG2 DOUBLE PRECISION A,B DOUBLE PRECISION HRMHGG,HRMHGL DOUBLE PRECISION ZL,ZG DOUBLE PRECISION F,F1,COEX1,COEX2,COEX4,M1,M2,M,P0 DOUBLE PRECISION XL0,DTNXL,DTNXG,DDNXL,DDNXG,DRMHXL,DRMHXG, $ DRMHTL,DRMHTG,DRMHVL,DRMHVG REAL A0(0:10),B0(0:10),XI(0:10) M1=0.018015268 M2=0.01703026 M=X*M2+(1.-X)*M1 R=8.314510 EPS=1.0E6 EPS01=1.0e6 EPS02=1.0e6 EPS4=1.E-2 EPS2=10 C ************** INITIAL PRESSURE ***************** DATA(A0(I),I=0,10)/24.61, 24.79, 25.00, 25.22, 25.43, $ 25.62, 25.88, 26.19, 26.55, 27.16, 23.06/ DATA(B0(I),I=0,10)/4932.7, 4996.2, 5020.6, 5045.8, $ 5059.5, 5064.3, 5077.3, 5087.4, 5083.1, 5080.7, 2768.9/ DATA(XI(I),I=0,10)/0., .1, .2, .3, .4, .5, .6, .7, .8, .9, 1./ I=10.*X IF(I.EQ.10) THEN A=A0(10) B=B0(10) ELSE IF(I.EQ.9) THEN A=A0(I)+10*(X-XI(I))*(A0(I)-A0(I-1)) B=B0(I)+10*(X-XI(I))*(B0(I)-B0(I-1)) ELSE A=A0(I)+10.*(X-XI(I))*(A0(I+1)-A0(I)) B=B0(I)+10.*(X-XI(I))*(B0(I+1)-B0(I)) ENDIF P0=DEXP(A-B/T) P=P0 C ************* END OF PRESSURE INITIALIZE *************** VL=1.5E-5 VG=R*T/P XG=X IF(DABS(X).LT.1.0E-5) THEN XL=0 ELSE IF(DABS(1.-X).LT.1.0E-5) THEN XL=1 ELSE XL=0.5*X*X ENDIF J=0 C ******************* TOTAL CYCLE ********************* DO 100 WHILE(EPS.GT.1.0E-5) J=J+1 IF(X.LT.1E-5) THEN COEX1=1 COEX2=0.8 ELSE IF(DABS(XG-1.).LT.1E-5) THEN IF(T.LT.200) THEN COEX1=100 ELSE COEX1=25*(400/T)**2 ENDIF COEX2=1 ELSE IF(J.EQ.1) THEN COEX1=400*X*X*(1-X)*273/T ELSE COEX1=100*X*(2-X)*273/T ENDIF COEX2=1 COEX3=4*X*(1-X)/J**0.3 COEX4=10*X*(1-X)/J**0.3 ENDIF I=0 C --------- THE VOLUME OF SATURATION LIQUID ----------- CALL TN_X(XL,TNL,VNL) F=R*T/VNL TLx=TNL/T DO 99 WHILE(EPS01.GT.0.8*EPS4) I=I+1 DLx=VNL/VL CALL D_R_MHV(DLx,TLx,XL,DRMHV) CALL D_R_MHVV(DLx,TLx,XL,DRMHVV) CALL COMP_Z(DLx,TLx,XL,ZL) P1=ZL*DLx*F F1=F*(1.+2*DLx*DRMHV+DLx*DLx*DRMHVV) F1=DABS(F1) DLx=DLx+COEX1*((P-P1)/F1) IF(T.GT.400) THEN IF(DLx.LT.1) DLx=3 ENDIF IF(DLx.LE.0.0) DLx=1.0E-5 EPS1=DABS((P-P1)/P) IF(EPS1.LT.EPS01) THEN COEX1=0.99*COEX1 ELSE IF(EPS1.LT.(EPS01+1.E-10).AND.EPS1.GT. $ (EPS01-1.E-10)) EPS1=1.0E-6 ENDIF IF(COEX1.LT.1) THEN IF(XL.GT.1.0E-5) THEN VL=3.*VNL COEX1=10 J=5 ENDIF ENDIF VL=VNL/DLx IF(I.GT.500) THEN EPS01=0 ELSE EPS01=EPS1 ENDIF 99 CONTINUE EPS01=1.0e6 I=0 C ------------ END OF SATURATION LIQUID VOLUM ------------- C ---------- THE VOLUME OF SATURATION GAS --------- CALL TN_X(XG,TNG,VNG) F=R*T/VNG TGx=TNG/T DO 999 WHILE(EPS02.GT.0.8*EPS4) I=I+1 DGx=VNG/VG CALL D_R_MHV(DGx,TGx,XG,DRMHV) CALL D_R_MHVV(DGx,TGx,XG,DRMHVV) CALL COMP_Z(DGx,TGx,XG,ZG) P1=ZG*DGx*F F1=F*(1.+2*DGx*DRMHV+DGx*DGx*DRMHVV) F1=DABS(F1) DGx=DGx+COEX2*(P-P1)/F1 IF(DGx.GT.1.) THEN DGx=0.8 ENDIF IF(DGx.LE.0.0) DGx=1.0E-5 EPS2=DABS((P-P1)/P) IF(EPS2.GE.EPS02) THEN COEX2=0.90*COEX2 ELSE IF(EPS2.LT.(EPS02+1.E-10).AND.EPS2.GT. $ (EPS02-1.E-10)) EPS2=1.0E-6 ENDIF VG=VNG/DGx IF(I.GT.500) THEN EPS02=0 ELSE EPS02=EPS2 ENDIF 999 CONTINUE C -------- END OF THE SATURATION GAS VOLUM -------- EPS02=1.e6 C ----------- PRESSURE CALCULATION ------------- CALL H_RMHG(DLx,TLx,XL,HRMHGL) CALL H_RMHG(DGx,TGx,XG,HRMHGG) IF((XL.LT.1.0E-5.AND.XL.GT.-1.0E-5).OR. $ (XL.LT.(1.+1.E-5).AND.XL.GT.(1.-1.E-5))) THEN CALL TN_X(XL,TNL,VNL) CALL TN_X(XG,TNG,VNG) P1=R*T*(HRMHGL-HRMHGG+DLOG(VG/VL))/(VG-VL) EPS4=ABS(P1-P)/P1 P=P1 ELSE CALL COMP_Z(DLx,TLx,XL,ZL) CALL COMP_Z(DGx,TGx,XG,ZG) CALL DTV_X(XG,DTNXG,DDNXG) CALL DTV_X(XL,DTNXL,DDNXL) CALL D_R_MHX(DLx,TLx,XL,DRMHXL) CALL D_R_MHX(DGx,TGx,XG,DRMHXG) CALL D_R_MHT(DLx,TLx,XL,DRMHTL) CALL D_R_MHT(DGx,TGx,XG,DRMHTG) CALL D_R_MHV(DLx,TLx,XL,DRMHVL) CALL D_R_MHV(DGx,TGx,XG,DRMHVG) F=DLOG((1.-XL)/(1.-XG))+ $ DLOG(DABS(ZG/ZL))+HRMHGL-HRMHGG $ -XL*(DRMHXL+DTNXL/T*DRMHTL+DDNXL/VL*DRMHVL) $ +XG*(DRMHXG+DTNXG/T*DRMHTG+DDNXG/VG*DRMHVG) P1=R*T*F/(VG-VL) IF(P1.LT.0) P1=P0 P=P+COEX4*(P1-P) EPS4=ABS(P1-P)/P ENDIF C ------------ END OF THE PRESSURE ------------- C ----------- THE MOLE FRACTION OF SATURATION LIQUID ------------ IF(DABS(X).LT.1.0E-5) THEN XL=0 EPS3=0 ELSE IF(DABS(1-X).LT.1.0E-5) THEN XL=1. EPS3=0 ELSE CALL FUGA_C(DGx,TGx,XG,FCG1,FCG2) CALL FUGA_C(DLx,TLx,XL,FCL1,FCL2) IF(FCL2.EQ.0) FCL2=1.0E-9 XL0=XL+COEX3*(XG*FCG2/FCL2-XL) IF(XL0.LT.0) XL0=0.01 IF(XL0.GT.1) XL0=0.99 EPS3=DABS(XL0-XL) XL=XL0 ENDIF C --------- END OF MOLE FRACTION OF SATURATION LIQUID ----------- EPS=MAX(EPS1,EPS2,EPS3,EPS4) EPS4=EPS IF(EPS4.LT.1.E-5) EPS4=1.E-5 IF(EPS4.GT.1.E-3) EPS4=1.E-3 IF(J.GT.30) EPS4=0 100 CONTINUE RETURN END C ========================================================== C ** TWO PHASE REGION: P,H,X ===> T,XL,XG,VL,VG ** C ========================================================== SUBROUTINE PHX_TVX(P,H,X,T,XL,XG,Y,VL,VG,TB,TD) DOUBLE PRECISION P,H,X,T,XL,XG,Y,VL,VG,X0,TB,TD DOUBLE PRECISION R,TNG,VNG,TNL,VNL,DRMHV,DRMHVV DOUBLE PRECISION DGx,DLx,TGx,TLx DOUBLE PRECISION FCL1,FCL2,FCG1,FCG2 DOUBLE PRECISION ZL,ZG,Y0 DOUBLE PRECISION F,F1,COEX1,COEX2 DOUBLE PRECISION DLx0,DGx0,XG0,XL0 DOUBLE PRECISION TO,DLO,DGO,H1,FHL,FHG,FH DOUBLE PRECISION DIMHL,DIMHG,DRMHTVL,DRMHTVG,DRMHTL,DRMHTG DOUBLE PRECISION DRMHTTL,DRMHTTG,DIHTL,DIHTG,DIHTTL,DIHTTG I=1 R=8.314510 EPS=0.1 EPS01=1.0E6 EPS02=1.0E6 DO 100 WHILE(EPS.GT.1.0E-5) I=I+1 IF(ABS(X).LT.1E-5) THEN COEX1=1 COEX2=1 COEX3=1 COEX4=1 ELSE IF(T.LT.300) THEN COEX1=200 COEX2=1 COEX3=0.5 COEX4=0.5 ELSE IF(T.LT.450) THEN COEX1=100 COEX2=1 COEX3=1 COEX4=1 ELSE COEX1=50 COEX2=1 COEX3=1 COEX4=1 ENDIF IF(COEX1.GT.50) COEX1=50 ENDIF IF(X.LT.1.AND.X.GT.0) THEN COEX5=.1/I**0.3 ELSE COEX5=1 ENDIF J=0 C --------- THE VOLUME OF SATURATION LIQUID ----------- DO 999 WHILE(EPS01.GT.0.5*EPS) CALL TN_X(XL,TNL,VNL) DLx=VNL/VL TLx=TNL/T CALL D_R_MHV(DLx,TLx,XL,DRMHV) CALL D_R_MHVV(DLx,TLx,XL,DRMHVV) CALL COMP_Z(DLx,TLx,XL,ZL) P1=ZL*DLx*R*T/VNL F=R*T/VNL F1=F*(1.+2*DLx*DRMHV+DLx*DLx*DRMHVV) F1=DABS(F1) DLx0=DLx+COEX1*((P-P1)/F1) IF(DLx0.LT.1.) THEN DLx0=3. ENDIF EPS1=DABS((P-P1)/P) IF(EPS1.GE.EPS01+1E-5) COEX1=0.99*COEX1 IF(ABS(EPS1-EPS01).LT.1.E-6) EPS1=1.1E-5 DLx=DLx0 VL=VNL/DLx EPS01=EPS1 IF(J.GT.400) EPS01=0 999 CONTINUE C ------------ END OF THE LIQUIT VOLUM ------------- EPS01=1E6 J=0 C ---------- THE VOLUME OF SATURATION GAS --------- DO 9999 WHILE(EPS02.GT.0.5*EPS) J=J+1 CALL TN_X(XG,TNG,VNG) DGx=VNG/VG TGx=TNG/T CALL D_R_MHV(DGx,TGx,XG,DRMHV) CALL D_R_MHVV(DGx,TGx,XG,DRMHVV) CALL COMP_Z(DGx,TGx,XG,ZG) P1=ZG*DGx*R*T/VNG F=R*T/VNG F1=F*(1.+2*DGx*DRMHV+DGx*DGx*DRMHVV) F1=DABS(F1) DGx0=DGx+COEX2*(P-P1)/F1 IF(DGx0.GT.1) THEN DGx0=0.8 ENDIF IF(DGx0.LT.0) DGx0=1.0E-5 EPS2=DABS((P-P1)/P) IF(ABS(EPS2-EPS02).LT.1.0E-6) EPS2=1.0E-5 IF(EPS2.GE.EPS02) COEX2=0.99*COEX2 IF(COEX2.LT.0.1) COEX2=0.1 EPS02=EPS2 IF(J.GT.400) EPS02=0 DGx=DGx0 VG=VNG/DGx 9999 CONTINUE C -------- END OF THE SATURATION GAS VOLUM -------- EPS02=1.0E6 C ----------- THE MOLE FRACTION OF GAS AND LIQUID ------------ IF(ABS(X).LT.1.0E-5.OR.ABS(1-X).LT.1.0E-5) THEN TO=500/T DLO=1/(15000*VL) DGO=1/(15000*VG) CALL DIM_H(DLx,TLx,DLO,TO,XL,DIMHL) CALL DIM_H(DGx,TGx,DGO,TO,XG,DIMHG) Y0=(H-R*T*DIMHL)/(R*T*(DIMHG-DIMHL)) IF(Y0.GT.1) Y0=1 IF(Y0.LT.0) Y0=0 ELSE CALL FUGA_C(DLx,TLx,XL,FCL1,FCL2) CALL FUGA_C(DGx,TGx,XG,FCG1,FCG2) IF(FCG1.EQ.0) FCG1=1.0E-10 IF(FCL2.EQ.0) FCL2=1.0E-10 XG0=XG+COEX3*((1-(1-XL)*FCL1/FCG1)-XG) XL0=XL+COEX4*(XG*FCG2/FCL2-XL) IF(T.LT.410) THEN X0=1 ELSE X0=0.95 ENDIF IF(XG0.GT.1) XG0=0.99*X0 IF(XG0.LT.0) XG0=0.005*XG IF(XL0.GT.1) XL0=0.99*X0 IF(XL0.LT.0) XL0=0.005*XL IF(XL0.GT.X) XL0=X IF(XG0.LT.X) XG0=X IF(XL0.GT.XG0) XL0=XG0*XG0 IF(XL0.LT.1.0E-8) XL0=1.0E-8 IF(XL0.LT.1.0E-10.AND.XG0.GT.0.1) XG0=0.1 IF(XL0.GT.(1-1.0E-10).AND.XG0.LT.0.9) XG0=1 IF(XG0.LT.1.0E-10.AND.XL0.GT.0.1) XL0=0.1 XG=XG0 XL=XL0 Y0=(X-XL)/(XG-XL) ENDIF EPS3=DABS(Y-Y0) Y=Y0 C ----------- END OF MOLE FRACTION OF GAS -------------- C --------- THE CALCULATION OF TEMPERATURE ------------- TO=500/T DLO=1/(15000*VL) DGO=1/(15000*VG) CALL DIM_H(DLx,TLx,DLO,TO,XL,DIMHL) CALL DIM_H(DGx,TGx,DGO,TO,XG,DIMHG) CALL D_R_MHTV(DLx,TLx,XL,DRMHTVL) CALL D_R_MHTV(DGx,TGx,XG,DRMHTVG) CALL D_R_MHT(DLx,TLx,XL,DRMHTL) CALL D_R_MHT(DGx,TGx,XG,DRMHTG) CALL D_R_MHTT(DLx,TLx,XL,DRMHTTL) CALL D_R_MHTT(DGx,TGx,XG,DRMHTTG) CALL D_I_HT(TO,XL,DIHTL) CALL D_I_HT(TO,XG,DIHTG) CALL D_I_HTT(TO,XL,DIHTTL) CALL D_I_HTT(TO,XG,DIHTTG) H1=R*((1-Y)*TNL/TLx*DIMHL+Y*TNG/TGx*DIMHG) FHL=-R*TNL*TNL/(T*T)*(-1.0/(TLx*TLx)*DIMHL $ +1.0/TLx*(DLx*DRMHTVL+DRMHTL+TLx*DRMHTTL $ +500./TNL*DIHTL+500*To/TNL*DIHTTL)) FHG=-R*TNG*TNG/(T*T)*(-1.0/(TGx*TGx)*DIMHG $ +1.0/TGx*(DGx*DRMHTVG+DRMHTG+TGx*DRMHTTG $ +500./TNG*DIHTG+500*To/TNG*DIHTTG)) FH=(1-Y)*FHL+Y*FHG T=T+COEX5*(H-H1)/FH IF(T.GT.TD) T=TD IF(T.LT.TB) T=TB EPS4=ABS(H-H1)/H C ------- END OF THE CALCULATION OF TEMPERATURE --------- EPS=MAX(EPS1,EPS2,EPS3,EPS4) IF(EPS.GT.0.001) EPS=0.001 IF(I.GT.30) EPS=0 100 CONTINUE RETURN END C ========================================================== C ** TWO PHASE REGION: P,S,X ===> T,XL,XG,VL,VG ** C ========================================================== SUBROUTINE PSX_TVX(P,S,X,T,XL,XG,Y,VL,VG,TB,TD) DOUBLE PRECISION P,S,X,T,XL,XG,Y,VL,VG,X0,TB,TD DOUBLE PRECISION R,TNG,VNG,TNL,VNL,DRMHV,DRMHVV DOUBLE PRECISION DGx,DLx,TGx,TLx DOUBLE PRECISION FCL1,FCL2,FCG1,FCG2 DOUBLE PRECISION ZL,ZG,Y0 DOUBLE PRECISION F,F1,COEX1,COEX2 DOUBLE PRECISION DLx0,DGx0,XG0,XL0 DOUBLE PRECISION TO,DLO,DGO,S1,FSL,FSG,FS DOUBLE PRECISION DIHTTL,DIHTTG,DRMHTTL,DRMHTTG,DIMSL,DIMSG I=1 R=8.314510 EPS=1.0e-3 EPS01=1.0E6 EPS02=1.0E6 C ----------------- TOTAL CYCLE ------------------ DO 100 WHILE(EPS.GT.1.0E-4) I=I+1 IF(ABS(X).LT.1E-5) THEN COEX1=1 COEX2=1 COEX3=1 COEX4=1 ELSE IF(T.LT.300) THEN COEX1=200 COEX2=1 COEX3=(X+0.05)*(1+X)/I**0.3 COEX4=(X+0.05)*(1+X)/I**0.3 ELSE IF(T.LT.500) THEN COEX1=20 COEX2=1 COEX3=(X+0.1)*(1+X)/I**0.3 COEX4=(X+0.1)*(1+X)/I**0.3 ELSE COEX1=5 COEX2=1 COEX3=0.5*(1+X)/I**0.3 COEX4=0.5*(1+X)/I**0.3 ENDIF ENDIF IF(X.LT.1.AND.X.GT.0) THEN COEX5=(X+0.1)*(1.1-X)/I**0.3 ELSE COEX5=1 ENDIF J=0 C --------- THE VOLUME OF SATURATION LIQUID ----------- DO 999 WHILE(EPS01.GT.0.8*EPS) CALL TN_X(XL,TNL,VNL) DLx=VNL/VL TLx=TNL/T CALL D_R_MHV(DLx,TLx,XL,DRMHV) CALL D_R_MHVV(DLx,TLx,XL,DRMHVV) CALL COMP_Z(DLx,TLx,XL,ZL) P1=ZL*DLx*R*T/VNL F=R*T/VNL F1=F*(1.+2*DLx*DRMHV+DLx*DLx*DRMHVV) F1=DABS(F1) DLx0=DLx+COEX1*((P-P1)/F1) IF(DLx0.LT.1.) THEN DLx0=3. ENDIF EPS1=DABS((P-P1)/P) IF(EPS1.GE.EPS01+1E-5) COEX1=0.99*COEX1 IF(ABS(EPS1-EPS01).LT.1.E-6) EPS1=1.1E-5 DLx=DLx0 VL=VNL/DLx EPS01=EPS1 IF(J.GT.400) EPS01=0 999 CONTINUE C ------------ END OF THE LIQUIT VOLUM ------------- EPS01=1E6 J=0 C ---------- THE VOLUME OF SATURATION GAS --------- DO 9999 WHILE(EPS02.GT.0.8*EPS) J=J+1 CALL TN_X(XG,TNG,VNG) DGx=VNG/VG TGx=TNG/T CALL D_R_MHV(DGx,TGx,XG,DRMHV) CALL D_R_MHVV(DGx,TGx,XG,DRMHVV) CALL COMP_Z(DGx,TGx,XG,ZG) P1=ZG*DGx*R*T/VNG F=R*T/VNG F1=F*(1.+2*DGx*DRMHV+DGx*DGx*DRMHVV) F1=DABS(F1) DGx0=DGx+COEX2*(P-P1)/F1 IF(DGx0.GT.1) THEN DGx0=0.8 ENDIF IF(DGx0.LT.0) DGx0=1.0E-5 EPS2=DABS((P-P1)/P) IF(ABS(EPS2-EPS02).LT.1.0E-6) EPS2=1.0E-5 IF(EPS2.GE.EPS02) COEX2=0.99*COEX2 IF(COEX2.LT.0.1) COEX2=0.1 EPS02=EPS2 IF(J.GT.400) EPS02=0 DGx=DGx0 VG=VNG/DGx 9999 CONTINUE C ------------ END OF THE SATURATION GAS VOLUM -------- EPS02=1.0E6 C ----------- THE MOLE FRACTION OF GAS AND LIQUID ------------ IF(ABS(X).LT.1.0E-5.OR.ABS(1-X).LT.1.0E-5) THEN TO=500/T DLO=1/(15000*VL) DGO=1/(15000*VG) CALL DIM_S(DLx,TLx,DLO,TO,XL,DIMSL) CALL DIM_S(DGx,TGx,DGO,TO,XG,DIMSG) Y0=(S-R*DIMSL)/(R*(DIMSG-DIMSL)) IF(Y0.GT.1) Y0=1 IF(Y0.LT.0) Y0=0 ELSE CALL FUGA_C(DLx,TLx,XL,FCL1,FCL2) CALL FUGA_C(DGx,TGx,XG,FCG1,FCG2) IF(FCG1.EQ.0) FCG1=1.0E-10 IF(FCL2.EQ.0) FCL2=1.0E-10 XG0=XG+COEX3*((1-(1-XL)*FCL1/FCG1)-XG) XL0=XL+COEX4*(XG*FCG2/FCL2-XL) IF(T.LT.410) THEN X0=1 ELSE X0=0.95 ENDIF IF(XG0.GT.1) XG0=0.99*X0 IF(XG0.LT.0) XG0=0.005*XG IF(XL0.GT.1) XL0=0.99*X0 IF(XL0.LT.0) XL0=0.005*XL IF(XL0.GT.X) XL0=X IF(XG0.LT.X) XG0=X IF(XL0.GT.XG0) XL0=XG0*XG0 IF(XL0.LT.1.0E-8) XL0=1.0E-8 IF(XL0.LT.1.0E-10.AND.XG0.GT.0.1) XG0=0.1 IF(XL0.GT.(1-1.0E-10).AND.XG0.LT.0.9) XG0=1 IF(XG0.LT.1.0E-10.AND.XL0.GT.0.1) XL0=0.1 XG=XG0 XL=XL0 Y0=(X-XL)/(XG-XL) ENDIF EPS3=DABS(Y-Y0) Y=Y0 C --------- END OF MOLE FRACTION OF GAS ----------- C -------- THE CALCULATION OF TEMPERATURE ---------- TO=500/T DLO=1/(15000*VL) DGO=1/(15000*VG) CALL D_I_HTT(TO,XL,DIHTTL) CALL D_R_MHTT(DLx,TLx,XL,DRMHTTL) CALL DIM_S(DLx,TLx,DLO,TO,XL,DIMSL) CALL D_I_HTT(TO,XG,DIHTTG) CALL D_R_MHTT(DGx,TGx,XG,DRMHTTG) CALL DIM_S(DGx,TGx,DGO,TO,XG,DIMSG) S1=R*((1-Y)*DIMSL+Y*DIMSG) FSL=-R*TNL/(T*T)*(TO*500./TNL*DIHTTL+TLx*DRMHTTL) FSG=-R*TNG/(T*T)*(TO*500./TNG*DIHTTG+TGx*DRMHTTG) FS=(1-Y)*FSL+Y*FSG T=T+COEX5*(S-S1)/FS IF(T.GT.TD) T=TD IF(T.LT.TB) T=TB EPS4=ABS(S-S1)/S C ------- END OF THE CALCULATION OF TEMPERATURE ----------- EPS=MAX(EPS1,EPS2,EPS3,EPS4) IF(EPS.GT.0.001) EPS=0.001 IF(I.GT.40) EPS=0 100 CONTINUE RETURN END C ========================================================== C ** TWO PHASE REGION: P,V,T ===> T,XL,XG,VL,VG ** C ========================================================== SUBROUTINE PVX_TVX(P,V,X,T,XL,XG,Y,VL,VG,TB,TD) DOUBLE PRECISION P,V,X,T,XL,XG,Y,VL,VG,X0,TB,TD DOUBLE PRECISION R,TNG,VNG,TNL,VNL,DRMHV,DRMHVV DOUBLE PRECISION DGx,DLx,TGx,TLx DOUBLE PRECISION FCL1,FCL2,FCG1,FCG2 DOUBLE PRECISION ZL,ZG,Y0 DOUBLE PRECISION F,F1,A1,A,B,C,T1 DOUBLE PRECISION DGx0,XG0,XL0 DOUBLE PRECISION DRMHVL,DRMHVG DOUBLE PRECISION DRMHTTG,DRMHTXG,DRMHVTG,DRMHXG,DTNXG,DDNXG DOUBLE PRECISION DRMHTTL,DRMHTXL,DRMHVTL,DRMHXL,DTNXL,DDNXL DOUBLE PRECISION DRMHTL,DRMHTG,HRMHGL,HRMHGG I=1 R=8.314510 EPS=1.0e-3 EPS01=1.0E6 EPS02=1.0E6 C ----------------- TOTAL CYCLE ------------------ DO 100 WHILE(EPS.GT.1.E-4) I=I+1 IF(ABS(X).LT.1E-5) THEN COEX1=1 COEX2=1 COEX3=1 COEX4=1 ELSE IF(T.GT.413.15) THEN IF(T.GT.500) THEN COEX1=10 ELSE COEX1=100*(1-X*X) ENDIF IF(T.GE.500) THEN COEX4=1 COEX2=1 ELSE IF(T.GE.470) THEN COEX4=5/I**0.3 COEX2=1 ELSE COEX4=10/((X*X+0.1)*I**0.3) COEX2=2 ENDIF IF(T.GE.433.15.AND.X.GE.0.9) COEX4=0.8 IF(T.GE.453.15.AND.X.GE.0.8) COEX4=0.8 IF(T.GE.473.15.AND.X.GE.0.7) COEX4=0.8 IF(T.GE.503.15.AND.X.GE.0.6) COEX4=0.5 IF(T.GE.523.15.AND.X.GE.0.5) COEX4=0.5 IF(T.GE.553.15.AND.X.GE.0.4) COEX4=0.5 IF(T.GE.573.15.AND.X.GE.0.3) COEX4=0.5 IF(T.GE.593.15.AND.X.GE.0.2) COEX4=0.5 IF(T.GE.613.15.AND.X.GE.0.8) COEX4=1 ELSE IF(T.GT.343.15) THEN COEX1=100*(1-0.1*X*X) COEX2=1*(1+X) COEX4=10/(X*X*X+0.1)/I**0.3 ELSE IF(T.GT.313.15) THEN COEX1=100 COEX2=1 COEX4=50/(X*X*X+0.1)/I**0.3 ELSE IF(T.GT.265) THEN COEX1=150 COEX2=1 COEX4=50/(X*X*X+0.01)/I**0.3 ELSE IF(T.GT.220) THEN COEX1=200 COEX2=1 COEX4=200/(X**3*I**0.3) ELSE COEX1=300 COEX2=1 COEX4=200/(X**5*(1-X)*I**0.3) ENDIF COEX3=1 ENDIF IF(X.LT.1.AND.X.GT.0) THEN COEX5=(X+0.1)*(1.1-X)/I**0.3 ELSE COEX5=1 ENDIF J=0 C -----------VOLUME CALCULATION OF LIQUID ------------ DO 999 WHILE(EPS01.GT.0.8*EPS) J=J+1 CALL TN_X(XL,TNL,VNL) DLx=VNL/VL TLx=TNL/T CALL D_R_MHV(DLx,TLx,XL,DRMHV) CALL D_R_MHVV(DLx,TLx,XL,DRMHVV) CALL COMP_Z(DLx,TLx,XL,ZL) P1=ZL*DLx*R*T/VNL F=R*T/VNL F1=F*(1.+2*DLx*DRMHV+DLx*DLx*DRMHVV) F1=DABS(F1) DLx0=DLx+COEX1*((P-P1)/F1) IF(DLx0.LT.1.) THEN DLx0=3. ENDIF EPS1=DABS((P-P1)/P) IF(EPS1.GE.EPS01+1E-5) COEX1=0.99*COEX1 IF(ABS(EPS1-EPS01).LT.1.E-6) EPS1=1.1E-5 DLx=DLx0 VL=VNL/DLx EPS01=EPS1 IF(J.GT.400) EPS01=0 999 CONTINUE EPS01=1.0E6 C ------------ END OF THE LIQUIT VOLUM ------------- J=0 C ---------- THE VOLUME OF SATURATION GAS --------- DO 9999 WHILE(EPS02.GT.0.8*EPS) J=J+1 CALL TN_X(XG,TNG,VNG) DGx=VNG/VG TGx=TNG/T CALL D_R_MHV(DGx,TGx,XG,DRMHV) CALL D_R_MHVV(DGx,TGx,XG,DRMHVV) CALL COMP_Z(DGx,TGx,XG,ZG) P1=ZG*DGx*R*T/VNG F=R*T/VNG F1=F*(1.+2*DGx*DRMHV+DGx*DGx*DRMHVV) F1=DABS(F1) DGx0=DGx+COEX2*(P-P1)/F1 IF(DGx0.GT.1) THEN DGx0=0.8 ENDIF IF(DGx0.LT.0) DGx0=1.0E-5 EPS2=DABS((P-P1)/P) IF(ABS(EPS2-EPS02).LT.1.0E-6) EPS2=1.0E-5 IF(EPS2.GE.EPS02) COEX2=0.99*COEX2 IF(COEX2.LT.0.1) COEX2=0.1 EPS02=EPS2 IF(J.GT.400) EPS02=0 DGx=DGx0 VG=VNG/DGx 9999 CONTINUE C -------- END OF THE SATURATION GAS VOLUM -------- EPS02=1.0E6 C --------- THE CALCULATION OF TEMPERATURE -------- IF(ABS(X).LT.1.0E-5.OR.ABS(1-X).LT.1.0E-5) THEN Y0=(V-VL)/(VG-VL) IF(Y0.GT.1) Y0=1 IF(Y0.LT.0) Y0=0 T=TB EPS4=0 ELSE CALL D_R_MHT(DLx,TLx,XL,DRMHTL) CALL D_R_MHT(DGx,TGx,XG,DRMHTG) CALL H_RMHG(DLx,TLx,XL,HRMHGL) CALL H_RMHG(DGx,TGx,XG,HRMHGG) CALL DTV_X(XL,DTNXL,DDNXL) CALL COMP_Z(DLx,TLx,XL,ZL) CALL D_R_MHV(DLx,TLx,XL,DRMHVL) CALL D_R_MHX(DLx,TLx,XL,DRMHXL) CALL D_R_MHTV(DLx,TLx,XL,DRMHVTL) CALL D_R_MHTX(DLx,TLx,XL,DRMHTXL) CALL D_R_MHTT(DLx,TLx,XL,DRMHTTL) CALL DTV_X(XG,DTNXG,DDNXG) CALL COMP_Z(DGx,TGx,XG,ZG) CALL D_R_MHV(DGx,TGx,XG,DRMHVG) CALL D_R_MHX(DGx,TGx,XG,DRMHXG) CALL D_R_MHTV(DGx,TGx,XG,DRMHVTG) CALL D_R_MHTX(DGx,TGx,XG,DRMHTXG) CALL D_R_MHTT(DGx,TGx,XG,DRMHTTG) A1=-DLOG(DABS(ZG/ZL))+ $ (HRMHGG-HRMHGL)+(DGx*DRMHVG-DLx*DRMHVL)+ $ (1-XG)*(DRMHXG+DDNXG/VG*DRMHVG+DTNXG/T*DRMHTG)- $ (1-XL)*(DRMHXL+DDNXL/VL*DRMHVL+DTNXL/T*DRMHTL) A=DLOG(XL/XG) B=TNG/(T*T)*(1/ZG*DGx-DRMHTG-DGx*DRMHVTG-(1-XG)* $ (DRMHTXG+DRMHTG+TGx*DRMHTTG+DGx*DRMHVTG)) C=TNL/(T*T)*(1/ZL*DLx-DRMHTL-DLx*DRMHVTL-(1-XL)* $ (DRMHTXL+DRMHTL+TLx*DRMHTTL+DLx*DRMHVTL)) T1=T+COEX4*(A-A1)/(B-C) EPS4=10*ABS(T1-T)/T1 T=T1 IF(T.GT.TD) T=TD IF(T.LT.TB) T=TB ENDIF C ---------- END OF THE CALCULATION OF TEMPERATURE ----------- C ----------- THE MOLE FRACTION OF GAS AND LIQUID ------------ IF(ABS(X).LT.1.0E-5.OR.ABS(1-X).LT.1.0E-5) THEN XL=X XG=X Y0=(V-VL)/(VG-VL) IF(Y0.GT.1) Y0=1 IF(Y0.LT.0) Y0=0 ELSE CALL FUGA_C(DLx,TLx,XL,FCL1,FCL2) CALL FUGA_C(DGx,TGx,XG,FCG1,FCG2) IF(FCG1.EQ.0) FCG1=1.0E-10 IF(FCL2.EQ.0) FCL2=1.0E-10 XG0=XG+COEX3*((1-(1-XL)*FCL1/FCG1)-XG) Y0=(V-VL)/(VG-VL) IF(ABS(1-Y0).LT.1.0E-5) THEN XL0=XL+COEX3*(XG*FCG2/FCL2-XL) ELSE XL0=(X-Y0*XG)/(1-Y0) ENDIF IF(T.LT.410) THEN X0=1 ELSE X0=0.95 ENDIF IF(XG0.GT.1) XG0=0.99*X0 IF(XG0.LT.0) XG0=0.005*XG IF(XL0.GT.1) XL0=0.99*X0 IF(XL0.LT.0) XL0=0.005*XL IF(XL0.GT.X) XL0=X IF(XG0.LT.X) XG0=X IF(XL0.GT.XG0) XL0=XG0*XG0 IF(XL0.LT.1.0E-8) XL0=1.0E-8 IF(XL0.LT.1.0E-10.AND.XG0.GT.0.1) XG0=0.1 IF(XL0.GT.(1-1.0E-10).AND.XG0.LT.0.9) XG0=1 IF(XG0.LT.1.0E-10.AND.XL0.GT.0.1) XL0=0.1 XG=XG0 XL=XL0 ENDIF EPS3=DABS(Y-Y0) Y=Y0 C --------- END OF MOLE FRACTION OF GAS ----------- EPS=MAX(EPS1,EPS2,EPS3,EPS4) IF(EPS.GT.0.001) EPS=0.001 IF(I.GT.40) EPS=0 100 CONTINUE RETURN END C ---------------------------------------------------- C ** THIS FUNCTION TRANSFERS THE UNIT OF FRACTION ** C ** FROM KMOL TO KG ** C ---------------------------------------------------- REAL FUNCTION AKG(X) REAL X AKG=17.03026*X/(18.015268-X*(18.015268-17.03026)) RETURN END C --------------------------------------------------------- C ** THIS FUNCTION TRANSFERS THE UNIT OF FRACTION ** C ** FROM KG TO KMOL ** C --------------------------------------------------------- REAL FUNCTION AKMOL(X) REAL X AKMOL=18.015268*X/(18.015268*X+17.03026*(1-X)) RETURN END C --------------------------------------------------- C ** The Identification of Ammonia and Water ** C ** I=1: NH3, I=2: WATER ** C --------------------------------------------------- CHARACTER*40 FUNCTION IDENTM(I,A) INTEGER I CHARACTER*1 A 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 ** FUNDAMENTAL PROPERTIES OF PURE COMPONENTS ** c ** I=1: AMMONIA ; I=2: WATER ** c ------------------------------------------------------- REAL FUNCTION FCM(I,A) INTEGER I,KPA,KAS,MESS,KSTAN CHARACTER*1 A CHARACTER*6 PRNAME COMMON/UNIT/ KPA,MESS,KSTAN,KAS PRNAME=' FCM ' IF(I.EQ.1) THEN IF(A.EQ.'M') THEN FCM=17.030260 ELSE IF(A.EQ.'R') THEN FCM=488.189 ELSE IF(A.EQ.'T') THEN FCM=405.40 IF(KPA.EQ.1.OR.KPA.EQ.3) THEN FCM=FCM-273.15 ENDIF ELSE IF(A.EQ.'P') THEN FCM=113.6 ELSE IF(A.EQ.'V') THEN IF(KAS.EQ.1) THEN FCM=0.0044444444 ELSE FCM=0.075690044 ENDIF ENDIF ELSE IF(I.EQ.2) THEN IF(A.EQ.'M') THEN FCM=18.015268 ELSE IF(A.EQ.'R') THEN FCM=461.51805 ELSE IF(A.EQ.'T') THEN FCM=647.096 IF(KPA.EQ.1.OR.KPA.EQ.3) THEN FCM=FCM-273.15 ENDIF ELSE IF(A.EQ.'P') THEN FCM=220.64 ELSE IF(A.EQ.'V') THEN IF(KAS.EQ.1) THEN FCM=0.0031055901 ELSE FCM=0.055948038 ENDIF ENDIF ENDIF RETURN END C ------------------------------------------------ C ** INDICTING THE REGION OF THE MIXTURE ** C ------------------------------------------------ INTEGER FUNCTION KPHASE(TT,PP,ZZ) INTEGER KPA,KAS,MESS,KSTAN DOUBLE PRECISION TB,TD,XLB,XLD,XGB,XGD,VLB,VLD,VGB,VGD REAL TT,PP,ZZ DOUBLE PRECISION T,P,Z COMMON/UNIT/ KPA,MESS,KSTAN,KAS CHARACTER*6 PRNAME PRNAME='KPHASE' IF(KPA.EQ.1) THEN P=1.0E5*PP T=273.15+TT ELSE IF(KPA.EQ.2) THEN P=1.0E5*PP T=TT ELSE IF(KPA.EQ.3) THEN T=273.15+TT P=PP ELSE T=TT P=PP ENDIF IF(KAS.EQ.1) THEN Z=AKMOL(ZZ) ELSE Z=ZZ ENDIF IF(T.LT.TRIP(Z).OR.T.LT.203.) THEN KPHASE=-2 RETURN ENDIF IF(P.GT.1.0E6) THEN KPHASE=-2 RETURN ENDIF IF(P.GT.PRM(Z)) THEN KPHASE=1 RETURN ENDIF CALL BPX_TXV(P,Z,TB,XLB,XGB,VLB,VGB) IF(Z.LT.1.AND.Z.GT.0) THEN CALL DPX_TXV(P,Z,TD,XLD,XGD,VLD,VGD) ELSE TD=TB ENDIF IF(T.LT.TB) THEN KPHASE=1 ELSE IF(T.GE.TB.AND.T.LE.TD) THEN KPHASE=2 ELSE IF(T.GT.TD) THEN KPHASE=3 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 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 KSTAN=STAND KAS=KASC RETURN END C ------------------------------------------------ C ** SATURATION PRESSURE OF PURE COMPONENT ** C ** I=1: AMMONIA, I=2: WATER ** C ------------------------------------------------ REAL FUNCTION PSTM(I,TT) DOUBLE PRECISION T,P,VL,VG REAL TT INTEGER I,KPA,MESS,KSTAN,KAS COMMON /UNIT/ KPA,MESS,KSTAN,KAS IF(KPA.EQ.1.OR.KPA.EQ.3) THEN T=273.15+TT ELSE T=TT ENDIF IF(I.EQ.1) THEN CALL BTX_PXV(T,1.D0,P,1.D0,1.D0,VL,VG) ELSE IF(I.EQ.2) THEN CALL BTX_PXV(T,0.D0,P,0.D0,0.D0,VL,VG) ENDIF IF(KPA.EQ.1.OR.KPA.EQ.2) THEN PSTM=P/1.0E5 ELSE PSTM=P ENDIF RETURN END c ------------------------------------------------------------------- c ** Bubble Point Pressure and other Thermal Properties of Mixture ** c ** [Input] ** c ** TT: Temperature [K],[C] ** c ** X : Composition [kmol/kmol],[kg/kg] ** c ** [OUTPUT] ** c ** J : Error Detection Code ** c ** P : Bubble Point Pressure [Pa],[bar] ** c ** V : Bubble Point Volume ** c ** H : Bubble Point Enthalpy ** c ** S : Bubble Point Entropy ** c ------------------------------------------------------------------- SUBROUTINE SUBPB(J,TT,P,X,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL TT,P,X,V,H,S DOUBLE PRECISION T,Z,PB,XL,XG,VL,VG,HB,SB DOUBLE PRECISION Tx,Dx,To,Do REAL M1,M2,M,R REAL DH1,DH2,DH,DS1,DS2,DS CHARACTER*6 PRNAME COMMON /UNIT/KPA,MESS,KSTAN,KAS M1=17.03026 M2=18.015268 R=8314.510 J=0 PRNAME='SUBPB' IF(KAS.EQ.1) THEN DH1=-143.19E3 DH2=200.0E3 DS1=-0.4718E3 DS2=1.0E3 DH=X*DH1+(1-X)*DH2 DS=X*DS1+(1-X)*DS2 ELSE DH1=-143.19E3*M1 DH2=200.0E3*M2 DS1=-0.4718E3*M1 DS2=1.0E3*M2 DH=X*DH1+(1-X)*DH2 DS=X*DS1+(1-X)*DS2 ENDIF IF((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ELSE T=TT ENDIF IF(KAS.EQ.1) THEN Z=AKMOL(X) ELSE Z=X ENDIF IF(Z.GT.1.OR.Z.LT.0) THEN J=-2 RETURN ENDIF IF(T.GT.TRM(Z).OR.T.LT.TRIP(Z)) THEN J=-2 RETURN ENDIF IF(P.GT.PRM(Z)) THEN J=-2 RETURN ENDIF M=Z*M1+(1.-Z)*M2 CALL BTX_PXV(T,Z,PB,XL,XG,VL,VG) CALL TRAN_TV(T,VL,Z,Tx,Dx,To,Do) CALL DIM_H(Dx,Tx,Do,To,Z,HB) CALL DIM_S(Dx,Tx,Do,To,Z,SB) IF((KPA.EQ.1).OR.(KPA.EQ.2)) THEN P=SNGL(PB/1.0E5) ELSE P=SNGL(PB) ENDIF IF(KAS.EQ.1) THEN HB=R*T*HB/M SB=R*SB/M V=SNGL(1000.*VL/M) ELSE HB=R*T*HB SB=R*SB V=SNGL(1000.*VL) ENDIF IF(KSTAN.EQ.1) THEN H=SNGL(HB)+DH S=SNGL(SB)+DS ELSE H=SNGL(HB) S=SNGL(SB) ENDIF RETURN END c ---------------------------------------------------------------- c ** Dew Point Pressure and other Thermal Properties of Mixture ** c ** [Input] ** c ** TT: Temperature [K],[C] ** c ** Y : Composition [kmol/kmol],[kg/kg] ** c ** [OUTPUT] ** c ** J : Error Detection Code ** c ** P : Dew Point Pressure [Pa],[bar] ** c ** V : Dew Point Volume ** c ** H : Dew Point Enthalpy ** c ** S : Dew Point Entropy ** c ---------------------------------------------------------------- SUBROUTINE SUBPD(J,TT,P,Y,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL TT,P,Y,V,H,S DOUBLE PRECISION T,Z,PD,XL,XV,VL,VV,HD,SD DOUBLE PRECISION Tx,Dx,To,Do REAL M1,M2,M,R REAL DH1,DH2,DS1,DS2,DH,DS CHARACTER*6 PRNAME COMMON /UNIT/KPA,MESS,KSTAN,KAS M1=17.03026 M2=18.015268 R=8314.510 J=0 PRNAME='SUBPD' IF(KAS.EQ.1) THEN DH1=-143.19E3 DH2=200.0E3 DS1=-0.4718E3 DS2=1.00E3 DH=Y*DH1+(1-Y)*DH2 DS=Y*DS1+(1-Y)*DS2 ELSE DH1=-143.19E3*M1 DH2=200.0E3*M2 DS1=-0.4718E3*M1 DS2=1.00E3*M2 DH=Y*DH1+(1-Y)*DH2 DS=Y*DS1+(1-Y)*DS2 ENDIF IF((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ELSE T=TT ENDIF IF(KAS.EQ.1) THEN Z=AKMOL(Y) ELSE Z=Y ENDIF IF(Z.GT.1.OR.Z.LT.0) THEN J=-2 RETURN ENDIF IF(T.GT.TRM(Z).OR.T.LT.TRIP(Z)) THEN J=-2 RETURN ENDIF IF(P.GT.PRM(Z)) THEN J=-2 RETURN ENDIF M=Z*M1+(1.-Z)*M2 IF(DABS(Z-1).LT.1.0E-5) THEN CALL BTX_PXV(T,Z,PD,XL,XV,VL,VV) ELSE CALL DTX_PXV(T,Z,PD,XL,XV,VL,VV) ENDIF CALL TRAN_TV(T,VV,Z,Tx,Dx,To,Do) CALL DIM_H(Dx,Tx,Do,To,Z,HD) CALL DIM_S(Dx,Tx,Do,To,Z,SD) IF((KPA.EQ.1).OR.(KPA.EQ.2)) THEN P=SNGL(PD/1.0E5) ELSE P=SNGL(PD) ENDIF IF(KAS.EQ.1) THEN HD=R*T*HD/M SD=R*SD/M V=SNGL(1000.*VV/M) ELSE HD=R*T*HD SD=R*SD V=SNGL(1000.*VV) ENDIF IF(KSTAN.EQ.1) THEN H=SNGL(HD)+DH S=SNGL(SD)+DS ELSE H=SNGL(HD) S=SNGL(SD) ENDIF RETURN END c ----------------------------------------------------------------- c ** Saturation Properties of Pure Component (Temperature Input) ** c ** Temperature => Pressure, Volume, Enthalpy, Entropy ** c ** [INPUT] ** c ** I : Component I=1: Ammonia; I=2: Water ** c ** TT: Temperature [K],[c] ** c ** [OUTPUT] ** c ** J : Error Detection Code ** c ** PS: Saturation Pressure [Pa],[bar] ** c ** VL: Liquid Volume [m**3/kmol],[m**3/kg] ** c ** VV: Vapor Volume [m**3/kmol],[m**3/kg] ** c ** HL: Liquid Enthalpy [J/kmol],[J/kg] ** c ** HV: Vapor Enthalpy [J/kmol],[J/kg] ** c ** SL: Liquid Entropy [J/(kmol.K)],[J/(kg.K)] ** c ** SV: Vapor Entropy [J/(kmol.K)],[J/(kg.K)] ** c ----------------------------------------------------------------- SUBROUTINE SUBPST(I,J,TT,PS,VL,VV,HL,HV,SL,SV) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL TT,PS,VL,VV,HL,HV,SL,SV,M,R REAL DH,DS DOUBLE PRECISION T,X,P,VLS,VVS,HLS,HVS,SLS,SVS DOUBLE PRECISION TLx,DLx,TLo,DLo,TVx,DVx,TVo,DVo CHARACTER*6 PRNAME COMMON /UNIT/ KPA,MESS,KSTAN,KAS R=8314.510 J=0 PRNAME='SUBPST' IF((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=TT+273.15 ELSE T=TT ENDIF IF(I.EQ.1) THEN M=17.03026 X=1 IF(KAS.EQ.1) THEN DH=-143.19E3 DS=-0.4718E3 ELSE DH=-143.19E3*M DS=-0.4718E3*M ENDIF IF(T.GT.405.65) THEN J=-2 RETURN ENDIF ELSE IF(I.EQ.2) THEN M=18.015268 X=0 IF(KAS.EQ.1) THEN DH=200.0E3 DS=1.0E3 ELSE DH=200.0E3*M DS=1.0E3*M ENDIF IF(T.GT.647.13) THEN J=-2 RETURN ENDIF ENDIF IF(T.LT.TRIP(X)) THEN J=-2 RETURN ENDIF CALL BTX_PXV(T,X,P,X,X,VLS,VVS) CALL TRAN_TV(T,VLS,X,TLx,DLx,TLo,DLo) CALL TRAN_TV(T,VVS,X,TVx,DVx,TVo,DVo) CALL DIM_H(DLx,TLx,DLo,TLo,X,HLS) CALL DIM_H(DVx,TVx,DVo,TVo,X,HVS) CALL DIM_S(DLx,TLx,DLo,TLo,X,SLS) CALL DIM_S(DVx,TVx,DVo,TVo,X,SVS) IF((KPA.EQ.1).OR.(KPA.EQ.2)) THEN PS=SNGL(P/1.0E5) ELSE PS=SNGL(P) ENDIF IF(KAS.EQ.1) THEN HLS=R*T*HLS/M SLS=R*SLS/M VL=SNGL(1000.*VLS/M) HVS=R*T*HVS/M SVS=R*SVS/M VV=SNGL(1000.*VVS/M) ELSE HLS=R*T*HLS SLS=R*SLS VL=SNGL(1000.*VLS) HVS=R*T*HVS SVS=R*SVS VV=SNGL(1000.*VVS) ENDIF IF(KSTAN.EQ.1) THEN HL=SNGL(HLS)+DH SL=SNGL(SLS)+DS HV=SNGL(HVS)+DH SV=SNGL(SVS)+DS ELSE HL=SNGL(HLS) SL=SNGL(SLS) HV=SNGL(HVS) SV=SNGL(SVS) ENDIF RETURN END c ------------------------------------------------------------- c ** Bubble Point Temperature and Other Thermal Properties ** c ** of Mixture (Pressure and X Input) ** c ** [INPUT] ** c ** PP: Pressure [Pa],[bar ** c ** X : Composition [kmol/kmol],[kg/kg] ** c ** [OUTPUT] ** c ** J : Error Detection Code ** c ** T : Bubble Point Temperature ** c ** V : Bubble Point Volume ** c ** H : Bubble Point Enthalpy ** c ** S : Bubble Point Entropy ** c ------------------------------------------------------------- SUBROUTINE SUBTB(J,T,PP,X,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,PP,X,V,H,S,M1,M2,M,R REAL DH1,DH2,DH,DS1,DS2,DS DOUBLE PRECISION Z,P,TB,XL,XG,VL,VG,HB,SB,Tx,Dx,To,Do CHARACTER*6 PRNAME COMMON /UNIT/ KPA,MESS,KSTAN,KAS J=0 PRNAME='SUBTB' M1=17.03026 M2=18.015268 R=8314.510 IF(KAS.EQ.1) THEN Z=AKMOL(X) ELSE Z=X ENDIF IF(KAS.EQ.1) THEN DH1=-143.19E3 DH2=200.0E3 DS1=-0.4718E3 DS2=1.00E3 DH=X*DH1+(1-X)*DH2 DS=X*DS1+(1-X)*DS2 ELSE DH1=-143.19E3*M1 DH2=200.0E3*M2 DS1=-0.4718E3*M1 DS2=1.00E3*M2 DH=X*DH1+(1-X)*DH2 DS=X*DS1+(1-X)*DS2 ENDIF M=Z*M1+(1.-Z)*M2 IF((KPA.EQ.1).OR.(KPA.EQ.2)) THEN P=1.0E5*PP ELSE P=PP ENDIF IF(Z.GT.1.OR.Z.LT.0) THEN J=-2 RETURN ENDIF IF(P.GT.PRM(Z)) THEN J=-2 RETURN ENDIF CALL BPX_TXV(P,Z,TB,XL,XG,VL,VG) IF(TB.LT.TRIP(Z)) TB=TRIP(Z) CALL TRAN_TV(TB,VL,Z,Tx,Dx,To,Do) CALL DIM_H(Dx,Tx,Do,To,Z,HB) CALL DIM_S(Dx,Tx,Do,To,Z,SB) IF(KAS.EQ.1) THEN HB=R*TB*HB/M SB=R*SB/M V=SNGL(1000.*VL/M) ELSE HB=R*TB*HB SB=R*SB V=SNGL(1000.*VL) ENDIF IF((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=SNGL(TB-273.15) ELSE T=SNGL(TB) ENDIF IF(KSTAN.EQ.1) THEN H=SNGL(HB)+DH S=SNGL(SB)+DS ELSE H=SNGL(HB) S=SNGL(SB) ENDIF RETURN END c ---------------------------------------------------------- c ** Dew Point Temperature and Other Thermal Properties ** c ** of Mixture (Pressure and X Input) ** c ** [INPUT] ** c ** PP: Pressure [Pa],[bar] ** c ** Y : Composition [kmol/kmol],[kg/kg] ** c ** [OUTPUT] ** c ** J : Error Detection Code ** c ** T : Dew Point Temperature [K],[C] ** c ** V : Dew Point Volume ** c ** H : Dew Point Enthalpy ** c ** S : Dew Point Entropy ** c ---------------------------------------------------------- SUBROUTINE SUBTD(J,T,PP,Y,V,H,S) INTEGER J,KPA,MESS,KSTAN,KAS REAL T,PP,Y,V,H,S,M1,M2,M,R REAL DH1,DH2,DH,DS1,DS2,DS DOUBLE PRECISION Z,P,TD,XL,XG,VL,VG,HD,SD,Tx,Dx,To,Do CHARACTER*6 PRNAME COMMON /UNIT/ KPA,MESS,KSTAN,KAS J=0 PRNAME='SUBTD' M1=17.03026 M2=18.015268 R=8314.510 IF(KAS.EQ.1) THEN DH1=-143.19E3 DH2=200.0E3 DS1=-0.4718E3 DS2=1.00E3 DH=Y*DH1+(1-Y)*DH2 DS=Y*DS1+(1-Y)*DS2 ELSE DH1=-143.19E3*M1 DH2=200.0E3*M2 DS1=-0.4718E3*M1 DS2=1.00E3*M2 DH=Y*DH1+(1-Y)*DH2 DS=Y*DS1+(1-Y)*DS2 ENDIF IF(KAS.EQ.1) THEN Z=AKMOL(Y) ELSE Z=Y ENDIF M=Z*M1+(1.-Z)*M2 IF((KPA.EQ.1).OR.(KPA.EQ.2)) THEN P=1.0E5*PP ELSE P=PP ENDIF IF(Z.GT.1.OR.Z.LT.0) THEN J=-2 RETURN ENDIF IF(P.GT.PRM(Z)) THEN J=-2 RETURN ENDIF IF(DABS(Z).LT.1.0E-5.OR.DABS(1.-Z).LT.1.E-5) THEN CALL BPX_TXV(P,Z,TD,XL,XG,VL,VG) ELSE CALL DPX_TXV(P,Z,TD,XL,XG,VL,VG) ENDIF IF(TD.LT.TRIP(Z)) TD=TRIP(Z) CALL TRAN_TV(TD,VG,Z,Tx,Dx,To,Do) CALL DIM_H(Dx,Tx,Do,To,Z,HD) CALL DIM_S(Dx,Tx,Do,To,Z,SD) IF(KAS.EQ.1) THEN HD=R*TD*HD/M SD=R*SD/M V=SNGL(1000.*VG/M) ELSE HD=R*TD*HD SD=R*SD V=SNGL(1000.*VG) ENDIF IF((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T=SNGL(TD-273.15) ELSE T=SNGL(TD) ENDIF IF(KSTAN.EQ.1) THEN H=SNGL(HD)+DH S=SNGL(SD)+DS ELSE H=SNGL(HD) S=SNGL(SD) ENDIF RETURN END c -------------------------------------------------------------- c ** Saturation Properties of Pure Component (Pressure Input) ** c ** [INPUT] ** c ** I : Component I=1: Ammonia, I=2: Water ** c ** PP: Pressure [Pa],[bar] ** c ** [OUTPUT] ** c ** J : Error Detection Code ** c ** TS: Saturation Temperature [K],[C] ** c ** VL: Saturation Liquid Volume ** c ** VV: Saturation Vapor Volume ** c ** HL: Saturation Liquid Enthalpy ** c ** HV: Saturation Vapor Enthalpy ** c ** SL: Saturation Liauid Entropy ** c ** SV: Saturation Vapor Entropy ** c -------------------------------------------------------------- SUBROUTINE SUBTSP(I,J,TS,PP,VL,VV,HL,HV,SL,SV) INTEGER I,J,KPA,MESS,KSTAN,KAS REAL TS,PP,VL,VV,HL,HV,SL,SV,M,R REAL DH,DS DOUBLE PRECISION T,Z,P,VLS,VVS,HLS,HVS,SLS,SVS DOUBLE PRECISION TLx,DLx,TLo,DLo,TVx,DVx,TVo,DVo CHARACTER*6 PRNAME COMMON /UNIT/ KPA,MESS,KSTAN,KAS R=8314.510 J=0 PRNAME='SUBTSP' IF((KPA.EQ.1).OR.(KPA.EQ.2)) THEN P=1.0E5*PP ELSE P=PP ENDIF IF(I.EQ.1) THEN M=17.03026 Z=1 IF(KAS.EQ.1) THEN DH=-143.19E3 DS=-0.4718E3 ELSE DH=-143.19E3*M DS=-0.4718E3*M ENDIF IF(P.GT.113.60E5) THEN J=-2 RETURN ENDIF ELSE IF(I.EQ.2) THEN M=18.015268 Z=0 IF(KAS.EQ.1) THEN DH=200.0E3 DS=1.00E3 ELSE DH=200.0E3*M DS=1.00E3*M ENDIF IF(P.GT.220.64E5) THEN J=-2 RETURN ENDIF ENDIF IF(Z.GT.1.OR.Z.LT.0) THEN J=-2 RETURN ENDIF CALL BPX_TXV(P,Z,T,Z,Z,VLS,VVS) CALL TRAN_TV(T,VLS,Z,TLx,DLx,TLo,DLo) CALL TRAN_TV(T,VVS,Z,TVx,DVx,TVo,DVo) CALL DIM_H(DLx,TLx,DLo,TLo,Z,HLS) CALL DIM_H(DVx,TVx,DVo,TVo,Z,HVS) CALL DIM_S(DLx,TLx,DLo,TLo,Z,SLS) CALL DIM_S(DVx,TVx,DVo,TVo,Z,SVS) IF(KAS.EQ.1) THEN HLS=R*T*HLS/M SLS=R*SLS/M VL=SNGL(1000.*VLS/M) HVS=R*T*HVS/M SVS=R*SVS/M VV=SNGL(1000.*VVS/M) ELSE HLS=R*T*HLS SLS=R*SLS VL=SNGL(1000.*VLS) HVS=R*T*HVS SVS=R*SVS VV=SNGL(1000.*VVS) ENDIF IF((KPA.EQ.1).OR.(KPA.EQ.3)) THEN TS=SNGL(T-273.15) ELSE TS=SNGL(T) ENDIF IF(KSTAN.EQ.1) THEN HL=SNGL(HLS)+DH SL=SNGL(SLS)+DS HV=SNGL(HVS)+DH SV=SNGL(SVS)+DS ELSE HL=SNGL(HLS) SL=SNGL(SLS) HV=SNGL(HVS) SV=SNGL(SVS) ENDIF RETURN END c ---------------------------------------------------------- c ** SUBROUTINE SUBXY ** c ** Thermal Properties of Mixture at VLE Region ** c ** [INPUT] ** c ** TT: Temperature [K],[C] ** c ** PP: Pressure [Pa],[bar] ** c ** [OUTPUT] ** c ** J : Error Detection Code ** c ** X : Composition of Liquid ** c ** Y : Composition of Vapor ** c ** VL: Liquid Volume [m**3/kmol],[m**3/kg] ** c ** VV: Vapor Volume [m**3/kmol],[m**3/kg] ** c ** HL: Liquid Enthalpy [J/kmol],[J/kg] ** c ** HV: Vapor Enthalpy [J/kmol],[J/kg] ** c ** SL: Liquid Entropy [J/(kmol.K)][J/(kg.K)** c ** SV: Vapor Entropy [J/(kmol.K)][J/(kg.K) ** c ---------------------------------------------------------- SUBROUTINE SUBXY(J,TT,PP,X,Y,VL,VV,HL,HV,SL,SV) INTEGER J,KPA,MESS,KSTAN,KAS REAL M1,M2,ML,MV,R REAL DH1,DH2,DS1,DS2,DHL,DHV,DSL,DSV REAL TT,PP,X,Y,VL,VV,HL,HV,SL,SV DOUBLE PRECISION T,P,XL,XV,VLE,VVE,HLE,HVE,SLE,SVE DOUBLE PRECISION DLx,TLx,DLo,TLo,DVx,TVx,DVo,TVo CHARACTER PRNAME COMMON /UNIT/ KPA,MESS,KSTAN,KAS J=0 PRNAME='SUBXY' M1=17.03026 M2=18.015268 R=8314.510 IF(KAS.EQ.1) THEN DH1=-143.19E3 DH2=200.0E3 DS1=-0.4718E3 DS2=1.00E3 ELSE DH1=-143.19E3*M1 DH2=200.0E3*M2 DS1=-0.4718E3*M1 DS2=1.00E3*M2 ENDIF IF(KPA.EQ.1) THEN P=1.0E5*PP T=273.15+TT ELSE IF(KPA.EQ.2) THEN P=1.0E5*PP T=TT ELSE IF(KPA.EQ.3) THEN P=PP T=273.15+TT ELSE P=PP T=TT ENDIF IF(T.GT.647.096) THEN J=-2 RETURN ENDIF IF(T.LT.203.15) THEN J=-2 RETURN ENDIF IF(P.GT.22.064E6) THEN J=-2 RETURN ENDIF CALL SPT_XV(P,T,XL,XV,VLE,VVE) CALL TRAN_TV(T,VLE,XL,TLx,DLx,TLo,DLo) CALL TRAN_TV(T,VVE,XV,TVx,DVx,TVo,DVo) CALL DIM_H(DLx,TLx,DLo,TLo,XL,HLE) CALL DIM_H(DVx,TVx,DVo,TVo,XV,HVE) CALL DIM_S(DLx,TLx,DLo,TLo,XL,SLE) CALL DIM_S(DVx,TVx,DVo,TVo,XV,SVE) ML=XL*M1+(1-XL)*M2 MV=XV*M1+(1-XV)*M2 IF(KAS.EQ.1) THEN HLE=R*T*HLE/ML SLE=R*SLE/ML VL=SNGL(1000.*VLE/ML) HVE=R*T*HVE/MV SVE=R*SVE/MV VV=SNGL(1000.*VVE/MV) X=AKG(SNGL(XL)) Y=AKG(SNGL(XV)) ELSE HLE=R*T*HLE SLE=R*SLE VL=SNGL(1000.*VLE) HVE=R*T*HVE SVE=R*SVE VV=SNGL(1000.*VVE) X=SNGL(XL) Y=SNGL(XV) ENDIF IF(KSTAN.EQ.1) THEN DHL=X*DH1+(1-X)*DH2 DHV=Y*DH1+(1-Y)*DH2 DSL=X*DS1+(1-X)*DS2 DSV=Y*DS1+(1-Y)*DS2 HL=SNGL(HLE)+DHL SL=SNGL(SLE)+DSL HV=SNGL(HVE)+DHV SV=SNGL(SVE)+DSV ELSE HL=SNGL(HLE) SL=SNGL(SLE) HV=SNGL(HVE) SV=SNGL(SVE) ENDIF RETURN END c -------------------------------------------------------- c ** Saturation Temperature Calculation of ** c ** Pure Component (Pressure INPUT) ** c ** I=1: Ammonia ** c ** I=2: Water ** c -------------------------------------------------------- REAL FUNCTION TSPM(I,PP) INTEGER I,KPA,MESS,KSTAN,KAS REAL PP CHARACTER*6 PRNAME DOUBLE PRECISION Z,P,T,XL,XG,VL,VG COMMON /UNIT/ KPA,MESS,KSTAN,KAS PRNAME=' TSPM ' IF(I.EQ.1) Z=1 IF(I.EQ.2) Z=0 IF(KPA.EQ.1.OR.KPA.EQ.2) THEN P=1.0E5*PP ELSE P=PP ENDIF CALL BPX_TXV(P,Z,T,XL,XG,VL,VG) IF(KPA.EQ.1.OR.KPA.EQ.3) THEN TSPM=SNGL(T-273.15) ELSE TSPM=SNGL(T) ENDIF RETURN END c ---------------------------------------------------------- c ** SUBROUTINE SUBPUR ** c ** Single Phase Thermal Properties of Pure Component ** c ** [INPUT] ** c ** I : Component, I=1,Ammonia;I=2,Water ** c ** TT: Temperature [K],[C] ** c ** PP: Pressure [Pa],[bar] ** c ** [OUTPUT] ** c ** J : Error Detection Code ** c ** V : Volume [m**3/kmol],[m**3/kg] ** c ** H : Enthalpy [J/kmol],[J/kg] ** c ** S : Entropy [J/(kmol.K)],[J/(kg.K)] ** c ---------------------------------------------------------- SUBROUTINE SUBPUR(I,J,TT,PP,V,H,S) INTEGER I,J,KPA,MESS,KSTAN,KAS,K REAL TT,PP,V,H,S,M,R REAL DH,DS DOUBLE PRECISION P,T,VS,HS,SS,Z,Dx,Tx,Do,To CHARACTER*6 PRNAME COMMON /UNIT/KPA,MESS,KSTAN,KAS J=0 R=8314.510 PRNAME='SUBPUR' IF(I.EQ.1) THEN M=17.03026 Z=1 IF(KAS.EQ.1) THEN DH=-143.19E3 DS=-0.4718E3 ELSE DH=-143.19E3*M DS=-0.4718E3*M ENDIF ELSE IF(I.EQ.2) THEN M=18.015268 Z=0 IF(KAS.EQ.1) THEN DH=200.0E3 DS=1.00E3 ELSE DH=200.E3*M DS=1.0E3*M ENDIF ENDIF IF(KPA.EQ.1) THEN P=1.0E5*PP T=273.15+TT ELSE IF(KPA.EQ.2) THEN P=1.0E5*PP T=TT ELSE IF(KPA.EQ.3) THEN P=PP T=273.15+TT ELSE P=PP T=TT ENDIF K=KPHASE(PP,TT,SNGL(Z)) IF(K.EQ.-1.OR.K.EQ.-2) THEN J=K RETURN ELSE J=K IF((K.EQ.1).OR.(K.EQ.2)) THEN CALL PTX_V(P,T,Z,VS,0) ELSE IF(K.EQ.3) THEN CALL PTX_V(P,T,Z,VS,1) ENDIF CALL TRAN_TV(T,VS,Z,Tx,Dx,To,Do) CALL DIM_H(Dx,Tx,Do,To,Z,HS) CALL DIM_S(Dx,Tx,Do,To,Z,SS) IF(KAS.EQ.1) THEN HS=R*T*HS/M SS=R*SS/M V=SNGL(1000.*VS/M) ELSE HS=R*T*HS SS=R*SS V=SNGL(1000.*VS) ENDIF IF(KSTAN.EQ.1) THEN H=SNGL(HS)+DH S=SNGL(SS)+DS ELSE H=SNGL(HS) S=SNGL(SS) ENDIF RETURN ENDIF END c ------------------------------------------------------------ c ** SUBROUTINE SUBMIX ** c ** Calculation of the Properties of Ammonia-Water Mixture ** c ** [Parameters] ** c ** I : Function Label: ** 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,H,S ** c ** I=4 : P,Z,V => T,H,S ** c ** J : Error Detection Code ** c ** TT: Temperature [K],[C] ** c ** PP: Pressure [Pa],[Bar] ** c ** ZZ: Composition [kg/kg],[kmol/kmol] ** c ** V : Volume [m**3/kg],[m**3/kmol] ** c ** H : Enthalpy [J/kmol],[J/kg] ** c ** S : Entropy [J/(kmol.K)],[J/(kmol.K)] ** c ------------------------------------------------------------ SUBROUTINE SUBMIX(I,J,TT,PP,ZZ,V,H,S) INTEGER I,J,KPA,MESS,KSTAN,KAS,PHASE REAL TT,PP,ZZ,V,H,S,M1,M2,M,R REAL DH1,DH2,DH,DS1,DS2,DS DOUBLE PRECISION T,P,Z,VS,HS,SS,XL,XV,VL,VV,Y,HL,HV,SL,SV DOUBLE PRECISION TB,TD,XLD,XLB,VLD,VLB,XVD,XVB,VVD,VVB DOUBLE PRECISION HB,HD,SB,SD,TBx,DBx,TBo,DBo,TDx,DDx,TDo,DDo DOUBLE PRECISION Tx,Dx,To,Do,TLx,DLx,TLo,DLo,TVx,DVx,TVo,DVo CHARACTER*6 PRNAME COMMON /UNIT/ KPA,MESS,KSTAN,KAS J=0 PRNAME='SUBMIX' M1=17.03026 M2=18.015268 R=8314.510 IF(KAS.EQ.1) THEN DH1=-143.19E3 DH2=200.0E3 DS1=-0.4718E3 DS2=1.00E3 DH=ZZ*DH1+(1-ZZ)*DH2 DS=ZZ*DS1+(1-ZZ)*DS2 ELSE DH1=-143.19E3*M1 DH2=200.0E3*M2 DS1=-0.4718E3*M1 DS2=1.00E3*M2 DH=ZZ*DH1+(1-ZZ)*DH2 DS=ZZ*DS1+(1-ZZ)*DS2 ENDIF IF(KPA.EQ.1) THEN P=1.0E5*PP ELSE IF(KPA.EQ.2) THEN P=1.0E5*PP ELSE IF(KPA.EQ.3) THEN P=PP ELSE P=PP ENDIF IF(KAS.EQ.1) THEN Z=AKMOL(ZZ) ELSE Z=ZZ ENDIF IF(Z.GT.1.OR.Z.LT.0) THEN J=-2 RETURN ENDIF IF(P.GT.40E6) THEN J=-2 RETURN ENDIF M=Z*M1+(1-Z)*M2 IF(I.EQ.1) THEN IF(KPA.EQ.1) THEN T=273.15+TT ELSE IF(KPA.EQ.2) THEN T=TT ELSE IF(KPA.EQ.3) THEN T=273.15+TT ELSE T=TT ENDIF IF(T.LT.TRIP(Z).OR.T.LT.203.) THEN J=-2 RETURN ENDIF IF(P.GT.PRM(Z)) THEN PHASE=1 ELSE CALL BPX_TXV(P,Z,TB,XLB,XVB,VLB,VVB) IF((Z.GT.0).AND.(Z.LT.1)) THEN CALL DPX_TXV(P,Z,TD,XLD,XVD,VLD,VVD) ELSE TD=TB XLD=XLB XVD=XVB VLD=VLB VVD=VVB ENDIF IF(T.LT.TB) THEN PHASE=1 ELSE IF(T.GT.TD) THEN PHASE=3 ELSE IF(T.EQ.TB) THEN PHASE=12 ELSE IF(T.EQ.TD) THEN PHASE=23 ELSE IF(T.GT.TB.AND.T.LT.TD) THEN PHASE=2 ENDIF ENDIF IF(PHASE.EQ.1) THEN CALL PTX_V(P,T,Z,VS,0) Y=2 ELSE IF(PHASE.EQ.3) THEN CALL PTX_V(P,T,Z,VS,1) Y=-1 ELSE IF(PHASE.EQ.12) THEN XL=XLB XV=XVB VL=VLB VV=VVB Y=0 ELSE IF(PHASE.EQ.23) THEN XL=XLD XV=XVD VL=VLD VV=VVD Y=1 ELSE IF(PHASE.EQ.2) THEN CALL SPT_XV(P,T,XL,XV,VL,VV) IF(Z.LT.1.AND.Z.GT.0) THEN Y=(Z-XL)/(XV-XL) ELSE Y=0 ENDIF ENDIF IF(PHASE.EQ.1.OR.PHASE.EQ.3) THEN CALL TRAN_TV(T,VS,Z,Tx,Dx,To,Do) CALL DIM_H(Dx,Tx,Do,To,Z,HS) CALL DIM_S(Dx,Tx,Do,To,Z,SS) HS=R*T*HS SS=R*SS ELSE IF(PHASE.EQ.12.OR.PHASE.EQ.23.OR.PHASE.EQ.2) THEN CALL TRAN_TV(T,VL,XL,TLx,DLx,TLo,DLo) CALL TRAN_TV(T,VV,XV,TVx,DVx,TVo,DVo) CALL DIM_H(DLx,TLx,DLo,TLo,XL,HL) CALL DIM_H(DVx,TVx,DVo,TVo,XV,HV) CALL DIM_S(DLx,TLx,DLo,TLo,XL,SL) CALL DIM_S(DVx,TVx,DVo,TVo,XV,SV) HL=R*T*HL HV=R*T*HV SL=R*SL SV=R*SV HS=Y*HV+(1-Y)*HL SS=Y*SV+(1-Y)*SL VS=Y*VV+(1-Y)*VL ENDIF ELSE IF(I.EQ.2) THEN IF(KSTAN.EQ.1) THEN H=H-DH ENDIF IF(KAS.EQ.1) THEN HS=H*M ELSE HS=H ENDIF IF(P.GT.PRM(Z)) THEN PHASE=1 ELSE CALL BPX_TXV(P,Z,TB,XLB,XVB,VLB,VVB) IF((Z.GT.0).AND.(Z.LT.1)) THEN CALL DPX_TXV(P,Z,TD,XLD,XVD,VLD,VVD) ELSE TD=TB XLD=XLB XVD=XVB VLD=VLB VVD=VVB ENDIF CALL TRAN_TV(TB,VLB,XLB,TBx,DBx,TBo,DBo) CALL TRAN_TV(TD,VVD,XVD,TDx,DDx,TDo,DBo) CALL DIM_H(DBx,TBx,DBo,TBo,XLB,HB) CALL DIM_H(DDx,TDx,DDo,TDo,XVD,HD) HB=R*TB*HB HD=R*TD*HD IF(HS.LT.HB) THEN PHASE=1 ELSE IF(HS.GT.HD) THEN PHASE=3 ELSE IF(HS.EQ.HB) THEN PHASE=12 ELSE IF(HS.EQ.HD) THEN PHASE=23 ELSE IF(HS.GT.HB.AND.HS.LT.HD) THEN PHASE=2 ENDIF ENDIF IF(PHASE.EQ.1) THEN HS=HS/1000 CALL PHX_VT(P,HS,Z,VS,T,0) HS=1000*HS Y=2 ELSE IF(PHASE.EQ.3) THEN HS=HS/1000 CALL PHX_VT(P,HS,Z,VS,T,1) HS=1000*HS Y=-1 ELSE IF(PHASE.EQ.12) THEN XL=XLB XV=XVB VL=VLB VV=VVB T=TB Y=0 ELSE IF(PHASE.EQ.23) THEN XL=XLD XV=XVD VL=VLD VD=VVD T=TD Y=1 ELSE IF(PHASE.EQ.2) THEN IF(Z.LT.1.AND.Z.GT.0) THEN T=TB+(HS-HB)*(TD-TB)/(HD-HB) XL=XLB+(HS-HB)*(XLD-XLB)/(HD-HB) XV=XVB+(HS-HB)*(XVD-XVB)/(HD-HB) VL=VLB+(HS-HB)*(VLD-VLB)/(HD-HB) VV=VVB+(HS-HB)*(VVD-VVB)/(HD-HB) Y=(Z-XL)/(XV-XL) HS=HS/1000 CALL PHX_TVX(P,HS,Z,T,XL,XV,Y, $ VL,VV,TB,TD) HS=1000*HS ELSE T=TB VL=VLB VV=VVB XL=Z XV=Z ENDIF ENDIF IF(PHASE.EQ.1.OR.PHASE.EQ.3) THEN CALL TRAN_TV(T,VS,Z,Tx,Dx,To,Do) CALL DIM_S(Dx,Tx,Do,To,Z,SS) SS=R*SS ELSE IF(PHASE.EQ.12.OR.PHASE.EQ.23.OR.PHASE.EQ.2) THEN CALL TRAN_TV(T,VL,XL,TLx,DLx,TLo,DLo) CALL TRAN_TV(T,VV,XV,TVx,DVx,TVo,DVo) CALL DIM_S(DLx,TLx,DLo,TLo,XL,SL) CALL DIM_S(DVx,TVx,DVo,TVo,XV,SV) CALL DIM_H(DLx,TLx,DLo,TLo,XL,HL) CALL DIM_H(DVx,TVx,DVo,TVo,XV,HV) SL=R*SL SV=R*SV HL=R*T*HL HV=R*T*HV IF((Z.EQ.1.OR.Z.EQ.0).AND.(PHASE.EQ.2)) THEN Y=(HS-HL)/(HV-HL) ENDIF VS=Y*VV+(1-Y)*VL SS=Y*SV+(1-Y)*SL ENDIF ELSE IF(I.EQ.3) THEN IF(KSTAN.EQ.1) THEN S=S-DS ENDIF IF(KAS.EQ.1) THEN SS=S*M ELSE SS=S ENDIF IF(P.GT.PRM(Z)) THEN PHASE=1 ELSE CALL BPX_TXV(P,Z,TB,XLB,XVB,VLB,VVB) IF((Z.GT.0).AND.(Z.LT.1)) THEN CALL DPX_TXV(P,Z,TD,XLD,XVD,VLD,VVD) ELSE TD=TB XLD=XLB XVD=XVB VLD=VLB VVD=VVB ENDIF CALL TRAN_TV(TB,VLB,XLB,TBx,DBx,TBo,DBo) CALL TRAN_TV(TD,VVD,XVD,TDx,DDx,TDo,DDo) CALL DIM_S(DBx,TBx,DBo,TBo,XLB,SB) CALL DIM_S(DDx,TDx,DDo,TDo,XVD,SD) SB=R*SB SD=R*SD IF(SS.LT.SB) THEN PHASE=1 ELSE IF(SS.GT.SD) THEN PHASE=3 ELSE IF(SS.EQ.SB) THEN PHASE=12 ELSE IF(SS.EQ.SD) THEN PHASE=23 ELSE IF(SS.GT.SB.AND.SS.LT.SD) THEN PHASE=2 ENDIF ENDIF IF(PHASE.EQ.1) THEN SS=SS/1000 CALL PSX_VT(P,SS,Z,VS,T) SS=1000*SS Y=2 ELSE IF(PHASE.EQ.3) THEN SS=SS/1000 CALL PSX_VT(P,SS,Z,VS,T) SS=1000*SS Y=-1 ELSE IF(PHASE.EQ.12) THEN XL=XLB XV=XVB VL=VLB VV=VVB T=TB Y=0 ELSE IF(PHASE.EQ.23) THEN XL=XLD XV=XVD VL=VLD VD=VVD T=TD Y=1 ELSE IF(PHASE.EQ.2) THEN IF(Z.LT.1.AND.Z.GT.0) THEN T=TB+(SS-SB)*(TD-TB)/(SD-SB) XL=XLB+(SS-SB)*(XLD-XLB)/(SD-SB) XV=XVB+(SS-SB)*(XVD-XVB)/(SD-SB) VL=VLB+(SS-SB)*(VLD-VLB)/(SD-SB) VV=VVB+(SS-SB)*(VVD-VVB)/(SD-SB) Y=(Z-XL)/(XV-XL) SS=SS/1000 CALL PSX_TVX(P,SS,Z,T,XL,XV,Y, $ VL,VV,TB,TD) SS=1000*SS ELSE T=TB VL=VLB VV=VVB XL=Z XV=Z ENDIF ENDIF IF(PHASE.EQ.1.OR.PHASE.EQ.3) THEN CALL TRAN_TV(T,VS,Z,Tx,Dx,To,Do) CALL DIM_H(Dx,Tx,Do,To,Z,HS) HS=R*T*HS ELSE IF(PHASE.EQ.12.OR.PHASE.EQ.23.OR.PHASE.EQ.2) THEN CALL TRAN_TV(T,VL,XL,TLx,DLx,TLo,DLo) CALL TRAN_TV(T,VV,XV,TVx,DVx,TVo,DVo) CALL DIM_H(DLx,TLx,DLo,TLo,XL,HL) CALL DIM_H(DVx,TVx,DVo,TVo,XV,HV) CALL DIM_S(DLx,TLx,DLo,TLo,XL,SL) CALL DIM_S(DVx,TVx,DVo,TVo,XV,SV) HL=R*T*HL HV=R*T*HV SL=R*SL SV=R*SV IF((Z.EQ.1.OR.Z.EQ.0).AND.(PHASE.EQ.2)) THEN Y=(SS-SL)/(SV-SL) ENDIF VS=Y*VV+(1-Y)*VL HS=Y*HV+(1-Y)*HL ENDIF ELSE IF(I.EQ.4) THEN IF(KAS.EQ.1) THEN VS=V*M/1000. ELSE VS=V/1000 ENDIF IF(P.GT.PRM(Z)) THEN PHASE=1 ELSE CALL BPX_TXV(P,Z,TB,XLB,XVB,VLB,VVB) IF((Z.GT.0).AND.(Z.LT.1)) THEN CALL DPX_TXV(P,Z,TD,XLD,XVD,VLD,VVD) ELSE TD=TB XLD=XLB XVD=XVB VLD=VLB VVD=VVB ENDIF CALL TRAN_TV(TB,VLB,XLB,TBx,DBx,TBo,DBo) CALL TRAN_TV(TD,VVD,XVD,TDx,DDx,TDo,DDo) IF(VS.LT.VLB) THEN PHASE=1 ELSE IF(VS.GT.VVD) THEN PHASE=3 ELSE IF(VS.EQ.VLB) THEN PHASE=12 ELSE IF(VS.EQ.VVD) THEN PHASE=23 ELSE IF(VS.GT.VLB.AND.VS.LT.VVD) THEN PHASE=2 ENDIF ENDIF IF(PHASE.EQ.1) THEN CALL PVX_T(P,VS,Z,T) Y=2 ELSE IF(PHASE.EQ.3) THEN CALL PVX_T(P,VS,Z,T) Y=-1 ELSE IF(PHASE.EQ.12) THEN XL=XLB XV=XVB VL=VLB VV=VVB T=TB Y=0 ELSE IF(PHASE.EQ.23) THEN XL=XLD XV=XVD VL=VLD VD=VVD T=TD Y=1 ELSE IF(PHASE.EQ.2) THEN IF(Z.LT.1.AND.Z.GT.0) THEN T=TB+(VS-VLB)*(TD-TB)/(VVD-VLB) XL=XLB+(VS-VLB)*(XLD-XLB)/(VVD-VLB) XV=XVB+(VS-VLB)*(XVD-XVB)/(VVD-VLB) VL=VLB+(VS-VLB)*(VLD-VLB)/(VVD-VLB) VV=VVB+(VS-VLB)*(VVD-VVB)/(VVD-VLB) Y=(Z-XL)/(XV-XL) CALL PVX_TVX(P,VS,Z,T,XL,XV,Y, $ VL,VV,TB,TD) Y=(Z-XL)/(XV-XL) ELSE T=TB VL=VLB VV=VVB XL=Z XV=Z ENDIF ENDIF IF(PHASE.EQ.1.OR.PHASE.EQ.3) THEN CALL TRAN_TV(T,VS,Z,Tx,Dx,To,Do) CALL DIM_H(Dx,Tx,Do,To,Z,HS) CALL DIM_S(Dx,Tx,Do,To,Z,SS) HS=R*T*HS SS=R*SS ELSE IF(PHASE.EQ.12.OR.PHASE.EQ.23.OR.PHASE.EQ.2) THEN CALL TRAN_TV(T,VL,XL,TLx,DLx,TLo,DLo) CALL TRAN_TV(T,VV,XV,TVx,DVx,TVo,DVo) CALL DIM_H(DLx,TLx,DLo,TLo,XL,HL) CALL DIM_H(DVx,TVx,DVo,TVo,XV,HV) CALL DIM_S(DLx,TLx,DLo,TLo,XL,SL) CALL DIM_S(DVx,TVx,DVo,TVo,XV,SV) HL=R*T*HL HV=R*T*HV SL=R*SL SV=R*SV IF((Z.EQ.1.OR.Z.EQ.0).AND.(PHASE.EQ.2)) THEN Y=(VS-VL)/(VV-VL) ENDIF VS=VS HS=Y*HV+(1-Y)*HL SS=Y*SV+(1-Y)*SL ENDIF ENDIF IF(KAS.EQ.1) THEN H=SNGL(HS/M) S=SNGL(SS/M) V=SNGL(1000.*VS/M) ELSE H=SNGL(HS) S=SNGL(SS) V=SNGL(1000.*VS) ENDIF IF(KSTAN.EQ.1) THEN H=H+DH S=S+DS ELSE H=H S=S ENDIF IF(KPA.EQ.1.OR.KPA.EQ.3) THEN TT=T-273.15 ELSE TT=T ENDIF RETURN END 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 ** Critical Temperature of Mixtures (Ammonia + Water) ** C ** INPUT X ** C -------------------------------------------------------- REAL FUNCTION TRM(X) DOUBLE PRECISION X REAL TC01,TC02,KT,ALFA,TC12 TC01=647.096 TC02=405.40 KT=0.9648407 ALFA=1.125455 TC12=0.5*KT*(TC01+TC02) TRM=(1-X)*(1-X)*TC01+X*X*TC02+2*X*(1-X**ALFA)*TC12 RETURN END C ----------------------------------------------------- C ** Critical Pressure of Mixtures (Ammonia + Water) ** C ** INPUT X ** C ----------------------------------------------------- REAL FUNCTION PRM(X) DOUBLE PRECISION X REAL PC01,PC02,PC12,KP,ALFA PC01=22.064E6 PC02=11.36E6 KP=2.65 ALFA=1.23 PC12=0.5*(PC01+PC02) PRM=(1-X)*(1-X)*PC01+X*X*PC02+KP*X*(1-X**ALFA)*PC12 RETURN END c ================================================== c ** The Triple-Point Line ** c ================================================== REAL FUNCTION TRIP(X) DOUBLE PRECISION X REAL C11,C12,C13,C21,C31,C32,C41,C42 C11=-0.3439823 C12=-1.3274271 C13=-274.973 C21=-4.987368 C31=-4.886151 C32=10.37298 C41=-0.323998 C42=-15.87560 IF((X.LE.0.33367).AND.(X.GE.0.0)) THEN TRIP=273.1599*(1.0+C11*X+C12*X**2+C13*X**7) ELSEIF((X.GT.0.33367).AND.(X.LE.0.58396)) THEN TRIP=193.549*(1+C21*(x-0.5)**2) ELSEIF((X.GT.0.58396).AND.(X.LE.0.81473)) THEN TRIP=194.38*(1+C31*(X-2.0/3.0)**2+C32*(x-2.0/3.0)**3) ELSEIF((X.GT.0.81473).AND.(X.LE.1.0)) THEN TRIP=195.495*(1+C41*(1-X)+C42*(1-x)**4) ELSE TRIP=-1000.0 ENDIF RETURN END C ------------------------------------------------ C ** INDICTING THE REGION OF THE MIXTURE ** C ------------------------------------------------ INTEGER FUNCTION IPHASE(TT,PP,ZZ) INTEGER KPA,KAS,MESS,KSTAN DOUBLE PRECISION TB,TD,XLB,XLD,XGB,XGD,VLB,VLD,VGB,VGD REAL TT,PP,ZZ DOUBLE PRECISION T,P,Z COMMON/UNIT/ KPA,MESS,KSTAN,KAS CHARACTER*6 PRNAME PRNAME='IPHASE' IF(KPA.EQ.1) THEN P=1.0E5*PP T=273.15+TT ELSE IF(KPA.EQ.2) THEN P=1.0E5*PP T=TT ELSE IF(KPA.EQ.3) THEN T=273.15+TT P=PP ELSE T=TT P=PP ENDIF IF(KAS.EQ.1) THEN Z=AKMOL(ZZ) ELSE Z=ZZ ENDIF IF(T.LT.TRIP(Z).OR.T.LT.203.) THEN IPHASE=-2 RETURN ENDIF IF(P.GT.1.0E6) THEN IPHASE=-2 RETURN ENDIF IF(P.GT.PRM(Z)) THEN IPHASE=1 RETURN ENDIF CALL BPX_TXV(P,Z,TB,XLB,XGB,VLB,VGB) IF(Z.LT.1.AND.Z.GT.0) THEN CALL DPX_TXV(P,Z,TD,XLD,XGD,VLD,VGD) ELSE TD=TB ENDIF IF(T.LT.TB) THEN IPHASE=1 ELSE IF(T.GE.TB.AND.T.LE.TD) THEN IPHASE=2 ELSE IF(T.GT.TD) THEN IPHASE=1 ENDIF RETURN END