[CHNOSZ-commits] r951 - in pkg/CHNOSZ: . src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 6 10:25:44 CET 2026
Author: jedick
Date: 2026-01-06 10:25:43 +0100 (Tue, 06 Jan 2026)
New Revision: 951
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/src/H2O92D.f
Log:
Remove unused Fortran subroutines
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2026-01-06 08:41:43 UTC (rev 950)
+++ pkg/CHNOSZ/DESCRIPTION 2026-01-06 09:25:43 UTC (rev 951)
@@ -1,6 +1,6 @@
Date: 2026-01-06
Package: CHNOSZ
-Version: 2.2.0-17
+Version: 2.2.0-18
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/src/H2O92D.f
===================================================================
--- pkg/CHNOSZ/src/H2O92D.f 2026-01-06 08:41:43 UTC (rev 950)
+++ pkg/CHNOSZ/src/H2O92D.f 2026-01-06 09:25:43 UTC (rev 951)
@@ -16,6 +16,12 @@
*** Abandoned: 8 November 1991
***
***********************************************************************
+***
+*** NOTE: Unused code paths in CHNOSZ were removed by jmd on 2026-01-06
+*** - See also previous changes in comments by jmd
+*** The original version of this file is in the inst/src directory
+*** of the CHNOSZ repo on GitHub.
+***
*
* specs - Input unit, triple point, saturation, and option specs:
*
@@ -87,15 +93,7 @@
END IF
- IF (useLVS) THEN
- Dens(1) = states(3)
- CALL LVSeqn(specs(6),specs(7),specs(5),
- 1 states(1),states(2),Dens,specs(9))
- Dens(1) = Dens(1) / 1.0d3
- IF (specs(6) .EQ. 1) THEN
- Dens(2) = Dens(2) / 1.0d3
- END IF
- ELSE
+ IF (.NOT. useLVS) THEN
Dens(1) = states(3) / 1.0d3
CALL HGKeqn(specs(6),specs(7),specs(5),
1 states(1),states(2),Dens,specs(9))
@@ -134,7 +132,7 @@
IMPLICIT DOUBLE PRECISION (a-h,o-z)
INTEGER useLVS, epseqn
- LOGICAL valspc, valTD, valTP
+ LOGICAL valspc, valTP
COMMON /tolers/ TTOL, PTOL, DTOL, XTOL, EXPTOL, FPTOL
COMMON /units/ ft, fd, fvd, fvk, fs, fp, fh, fst, fc
@@ -164,9 +162,7 @@
Pcrit = Pc * 1.0d1
IF (isat .EQ. 0) THEN
- IF (iopt .EQ. 1) THEN
- valid = valTD(T,D,isat,epseqn)
- ELSE
+ IF (iopt .NE. 1) THEN
valid = valTP(T,P)
END IF
ELSE
@@ -211,91 +207,6 @@
*****************************************************************
-*** valTD - Returns "TRUE" if T-D defines liquid or vapor H2O
-* within validity limits of the HGK equation of state;
-* returns "FALSE" otherwise.
-
- LOGICAL FUNCTION valTD(T,D,isat,epseqn)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- INTEGER epseqn
-
- COMMON /tolers/ TTOL, PTOL, DTOL, XTOL, EXPTOL, FPTOL
- COMMON /aconst/ wm, gascon, tz, aa, zb, dzb, yb, uref, sref
- COMMON /RTcurr/ rt
- COMMON /crits/ Tc, rhoC, Pc, Pcon, Ucon, Scon, dPcon
- COMMON /tpoint/ Utr, Str, Htr, Atr, Gtr,
- 1 Ttr, Ptripl, Dltrip, Dvtrip
- COMMON /HGKbnd/ Ttop, Tbtm, Ptop, Pbtm, Dtop, Dbtm
- COMMON /liqice/ sDli1, sPli1, sDli37, sPli37, sDIB30,
- 1 Tli13, Pli13, Dli13, TnIB30, DnIB30
- COMMON /coefs/ a(20), q(20), x(11)
- COMMON /satur/ Dliq, Dvap, DH2O, iphase
-
- SAVE
-
- EQUIVALENCE (TmnLVS, x(1))
-
-
- IF ((T-FPTOL .GT. Ttop) .OR. (T+FPTOL .LT. Tbtm) .OR.
- 1 (D-FPTOL .GT. Dtop) .OR. (D+FPTOL .LT. Dbtm)) THEN
- valTD = .FALSE.
- RETURN
- END IF
-
- Tcrit = Tc - 273.15d0
- Ttripl = Ttr - 273.15d0
-
- IF ((T+FPTOL .GE. Tcrit) .OR.
- 1 ((T .GE. TnIB30) .AND. (D .GE. Dltrip))) THEN
- Dlimit = sDIB30 * (T-TnIB30) + Dtop
- valTD = (D-FPTOL .LE. Dlimit)
- ELSE
- IF (D-FPTOL .LE. Dltrip) THEN
- IF (T .GE. Ttripl) THEN
- valTD = .TRUE.
- Tk = T + 273.15d0
- IF (Tk .LT. TmnLVS) THEN
- rt = gascon * Tk
- CALL pcorr(0,Tk,Ps,Dl,Dv,epseqn)
- ELSE
- istemp = 1
- DH2O = 0.0d0
- P = Pfind(istemp,Tk,DH2O)
- CALL denLVS(istemp,Tk,P)
- Dv = Dvap / 1.0d3
- Dl = Dliq / 1.0d3
- END IF
- IF ((D .GE. Dv) .AND. (D. LE. Dl)) THEN
- isat = 1
- END IF
- ELSE
- P = Psublm(T)
- PMPa = P / 1.0d1
- Tk = T + 273.15d0
- Dguess = PMPa / Tk / 0.4d0
- rt = gascon * Tk
- CALL bb(Tk)
- CALL denHGK(Dsublm,PMPa,Dguess,Tk,dPdD)
- valTD = (D-FPTOL .LE. Dsublm)
- END IF
- ELSE
- IF (D .LE. Dli13) THEN
- Dlimit = sDli1 * (T-Tli13) + Dli13
- valTD = (D+FPTOL .GE. Dlimit)
- ELSE
- Dlimit = sDli37 * (T-Tli13) + Dli13
- valTD = (D-FPTOL .LE. Dlimit)
- END IF
- END IF
- END IF
-
- RETURN
- END
-
-*****************************************************************
-
*** valTP - Returns "TRUE" if T-P defines liquid or vapor H2O
* within validity limits of the HGK equation of state;
* returns "FALSE" otherwise.
@@ -332,9 +243,6 @@
valTP = (P+FPTOL .GE. Plimit)
*** a hack?, keep valTP TRUE at temperatures below the 0.01 deg C triple point jmd 20130202
valTP = .TRUE.
- ELSE
- Psubl = Psublm(T)
- valTP = (P-FPTOL .LE. Psubl)
END IF
END IF
@@ -343,30 +251,6 @@
*****************************************************************
-*** Psublm - Returns Psublimation(T) computed from the
-* equation given by Washburn (1924): Monthly
-* Weather Rev., v.52, pp.488-490.
-
- DOUBLE PRECISION FUNCTION Psublm(Temp)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- SAVE
-
-
- T = Temp + 2.731d2
-
- PmmHg = power(1.0d1, (-2.4455646d3/T + 8.2312d0*DLOG10(T) -
- 1 1.677006d-2*T + 1.20514d-5*T*T - 6.757169d0))
-
-*** convert mmHg to bars ***
- Psublm = PmmHg * 1.33322d-3
-
- RETURN
- END
-
-************************************************************************
-
*** HGKcon - Constant parameters for the H2O equation of state
* given by Haar, Gallagher, & Kell (1984):
* bp, bq = b(j), B(j) from Table A.1, p.272
@@ -639,14 +523,7 @@
IF (llim .AND. ulim) THEN
crtreg = .TRUE.
ELSE
- IF (llim .AND. (T .LE. Tmin2)) THEN
- isat1 = 1
- ddummy = 0.0d0
- Pstest = Pfind(isat1,T,ddummy)
- crtreg = (P .LE. Pstest)
- ELSE
- crtreg = .FALSE.
- END IF
+ crtreg = .FALSE.
END IF
END IF
END IF
@@ -733,8 +610,6 @@
IF (isat .EQ. 1) THEN
IF (iopt .EQ. 1) THEN
CALL pcorr(itripl,Temp,Pres,Dens(1),Dens(2),epseqn)
- ELSE
- CALL tcorr(itripl,Temp,Pres,Dens(1),Dens(2),epseqn)
END IF
ELSE
IF ((Temp .GT. Tc) .OR. (Temp .LT. Ttripl) .OR.
@@ -1436,82 +1311,6 @@
***********************************************************************
-*** TsHGK - Returns Tsaturation(P).
-
- DOUBLE PRECISION FUNCTION TsHGK(p)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- SAVE
-
-
- TsHGK = 0.0d0
-
- IF (p .GT. 22.05d0) RETURN
-
- k = 0
- pl = 2.302585d0 + DLOG(p)
- tg = 372.83d0 +
- 1 pl*(27.7589d0 + pl*(2.3819d0 + pl*(0.24834d0 +
- 2 pl*0.0193855d0)))
-
- 1 IF (tg .LT. 273.15d0) tg = 273.15d0
- IF (tg .GT. 647.00d0) tg = 647.00d0
-
- IF (k .GE. 8) THEN
- TsHGK = tg
- ELSE
- k = k + 1
- pp = PsHGK(tg)
- dp = TdPsdT(tg)
- IF (ABS(1.0d0 - pp/p) .LT. 1.0d-5) THEN
- TsHGK = tg
- ELSE
- tg = tg * (1.0d0 + (p - pp)/dp)
- GO TO 1
- END IF
- END IF
-
- RETURN
- END
-
-***********************************************************************
-
-*** TdPsdT - Returns T*(dPsat/dT).
-
- DOUBLE PRECISION FUNCTION TdPsdT(t)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- DOUBLE PRECISION a(8)
-
- SAVE
-
- DATA a /-.78889166d1, .25514255d1, -.6716169d1, .33239495d2,
- 1 -.10538479d3, .17435319d3, -.14839348d3, .48631602d2/
-
-
- v = t / 647.25d0
- w = 1.0 - v
- b = 0.0d0
- c = 0.0d0
-
- DO i=1,8
- z = i
- y = a(i) * power(w,(z + 1.0d0)/2.0d0)
- c = c + y/w*(0.5d0 - 0.5d0*z - 1.0d0/v)
- b = b + y
- END DO
-
- q = b / v
- TdPsdT = 22.093d0 * DEXP(q) * c
-
- RETURN
-
- END
-
-***********************************************************************
-
*** corr - Computes liquid and vapor densities (dliq & dvap)
* and (Gl-Gv)/RT (delg) for T-P conditions on or
* near the saturation surface.
@@ -1591,1148 +1390,6 @@
RETURN
END
-************************************************************
-
-*** tcorr - Computes Tsaturation(P) (t) and liquid and vapor
-* densities (dl & dv) from refinement of an initial
-* approximation (TsHGK(p)) in accord with Gl = Gv.
-
- SUBROUTINE tcorr(itripl,t,p,dl,dv,epseqn)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- INTEGER epseqn
-
- COMMON /aconst/ wm, gascon, tz, aa, zb, dzb, yb, uref, sref
- COMMON /RTcurr/ rt
-
- SAVE
-
-
- t = TsHGK(p)
- IF (t .EQ. 0.0d0) RETURN
- dl = 0.0d0
- dv = 0.0d0
-
- 1 rt = t * gascon
- CALL corr(itripl,t,p,dl,dv,delg,epseqn)
-
- dp = delg * gascon * t / (1.0d0/dv - 1.0d0/dl)
- t = t * (1.0d0 - dp/TdPsdT(t))
-
- IF (DABS(delg) .GT. 1.0d-4) GO TO 1
-
- RETURN
- END
-
-***************************************************************
-
-*** LVSeqn - Computes thermodynamic and transport properties of
-* critical region H2O (369.85-419.85 degC,
-* 0.20-0.42 gm/cm3) from the fundamental equation given
-* by Levelt Sengers, et al (1983): J.Phys.Chem.Ref.Data,
-* V.12, No.1, pp.1-28.
-
- SUBROUTINE LVSeqn(isat,iopt,itripl,T,P,Dens,epseqn)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- PARAMETER (NPROP = 23)
-
- DOUBLE PRECISION wprops(NPROP), wpliq(NPROP), Dens(2)
- LOGICAL cpoint
- INTEGER epseqn
-
- COMMON /coefs/ a(20), q(20), x(11)
- COMMON /crits/ Tc, rhoc, Pc, Pcon, Ucon, Scon, dPcon
- COMMON /therm/ AE, GE, U, H, Entrop, Cp, Cv, betaw, alphw,
- 1 heat, Speed
- COMMON /satur/ Dliq, Dvap, DH2O, iphase
- COMMON /param/ r1, th1
- COMMON /units/ ft, fd, fvd, fvk, fs, fp, fh, fst, fc
- COMMON /wpvals/ wprops, wpliq
-
- SAVE
-
-
- cpoint = .FALSE.
- DH2O = Dens(1)
-
- 10 CALL LVSsat(iopt,isat,T,P,DH2O)
-
- IF ((isat .NE. 0) .OR. (iopt .NE. 1)) CALL denLVS(isat,T,P)
-
- IF (isat .EQ. 0) THEN
- Dens(1) = DH2O
- ELSE
- Dens(1) = Dliq
- Dens(2) = Dvap
- END IF
-
- IF (isat .EQ. 0) THEN
- CALL thmLVS(isat,T,r1,th1)
- CALL dimLVS(isat,itripl,th1,T,P*1.0d1,dl,dv,wprops,epseqn)
- IF (cpoint) THEN
- CALL cpswap
- Dens(1) = cdens
- Dens(2) = cdens
- isat = 1
- iopt = ioptsv
- END IF
- ELSE
- th1 = -1.0d0
- CALL thmLVS(isat,T,r1,th1)
- CALL dimLVS(isat,itripl,th1,T,P*1.0d1,dl,dv,wprops,epseqn)
- th1 = 1.0d0
- CALL thmLVS(isat,T,r1,th1)
- CALL dimLVS(isat,itripl,th1,T,P*1.0d1,dl,dv,wpliq,epseqn)
- IF (dl .EQ. dv) THEN
- cpoint = .TRUE.
- cdens = dl
- T = 647.0670000003d0
- P = 22.0460000008d0
- ioptsv = iopt
- iopt = 2
- isat = 0
- GO TO 10
- END IF
- END IF
-
- END
-
-*********************************************************************
-
-*** cpswap - Load critical point A, G, U, H, S, Vs, Di, ZB,
-* albe values from wpliq into wprops and
-* approximations to critical Cv, Cp, alpha, beta,
-* visc, tcond, Prndtl, tdiff, visck, YB, QB, XB,
-* daldT, st values from wprops into wpliq.
-
- SUBROUTINE cpswap
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- PARAMETER (NPROP = 23)
-
- INTEGER aw, gw, sw, uw, hw, cvw, cpw, vsw, alw, bew,
- 1 diw, viw, tcw, stw, tdw, Prw, vikw, albew,
- 2 ZBw, YBw, QBw, dalwdT, XBw
- DOUBLE PRECISION wprops(NPROP), wpliq(NPROP)
-
- COMMON /wpvals/ wprops, wpliq
- COMMON /units/ ft, fd, fvd, fvk, fs, fp, fh, fst, fc
-
- SAVE
-
- DATA aw, gw, sw, uw, hw, cvw, cpw, vsw, alw, bew, diw, viw,
- 1 tcw, stw, tdw, Prw, vikw, albew, ZBw, YBw, QBw,
- 2 dalwdT, XBw
- 2 / 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
- 3 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23 /
-
-
- wprops(aw) = wpliq(aw)
- wprops(gw) = wpliq(gw)
- wprops(sw) = wpliq(sw)
- wprops(uw) = wpliq(uw)
- wprops(hw) = wpliq(hw)
- wprops(diw) = wpliq(diw)
- wprops(ZBw) = wpliq(ZBw)
- wprops(stw) = wpliq(stw)
-
- wpliq(cvw) = wprops(cvw)
- wpliq(cpw) = wprops(cpw)
- wpliq(alw) = wprops(alw)
- wpliq(bew) = wprops(bew)
- wpliq(YBw) = wprops(YBw)
- wpliq(QBw) = wprops(QBw)
- wpliq(XBw) = wprops(XBw)
- wpliq(tcw) = wprops(tcw)
- wpliq(tdw) = wprops(tdw)
- wpliq(Prw) = wprops(Prw)
- wpliq(dalwdT) = wprops(dalwdT)
- wpliq(albew) = wprops(albew)
-
- wprops(vsw) = 0.429352766443498d2 * fs
- wprops(viw) = 1.0d6
- wprops(vikw) = 1.0d6
-
- wpliq(vsw) = wprops(vsw)
- wpliq(viw) = wprops(viw)
- wpliq(vikw) = wprops(vikw)
-
- END
-
-*********************************************************************
-
-*** LVSsat - If isat=1, computes Psat(T) or Tsat(P) (iopt=1,2).
-* If isat=0, checks whether T-D or T-P (iopt=1,2)
-* falls on or within TOL of the liq-vap surface; if so,
-* isat <- 1 and T <- Tsat.
-
- SUBROUTINE LVSsat(iopt,isat,T,P,D)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- COMMON /tolers/ TTOL, PTOL, DTOL, XTOL, EXPTOL, FPTOL
- COMMON /crits/ Tc, rhoc, Pc, Pcon, Ucon, Scon, dPcon
-
- SAVE
-
- DATA ERRTOL, TCTOL / 1.0d-12, 1.0d-2 /
-
-
- IF (isat .EQ. 1) THEN
- IF (iopt .EQ. 1) THEN
- P = Pfind(isat,T,D)
- END IF
- T = TsLVS(isat,P)
- ELSE
- IF (iopt .EQ. 1) THEN
- P = Pfind(isat,T,D)
- END IF
- IF (P-ERRTOL .GT. Pc) THEN
- RETURN
- ELSE
- CALL backup
- Tsat = TsLVS(isat,P)
- IF (DABS(Tsat-T) .LT. TCTOL) THEN
- T = Tsat
- isat = 1
- ELSE
- CALL restor
- END IF
- END IF
- END IF
-
- RETURN
- END
-
-*********************************************************************
-
-*** denLVS - Calculates DH2O(T,P) or Dvap,Dliq(T,P) from the
-* Levelt Sengers, et al (1983) critical region
-* equation of state.
-
- SUBROUTINE denLVS(isat,T,P)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- DOUBLE PRECISION s(2), sd(2)
-
- COMMON /coefs/ a(20), q(20), x(11)
- COMMON /crits/ Tc, rhoc, Pc, Pcon, Ucon, Scon, dPcon
- COMMON /satur/ Dliq, Dvap, DH2O, iphase
- COMMON /therm/ AE, GE, U, H, Entrop, Cp, Cv, betaw, alphw,
- 1 heat, Speed
- COMMON /param/ r1, th1
- COMMON /deri2/ dPdD, dPdT
-
- SAVE
-
- EQUIVALENCE (Dmin, x(4)), (Dmax, x(5)), (pw11, q(9)),
- 1 (xk0, a(7)), (xk1, a(12))
-
- IF (isat .EQ. 0) THEN
- DH2O = rhoc
- DO 10 i=1,20
- Pnext = Pfind(isat,T,DH2O)
- Pdif = Pnext - P
- IF (iphase .EQ. 2) THEN
- IF (DABS(Pdif) .LE. 0.0d0) THEN
- RETURN
- ELSE
- END IF
- IF (Pdif .LT. 0.0d0) THEN
- DH2O = Dmax
- ELSE
- DH2O = Dmin
- END IF
- ELSE
- delD = -Pdif/dPdD
- DH2O = DH2O + delD
- IF (DH2O .LT. Dmin) DH2O = Dmin
- IF (DH2O .GT. Dmax) DH2O = Dmax
- IF (DABS(delD/DH2O) .LT. 1.0d-6) RETURN
- END IF
- 10 CONTINUE
- ELSE
- Tw = -Tc/T
- dTw = 1.0d0 + Tw
-
- CALL ss(r1,th1,s,sd)
- rho1 = 1.0d0+pw11*dTw+a(1)*(s(1)+s(2))
- rho2 = xk0*power(r1,a(6)) + xk1*power(r1,q(16))
-
- Dvap = rhoc * (rho1 - rho2)
- Dliq = rhoc * (rho1 + rho2)
-
- RETURN
- END IF
-
- RETURN
- END
-
-*********************************************************************
-
-*** TsLVS - Returns saturation T(P)
-
- DOUBLE PRECISION FUNCTION TsLVS(isat,P)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- COMMON /therm/ AE, GE, U, H, Entrop, Cp, Cv, betaw, alphw,
- 1 heat, Speed
- COMMON /satur/ Dliq, Dvap, DH2O, iphase
- COMMON /crits/ Tc, rhoc, Pc, Pcon, Ucon, Scon, dPcon
- COMMON /deri2/ dPdD, dPdT
-
- SAVE
-
-
- TsLVS2 = Tc - 1.0d0
- D = rhoc
-
- DO 10 i=1,20
- Pnext = Pfind(isat,TsLVS2,D)
- dT = (Pnext - P)/dPdT
- TsLVS2 = TsLVS2 - dT
- IF (TsLVS2 .GT. Tc) THEN
- TsLVS2 = Tc
- ELSE
- IF (DABS(dT/TsLVS2) .LT. 1.0d-8) THEN
- GO TO 20
- ELSE
- END IF
- END IF
- 10 CONTINUE
-
- 20 TsLVS = TsLVS2
-
- RETURN
- END
-
-*********************************************************************
-
-*** Pfind - Returns P(T,D). Computes (dP/dD)T when invoked by SUB
-* Dfind (isat=0) and (dP/dT)D when invoked by SUB TsLVS
-* (isat=1). Also computes 1st & 2nd partial derivatives
-* the singular part of the potential (Delta P tilde) that
-* are used in SUB thmLVS.
-
- DOUBLE PRECISION FUNCTION Pfind(isat,T,D)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- DOUBLE PRECISION s(2), xk(2), sd(2)
-
- COMMON /coefs/ a(20), q(20), x(11)
- COMMON /crits/ Tc, rhoc, Pc, Pcon, Ucon, Scon, dPcon
- COMMON /satur/ Dliq, Dvap, DH2O, iphase
- COMMON /therm/ AE, GE, U, H, Entrop, Cp, Cv, betaw, alphw,
- 1 heat, Speed
- COMMON /param/ r1, th1
- COMMON /tolers/ TTOL, PTOL, DTOL, XTOL, EXPTOL, FPTOL
- COMMON /deriv/ amu, s, Pw, Tw, dTw, dM0dT, dP0dT,
- 1 d2PdM2, d2PdMT, d2PdT2, p0th, p1th, xk
- COMMON /deri2/ dPdD, dPdT
-***************************************
- COMMON /abc2/ r, th
-***************************************
-
- SAVE
-
- EQUIVALENCE (Pw1, a(5)), (Pw2, a(4)), (Pw3, a(2)),
- 1 (amc, a(13)), (am1, a(14)), (am2, a(15)),
- 2 (am3, a(16)), (p00, q(11)), (p20, q(12)),
- 3 (p40, q(13)), (p01, q(18)), (p21, q(19)),
- 4 (p41, q(20)), (aa, a(10)), (xk0, a(7)),
- 5 (xk1, a(12)), (pw11,q(9)), (alpha,q(10)),
- 6 (alhi,q(15)), (besq,a(9))
-
- xk(1) = xk0
- xk(2) = xk1
- IF (DABS(T-Tc) .LT. FPTOL) T = Tc
- Tee = (T-Tc)/Tc
- Tw = -Tc/T
- dTw = 1.0d0 + Tw
-
- IF (isat .EQ. 0) THEN
- rho = D / rhoc
- CALL conver(rho,Tee,amu,th1,r1,rho1,s,rhodi,err)
- ELSE
- th1 = -1.0d0
- th = th1
- r1 = dTw/(1.0d0-besq)
- r = r1
- CALL ss(r1,th1,s,sd)
- rho = th1 * (xk0*power(r1,a(6)) +
- 1 xk1*power(r1,q(16))) +
- 2 a(1)*(s(1)+s(2))
- rho = 1.0d0+pw11*dTw+rho
- amu = 0.0d0
- D = rho * rhoc
- END IF
-
- tt1 = th1*th1
- tt2 = tt1*tt1
-
- Pw0 = 1.0d0+dTw*(Pw1+dTw*(Pw2+dTw*Pw3))
-
- IF (isat .EQ. 0) THEN
- Pwmu = amu*rhodi
- ELSE
- Pwmu = 0.0d0
- END IF
-
- p0th = p00+p20*tt1+p40*tt2
- p1th = p01+p21*tt1+p41*tt2
-
- dPw0 = xk0*p0th*power(r1,2.0d0-alpha)
- dPw1 = xk1*p1th*power(r1,2.0d0-alhi)
- dPw = aa*(dPw0+dPw1)
-
- Pw = Pw0 + Pwmu + dPw
-
- Pfind = Pw * Pcon * T
-
- IF (DABS(th1) .LT. 1.0d0) THEN
- iphase = 1
- ELSE
- iphase = 2
-
- dP0dT = Pw1+dTw*(2.0d0*Pw2+3.0d0*Pw3*dTw)
- dM0dT = am1+dTw*(2.0d0*am2+3.0d0*am3*dTw)
- Uw = dP0dT-rho*dM0dT+pw11*amu+s(1)+s(2)
-
- dPdTcd = Uw + rho*dM0dT
- dPwdTw = Pw - Tw*dPdTcd
-
- dPdT = Pcon * dPwdTw
-
- END IF
-
- CALL aux(r1,th1,d2PdT2,d2PdMT,d2PdM2,aa,xk,sd,Cvcoex)
-
- IF (iphase .EQ. 1) dPdD = dPcon * D * T / d2PdM2
-
- RETURN
- END
-
-***************************************************************
-
-*** aux - Calculates some second derivatives of the
-* anomalous part of the equation of state.
-
- SUBROUTINE aux(r1,th1,d2PdT2,d2PdMT,d2PdM2,aa,xk,sd,Cvcoex)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- DOUBLE PRECISION xk(2), s(2), sd(2), w(2), y(2), z(2), coex(2)
-
- COMMON /coefs/ a(20), q(20), x(11)
-
- SAVE
-
- EQUIVALENCE (cc, a(1)), (beta, a(6)), (besq,a(9)),
- 1 (delta,a(11)), (alpha,q(10)), (s00, a(17)),
- 2 (s20, a(18)), (s01, a(19)), (s21, a(20))
-
-
- deli = 0.0d0
- s(1) = s00+s20*th1*th1
- s(2) = s01+s21*th1*th1
- sd(1) = 2.0*th1*s20
- sd(2) = 2.0*th1*s21
- ww = 0.0d0
- yy = 0.0d0
- zz = 0.0d0
- gamma = beta*(delta-1.0d0)
- tt1 = th1*th1
- ter = 2.0d0*beta*delta-1.0d0
- g = (1.0+(besq*ter-3.0)*tt1 - besq*(ter-2.0)*tt1*tt1)
- Cvcoex = 0.0d0
-
- DO 30 i=1,2
- alhi = alpha - deli
- beti = beta + deli
- gami = gamma - deli
- IF (r1 .NE. 0.0d0) THEN
- w(i) = (1.0-alhi)*(1.0-3.0*tt1)*s(i) -
- 1 beta*delta*(1.0-tt1)*th1*sd(i)
- w(i) = (w(i)*power(r1,-alhi))/g
- w(i) = w(i) * xk(i)
- ww = ww + w(i)
-
- y(i) = beti*(1.0d0-3.0d0*tt1)*th1 -
- 1 beta*delta*(1.0d0-tt1)*th1
- y(i) = (y(i)*power(r1,beti-1.0d0)) * xk(i) / g
- yy = yy + y(i)
-
- z(i) = 1.0d0-besq*(1.0d0-(2.0d0*beti))*tt1
- z(i) = (z(i)*power(r1,-gami)) * xk(i) / g
- zz = zz + z(i)
-
- a1 = (beta*(delta-3.0d0)-3.0d0*deli-besq*alhi*gami) /
- 1 (2.0d0*besq*besq*(2.0d0-alhi)*(1.0d0-alhi)*alhi)
- a2 = 1+((beta*(delta-3.0d0)-3.0d0*deli-besq*alhi*ter) /
- 1 (2.0d0*besq*(1.0d0-alhi)*alhi))
- a2 = -a2
-
- a4 = 1.0d0+((ter-2.0d0)/(2.0d0*alhi))
- f1 = a1 + a2 + a4
-
- coex(i) = ((2.0d0-alhi)*(1.0d0-alhi)*power(r1,-alhi) *
- 1 f1*xk(i))
- Cvcoex = Cvcoex + coex(i)
- END IF
- deli = 0.5d0
- 30 CONTINUE
-
- d2PdT2 = aa * ww
- d2PdMT = yy + aa*cc*ww
- d2PdM2 = zz/aa + 2.0d0*cc*yy + cc*cc*aa*ww
-
- RETURN
- END
-
-***************************************************************
-
-*** conver - Transforms T,D to parametric variables r,theta
-* according to the revised and scaled equations.
-
- SUBROUTINE conver(rho,Tee,amu,th1,r1,rho1s,s1,rhodi,error1)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- DOUBLE PRECISION s1(2), sd(2)
-
- COMMON /coefs/ a(20), q(20), x(11)
- COMMON /crits/ Tc, rhoc, Pc, Pcon, Ucon, Scon, dPcon
-**************************************************************
- COMMON /abc2/ r, th
-**************************************************************
-
- SAVE
-
- EQUIVALENCE (beta,a(6)), (delta,a(11)), (xk1, a(12)),
- 1 (cc, a(1)), (alhi, q(15)), (alpha,q(10)),
- 2 (besq,a(9)), (p11, q(9)), (deli, q(14)),
- 3 (p1w, q(18)), (p2w, q(19)), (p4w, q(20)),
- 4 (aa, a(10)), (xk0, a(7)), (s00, a(17)),
- 5 (s20, a(18)), (betai,q(16))
-
-
- Tstar = Tee + 1.0d0
- dtstin = 1.0d0 - (1.0d0 / Tstar)
- r1 = dtstin
-
- IF (dtstin .LE. 0.0d0) THEN
- r1 = dtstin/(1.0d0-besq)
- th1 = 1.0d0
- ELSE
- th1 = 0.0d0
- END IF
-
- CALL ss(r1,th1,s1,sd)
-
- rhodi = 1.0d0 + p11*dtstin
- rhodit = rhodi + cc*s1(1) + cc*s1(2)
- drho = rho - rhodit
- amu = 0.0d0
-
- IF (dtstin .LE. 0.0d0) THEN
- rho1co = xk0*power(r1,beta) + xk1*power(r1,betai)
- twofaz = rho1co
- IF (DABS(drho) .LE. twofaz) THEN
- rho1s = DSIGN(rho1co,drho) + cc*s1(1)
- th1 = DSIGN(1.00d0,drho)
- error1 = 1.0d0
- r = r1
- th = th1
- RETURN
- END IF
- END IF
-
- IF (drho .EQ. 0.0d0) THEN
- th1 = 0.0d0
- r1 = dtstin
- rho1s = cc*s1(1)
- END IF
-
-*** rule for first pass ***
-
- y1 = dtstin
- den1 = rho - rhodit
-
- CALL rtheta(r1,th1,den1,y1)
-
- tt = th1*th1
- amu = aa*power(r1,beta*delta)*th1*(1.0d0-tt)
- y1 = dtstin + cc*amu
-
- CALL ss(r1,th1,s1,sd)
-
- rhoweg = xk1*power(r1,betai)*th1 + cc*s1(2)
- rho1s = den1 + cc*s1(1) + rhoweg
- error1 = rho - rhodi - rho1s
- r = r1
- th = th1
-
- IF (DABS(error1) .LT. 1.0d-5) THEN
- RETURN
- END IF
-
-*** rule for second pass ***
-
- den12 = rho - rhodi - cc*s1(1) + rhoweg
-
- IF (den12 .EQ. den1) den12 = den1 - 1.0d-6
-
- CALL rtheta(r1,th1,den12,y1)
-
- tt = th1*th1
- amu = aa*power(r1,beta*delta)*th1*(1.0d0-tt)
- y1 = dtstin + cc*amu
-
- CALL ss(r1,th1,s1,sd)
-
- rhoweg = xk1*power(r1,betai)*th1 + cc*s1(2)
- rho1s2 = den12 + cc*s1(1) + rhoweg
- error2 = rho - rhodi - rho1s2
-
- IF (DABS(error2) .LE. 1.0d-5) THEN
- r = r1
- th = th1
- error1 = error2
- rho1s = rho1s2
- RETURN
- END IF
-
-*** rule for nth pass ***
-
- den2 = den12
-
- DO 44 isig=1,10
- slope = (error2-error1)/(den2-den1)
- hold = den2
- den2 = den1 - (error1/slope)
-
- CALL rtheta(r1,th1,den2,y1)
-
- tt = th1*th1
- amu = aa*power(r1,beta*delta)*th1*(1.0d0-tt)
- y1 = dtstin + cc*amu
-
- CALL ss(r1,th1,s1,sd)
-
- rhoweg = xk1*power(r1,betai)*th1 + cc*s1(2)
- rho1s = den2 + cc*s1(1) + rhoweg
- error1 = error2
- error2 = rho - rhodi - rho1s
- r = r1
- th = th1
-
- IF (DABS(error2) .LT. 1.0d-6) RETURN
-
- den1 = hold
-
- 44 CONTINUE
-
- RETURN
- END
-
-*********************************************************************
-
-*** rtheta - Fits data for 1.0 < theta < 1.000001.
-* Solves:
-* rho = em*theta*(r**beta)
-* Tee = r*(1.0d0-besq*theta*theta)
-*
-* Routine given by Moldover (1978): Jour. Res. NBS, v. 84, n. 4,
-* p. 329 - 334.
-
-
- SUBROUTINE rtheta(r,theta,rho,Tee)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- COMMON /coefs/ a(20), q(20), x(11)
-
- SAVE
-
- EQUIVALENCE (beta,a(6)), (em,a(7)), (besq,a(9))
-
-
- IF (em .LE. 0.0d0 .OR. besq .LE. 1.0d0) GO TO 600
-
- absrho = DABS(rho)
-
- IF (absrho .LT. 1.0d-12) GO TO 600
-
- bee = DSQRT(besq)
-
- IF (DABS(Tee) .LT. 1.0d-12) GO TO 495
- IF (Tee .LT. 0.0d0) THEN
- z = 1.0d0-(1.0d0-bee)*Tee/(1.0d0-besq) *
- 1 power(em/absrho,1.0d0/beta)
- ELSE
- z = power(1.0d0+Tee*power(em/bee/absrho,1.0d0/beta),
- 1 -beta)
- END IF
- IF (z .GT. 1.00234d0*bee) GO TO 496
-
- c = -rho*bee/em/power(DABS(Tee),beta)
- z = DSIGN(z,rho)
-
- DO 500 n=1,16
- z2 = z*z
- z3 = 1.0d0 - z2
- dz = z3*(z+c*power(DABS(z3),beta))/(z3+2.0d0*beta*z2)
- z = z - dz
-
- IF (DABS(dz/z) .LT. 1.0d-12) GO TO 498
-
- 500 CONTINUE
-
- 601 IF (DABS(theta) .GT. 1.0001d0) theta = theta/DABS(theta)
- RETURN
-
- 498 theta = z/bee
- r = Tee/(1.0d0-z*z)
- r = DABS(r)
- RETURN
-
- 495 theta = DSIGN(1.0d0,rho)/bee
- r = power(rho/(em*theta),1.0d0/beta)
- RETURN
-
- 496 theta = DSIGN(1.0d0,rho)
- r = Tee/(1.0d0-besq)
- r = DABS(r)
- RETURN
-
- 600 IF (DABS(Tee) .LT. 1.0d-12) GO TO 601
-
- IF (Tee .LT. 0.0d0) GO TO 496
-
- theta = 1.0d-12
- r = Tee
- RETURN
-
- END
-
-*********************************************************************
-
-*** ss - Computes terms of the summation that defines dPotl/dT
-* and the 1st derivative of the theta (s) square polynomial.
-
- SUBROUTINE ss(r,th,s,sd)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- DOUBLE PRECISION s(2), sd(2), sx(2)
-
- COMMON /coefs/ a(20), q(20), x(11)
-***************************************************************
- COMMON /abc1/ dPdM
-***************************************************************
-
- SAVE
-
- EQUIVALENCE (alpha,q(10)), (beta,a(6)), (besq,a(9)),
- 1 (delta,a(11)), (deli,q(14)), (alhi,q(15)),
- 2 (beti, q(16)), (gami,q(17)), (p00, q(11)),
- 3 (p01, q(18)), (s00, a(17)), (s20, a(18)),
- 4 (s01, a(19)), (s21, a(20))
-
- tt = th*th
- sx(1) = s00 + s20*tt
- sd(1) = 2.0d0*s20*th
- sx(2) = s01 + s21*tt
- sd(2) = 2.0d0*s21*th
- s(1) = sx(1)*a(10)*a(7)*power(r,1.0d0-alpha)
- s(2) = sx(2)*a(10)*a(12)*power(r,1.0d0-alhi)
-
- dPdM = power(r,beta)*a(7)*th + a(1)*power(r,1.0d0-alpha)*
- 1 a(10)*a(7)*sx(1) +
- 2 power(r,beti)*a(12)*th + a(1)*power(r,1.0d0-alhi)*
- 3 a(10)*a(12)*sx(2)
-
- RETURN
- END
-
-*****************************************************************
-
-*** thmLVS - Calculates thermodynamic and transport properties
-* of critical region H2O using the Levelt Sengers, et al
-* (1983) equation of state.
-
- SUBROUTINE thmLVS(isat,T,r1,th1)
-
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-
- DOUBLE PRECISION s(2), xk(2), sd(2)
-
- COMMON /coefs/ a(20), q(20), x(11)
- COMMON /crits/ Tc, rhoc, Pc, Pcon, Ucon, Scon, dPcon
- COMMON /therm/ AE, GE, U, H, Entrop, Cp, Cv, betaw, alphw,
- 1 heat, Speed
- COMMON /satur/ Dliq, Dvap, DH2O, iphase
- COMMON /deriv/ amu, s, Pw, Tw, dTw, dM0dT, dP0dT,
- 1 d2PdM2, d2PdMT, d2PdT2, p0th, p1th, xk
- COMMON /deri2/ dPdD, dPdT
-*************************************************************
- COMMON /abc1/ dPdM
- COMMON /abc3/ dPdTcd
-*************************************************************
-
- SAVE
-
- EQUIVALENCE (pw2, a(4)), (pw3, a(2)), (besq, a(9)),
- 1 (amc, a(13)), (am1, a(14)), (am2, a(15)),
- 2 (aa, a(10)), (xk0, a(7)), (am3, a(16)),
- 3 (xk1, a(12)), (pw11,q(9)), (alpha, q(10)),
- 4 (alhi,q(15)), (pw1, a(5))
-
- d2P0dT = 2.0d0*pw2 + 6.0d0*pw3*dTw
- d2M0dT = 2.0d0*am2 + 6.0d0*am3*dTw
-
- dP0dT = pw1+dTw*(2.0d0*pw2+3.0d0*pw3*dTw)
- dM0dT = am1+dTw*(2.0d0*am2+3.0d0*am3*dTw)
-
- IF (isat .EQ. 0) THEN
- rho = DH2O / rhoc
- Uw = dP0dT-rho*dM0dT+pw11*amu+s(1)+s(2)
- ELSE
- rho = th1 * (xk0*power(r1,a(6)) + xk1*power(r1,q(16)))
- 1 + a(1)*(s(1)+s(2))
- rho = 1.0d0+pw11*dTw+rho
- Uw = dP0dT-rho*dM0dT+pw11*amu+s(1)+s(2)
- DH2O = rho * rhoc
- dPdT2 = Pw - Tw*(Uw+rho*dM0dT)
- heat = 1.0d3*T*(Pcon*dPdT2)*(1.0d0/Dvap-1.0d0/Dliq)
-
- CALL ss(r1,th1,s,sd)
- CALL aux(r1,th1,d2PdT2,d2PdMT,d2PdM2,aa,xk,sd,Cvcoex)
- IF (r1 .NE. 0.0d0) THEN
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 951
More information about the CHNOSZ-commits
mailing list