[CHNOSZ-commits] r708 - in pkg/CHNOSZ: . inst/extdata inst/extdata/src man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 23 07:40:09 CET 2022
Author: jedick
Date: 2022-03-23 07:40:09 +0100 (Wed, 23 Mar 2022)
New Revision: 708
Added:
pkg/CHNOSZ/inst/extdata/src/
pkg/CHNOSZ/inst/extdata/src/H2O92D.f.orig
pkg/CHNOSZ/inst/extdata/src/README.txt
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/man/ionize.aa.Rd
Log:
Add extdata/src/H2O92D.f.orig
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-02-20 09:03:54 UTC (rev 707)
+++ pkg/CHNOSZ/DESCRIPTION 2022-03-23 06:40:09 UTC (rev 708)
@@ -1,6 +1,6 @@
-Date: 2022-02-20
+Date: 2022-03-23
Package: CHNOSZ
-Version: 1.4.3
+Version: 1.4.3-1
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Added: pkg/CHNOSZ/inst/extdata/src/H2O92D.f.orig
===================================================================
--- pkg/CHNOSZ/inst/extdata/src/H2O92D.f.orig (rev 0)
+++ pkg/CHNOSZ/inst/extdata/src/H2O92D.f.orig 2022-03-23 06:40:09 UTC (rev 708)
@@ -0,0 +1,3457 @@
+*** H2O92 - Computes state, thermodynamic, transport, and electroststic
+*** properties of fluid H2O at T,[P,D] using equations and data
+*** given by Haar et al. (1984), Levelt Sengers et al. (1983),
+*** Johnson and Norton (1991), Watson et al. (1980), Sengers and
+*** Kamgar-Parsi (1984), Sengers et al. (1984), Helgeson and Kirkham
+*** (1974), Uematsu and Franck (1980), and Pitzer (1983).
+***
+***********************************************************************
+***
+*** Author: James W. Johnson
+*** Earth Sciences Dept., L-219
+*** Lawrence Livermore National Laboratory
+*** Livermore, CA 94550
+*** johnson at s05.es.llnl.gov
+***
+*** Abandoned: 8 November 1991
+***
+***********************************************************************
+*
+* specs - Input unit, triple point, saturation, and option specs:
+*
+***** it, id, ip, ih, itripl, isat, iopt, useLVS, epseqn, icrit;
+*
+* note that the returned value of isat may differ from
+* its input value and that icrit need not be specified
+* prior to invocation.
+*
+*
+* states - State variables:
+*
+***** temp, pres, dens(1), dens(2);
+*
+* note that the first three of these must be specified prior
+* to invocation and that, in the case of saturation, vapor
+* density is returned in dens(1), liquid in dens(2).
+*
+*
+* props - Thermodynamic, transport, electrostatic, and combined
+* property values:
+*
+***** A, G, S, U, H, Cv, Cp, Speed, alpha, beta, diel, visc,
+***** tcond, surten, tdiff, Prndtl, visck, albe,
+***** ZBorn, YBorn, QBorn, daldT, XBorn
+*
+*
+* error - LOGICAL argument that indicates success ("FALSE") or
+* failure ("TRUE") of the call, the latter value in
+* response to out-of-bounds specs or states variables.
+*
+***********************************************************************
+
+ SUBROUTINE H2O92(specs,states,props,error)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ PARAMETER (NPROP = 23, NPROP2 = 46)
+
+ INTEGER specs(10)
+ DOUBLE PRECISION states(4), props(NPROP2), Dens(2),
+ 1 wpliq(NPROP), wprops(NPROP)
+ LOGICAL crtreg, valid, error, useLVS
+
+ COMMON /units/ ft, fd, fvd, fvk, fs, fp, fh, fst, fc
+ COMMON /wpvals/ wprops, wpliq
+
+ SAVE
+
+
+ CALL unit(specs(1),specs(2),specs(3),specs(4),specs(5))
+
+ IF (.NOT. (valid(specs(1),specs(2),specs(3),specs(4),specs(5),
+ 1 specs(6),specs(7),specs(8),specs(9),
+ 2 states(1),states(2),states(3)))) THEN
+ error = .TRUE.
+ RETURN
+ ELSE
+ error = .FALSE.
+ END IF
+
+ IF (crtreg(specs(6),specs(7),specs(1),
+ 1 states(1),states(2),states(3))) THEN
+ specs(10) = 1
+ useLVS = (specs(8) .EQ. 1)
+ ELSE
+ specs(10) = 0
+ useLVS = .FALSE.
+ 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
+ Dens(1) = states(3) / 1.0d3
+ CALL HGKeqn(specs(6),specs(7),specs(5),
+ 1 states(1),states(2),Dens,specs(9))
+ END IF
+
+ CALL load(1,wprops,props)
+
+ IF (specs(6) .EQ. 1) THEN
+ tempy = Dens(1)
+ Dens(1) = Dens(2)
+ Dens(2) = tempy
+ CALL load(2,wpliq,props)
+ END IF
+
+ states(1) = TdegUS(specs(1),states(1))
+ states(2) = states(2) * fp
+ states(3) = Dens(1) / fd
+
+ IF (specs(6) .EQ. 1) THEN
+ states(4) = Dens(2) / fd
+ END IF
+
+ RETURN
+ END
+
+************************************************************************
+
+*** valid - Returns "TRUE" if unit and equation specifications
+* are valid and input state conditions fall within
+* the HGK equation's region of validity;
+* returns "FALSE" otherwise.
+
+ LOGICAL FUNCTION valid(it,id,ip,ih,itripl,isat,iopt,
+ 1 useLVS,epseqn,Temp,Pres,Dens)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ INTEGER useLVS, epseqn
+ LOGICAL valspc, valTD, valTP
+
+ COMMON /tolers/ TTOL, PTOL, DTOL, XTOL, EXPTOL, FPTOL
+ COMMON /units/ ft, fd, fvd, fvk, fs, fp, fh, fst, fc
+ 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
+
+ SAVE
+
+
+*** ensure validity of input specifications
+ IF (.NOT. valspc(it,id,ip,ih,itripl,isat,iopt,
+ 1 useLVS,epseqn)) THEN
+ valid = .FALSE.
+ RETURN
+ END IF
+
+*** convert to degC, bars, g/cm3 ***
+ T = TdegK(it,Temp) - 273.15d0
+ D = Dens * fd
+ P = Pres / fp * 1.0d1
+ Ttripl = Ttr - 273.15d0
+ Tcrit = Tc - 273.15d0
+ Pcrit = Pc * 1.0d1
+
+ IF (isat .EQ. 0) THEN
+ IF (iopt .EQ. 1) THEN
+ valid = valTD(T,D,isat,epseqn)
+ ELSE
+ valid = valTP(T,P)
+ END IF
+ ELSE
+ IF (iopt .EQ. 1) THEN
+ valid = ((T+FPTOL .GE. Ttripl) .AND.
+ 1 (T-FPTOL .LE. Tcrit))
+ ELSE
+ valid = ((P+FPTOL .GE. Ptripl) .AND.
+ 1 (P-FPTOL .LE. Pcrit))
+ END IF
+ END IF
+
+ RETURN
+ END
+
+*****************************************************************
+
+*** valspc - Returns "TRUE" if it, id, ip, ih, itripl, isat, iopt,
+* useLVS, and epseqn values all define valid input;
+* returns "FALSE" otherwise.
+
+ LOGICAL FUNCTION valspc(it,id,ip,ih,itripl,isat,iopt,
+ 1 useLVS,epseqn)
+
+ INTEGER useLVS, epseqn
+
+ SAVE
+
+
+ valspc = (1 .LE. it) .AND. (it .LE. 4) .AND.
+ 1 (1 .LE. id) .AND. (id .LE. 4) .AND.
+ 2 (1 .LE. ip) .AND. (ip .LE. 5) .AND.
+ 3 (1 .LE. ih) .AND. (ih .LE. 6) .AND.
+ 4 (0 .LE. itripl) .AND. (itripl .LE. 1) .AND.
+ 5 (0 .LE. isat) .AND. (isat .LE. 1) .AND.
+ 6 (1 .LE. iopt) .AND. (iopt .LE. 2) .AND.
+ 7 (0 .LE. useLVS) .AND. (useLVS .LE. 1) .AND.
+ 8 (1 .LE. epseqn) .AND. (epseqn .LE. 5)
+
+ RETURN
+ END
+
+*****************************************************************
+
+*** 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.
+
+ LOGICAL FUNCTION valTP(T,P)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ COMMON /tolers/ TTOL, PTOL, DTOL, XTOL, EXPTOL, FPTOL
+ COMMON /crits/ Tcrit, 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
+
+ SAVE
+
+
+ IF ((T-FPTOL .GT. Ttop) .OR. (T+FPTOL .LT. Tbtm) .OR.
+ 1 (P-FPTOL .GT. Ptop) .OR. (P+FPTOL .LT. Pbtm)) THEN
+ valTP = .FALSE.
+ RETURN
+ ELSE
+ valTP = .TRUE.
+ END IF
+
+ IF (P .GE. Pli13) THEN
+ Plimit = sPli37 * (T-Tli13) + Pli13
+ valTP = (P-FPTOL .LE. Plimit)
+ ELSE
+ IF (P .GE. Ptripl) THEN
+ Plimit = sPli1 * (T-Tli13) + Pli13
+ valTP = (P+FPTOL .GE. Plimit)
+ ELSE
+ Psubl = Psublm(T)
+ valTP = (P-FPTOL .LE. Psubl)
+ END IF
+ END IF
+
+ RETURN
+ END
+
+*****************************************************************
+
+*** 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
+* g1, g2, gf = alpha, beta, gamma from eq (A.2), p.272
+* g, ii, jj = g(i), k(i), l(i) from eq (A.5), p.272.
+* Note that tz < tcHGK.
+* Tolerence limits required in various real & inexact
+* comparisons are set and stored in COMMON /tolers/.
+
+
+ BLOCK DATA HGKcon
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ COMMON /aconst/ wm, gascon, tz, aa, zb, dzb, yb, uref, sref
+ COMMON /nconst/ g(40), ii(40), jj(40), nc
+ COMMON /ellcon/ g1, g2, gf, b1, b2, b1t, b2t, b1tt, b2tt
+ COMMON /bconst/ bp(10), bq(10)
+ COMMON /addcon/ atz(4), adz(4), aat(4), aad(4)
+ COMMON /HGKcrt/ tcHGK, dcHGK, pcHGK
+ COMMON /tolers/ TTOL, PTOL, DTOL, XTOL, EXPTOL, FPTOL
+ COMMON /HGKbnd/ Ttop, Tbtm, Ptop, Pbtm, Dtop, Dbtm
+ COMMON /liqice/ sDli1, sPli1, sDli37, sPli37, sDIB30,
+ 1 Tli13, Pli13, Dli13, TnIB30, DnIB30
+ COMMON /tpoint/ Utripl, Stripl, Htripl, Atripl, Gtripl,
+ 1 Ttripl, Ptripl, Dltrip, Dvtrip
+
+ SAVE
+
+ DATA Ttripl, Ptripl, Dltrip, Dvtrip
+ 1 / 2.7316d2,
+ 2 0.611731677193563186622762580414d-2,
+ 3 0.999778211030936587977889295063d0,
+ 4 0.485467583448287303988319166423d-5 /
+
+ DATA Ttop, Tbtm, Ptop, Pbtm, Dtop, Dbtm
+ 1 / 2.25d3, -2.0d1, 3.0d4, 1.0d-3,
+ 2 0.138074666423686955066817336896d1,
+ 3 0.858745555396173972667420987465d-7 /
+
+ DATA sDli1, sPli1, sDli37, sPli37, sDIB30,
+ 1 Tli13, Pli13, Dli13, TnIB30, DnIB30
+ 2 / -0.584797401732178547634910059828d-2,
+ 3 -0.138180804975562958027981345769d3,
+ 4 0.183244000000000000000000000007d-2,
+ 5 0.174536874999999999999999999995d3,
+ 6 -0.168375439429928741092636579574d-3,
+ 7 -0.15d2,
+ 8 0.20741d4,
+ 9 0.108755631570602617113573577945d1,
+ 1 0.145d3,
+ 1 0.102631640581853166397515716306d1 /
+
+
+ DATA TTOL, PTOL, DTOL, XTOL, EXPTOL, FPTOL
+ 1 / 1.0d-6, 1.0d-6, 1.0d-9, 1.0d-5, -673.5d0, 1.0d-7 /
+
+ DATA tcHGK, dcHGK, pcHGK / .647126d3, .322d3, .22055d2 /
+
+ DATA atz /.64d3, .64d3, .6416d3, .27d3/
+ DATA adz /.319d0, .319d0, .319d0, .155d1/
+ DATA aat /.2d5, .2d5, .4d5, .25d2/
+ DATA aad /.34d2, .4d2, .3d2, .105d4/
+
+ DATA wm, gascon, tz, aa, uref, sref
+ 1 / .1801520000d2, .46152200d0, .647073d3, .1d1,
+ 2 -.4328455039d4, .76180802d1 /
+
+ DATA g1, g2, gf /.11d2, .44333333333333d2, .35d1/
+
+ DATA bp / .7478629d0, -.3540782d0, 2*.0d0, .7159876d-2,
+ 1 .0d0, -.3528426d-2, 3*.0d0/
+
+ DATA bq / .11278334d1, .0d0, -.5944001d0, -.5010996d1, .0d0,
+ 1 .63684256d0, 4*.0d0/
+
+ DATA nc / 36 /
+
+ DATA g /-.53062968529023d3, .22744901424408d4, .78779333020687d3
+ 1, -.69830527374994d2, .17863832875422d5, -.39514731563338d5
+ 2, .33803884280753d5, -.13855050202703d5, -.25637436613260d6
+ 3, .48212575981415d6, -.34183016969660d6, .12223156417448d6
+ 4, .11797433655832d7, -.21734810110373d7, .10829952168620d7
+ 5, -.25441998064049d6, -.31377774947767d7, .52911910757704d7
+ 6, -.13802577177877d7, -.25109914369001d6, .46561826115608d7
+ 7, -.72752773275387d7, .41774246148294d6, .14016358244614d7
+ 8, -.31555231392127d7, .47929666384584d7, .40912664781209d6
+ 9, -.13626369388386d7, .69625220862664d6, -.10834900096447d7
+ a, -.22722827401688d6, .38365486000660d6, .68833257944332d4
+ b, .21757245522644d5, -.26627944829770d4, -.70730418082074d5
+ c, -.22500000000000d0, -.16800000000000d1
+ d, .5500000000000d-1, -.93000000000000d2/
+
+ DATA ii / 4*0, 4*1, 4*2, 4*3, 4*4, 4*5, 4*6, 4*8, 2*2, 0, 4,
+ 1 3*2, 4/
+
+ DATA jj / 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7,
+ 1 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7,
+ 2 1, 4, 4, 4, 0, 2, 0, 0/
+
+ END
+
+*********************************************************************
+
+*** LVScon - Constant parameters for the H2O critical region equation
+* of state given by Levelt Sengers, Kamgar-Parsi, Balfour,
+* & Sengers (1983).
+
+ BLOCK DATA LVScon
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ COMMON /crits/ Tc, rhoC, Pc, Pcon, Ucon, Scon, dPcon
+ COMMON /coefs/ a(20), q(20), x(11)
+
+ SAVE
+
+ DATA Tc, rhoC, Pc, Pcon, Ucon, Scon, dPcon
+ 1 / 647.067d0, 322.778d0, 22.046d0,
+ 2 0.034070660379837018423130834983d0, 22046.0d0,
+ 3 0.034070660379837018423130834983d3,
+ 4 0.000000327018783663660700780197d0 /
+
+ DATA a / -0.017762d0, 5.238000d0, 0.000000d0, -2.549150d1,
+ 1 6.844500d0, 0.325000d0, 1.440300d0, 0.000000d0,
+ 2 1.375700d0, 2.366660d1, 4.820000d0, 0.294200d0,
+ 3 -1.123260d1, -2.265470d1, -1.788760d1, -4.933200d0,
+ 4 1.109430391161019373812391218008d0,
+ 5 -1.981395981400671095301629432211d0,
+ 6 0.246912528778663959151808173743d0,
+ 7 -0.843411332867484343974055795059d0 /
+
+ DATA q / -0.006000d0, -0.003000d0, 0.000000d0, 6.470670d2,
+ 1 3.227780d2, 2.204600d1, 0.267000d0, -1.600000d0,
+ 2 0.491775937675717720291497417773d0, 0.108500d0,
+ 3 0.586534703230779473334597524774d0,
+ 4 -1.026243389120214352553706598564d0,
+ 5 0.612903225806451612903225804745d0, 0.500000d0,
+ 6 -0.391500d0, 0.825000d0, 0.741500d0,
+ 7 0.103245882826119154987166286332d0,
+ 8 0.160322434159191991394857495360d0,
+ 9 -0.169859514687100893997445721324d0 /
+
+ DATA x / 6.430000d2, 6.453000d2, 6.950000d2,
+ 1 1.997750d2, 4.200400d2,
+ 2 2.09945691135940719075293945960d1,
+ 3 2.15814057875264119875397458907d1,
+ 4 3.0135d1, 4.0484d1,
+ 5 .175777517046267847932127026995d0,
+ 6 .380293646126229135059562456934d0 /
+
+
+* EQUIVALENCE (cc, a(1) ), (pointA, q(1) ), (Tmin1, x(1)),
+* 1 (p3, a(2) ), (pointB, q(2) ), (Tmin2, x(2)),
+* 2 (delroc, a(3) ), (delpc, q(3) ), (Tmax, x(3)),
+* 3 (p2, a(4) ), (Tc, q(4) ), (Dmin, x(4)),
+* 4 (p1, a(5) ), (rhoc, q(5) ), (Dmax, x(5)),
+* 5 (beta, a(6) ), (Pc, q(6) ), (Pmin1, x(6)),
+* 6 (xko, a(7) ), (dPcdTc, q(7) ), (Pmin2, x(7)),
+* 7 (delTc, a(8) ), (slopdi, q(8) ), (Pmax1, x(8)),
+* 8 (besq, a(9) ), (p11, q(9) ), (Pmax2, x(9)),
+* 9 (aa, a(10)), (alpha, q(10)), (sl1, x(10)),
+* 0 (delta, a(11)), (p00, q(11)), (sl2, x(11)),
+* 1 (k1, a(12)), (p20, q(12)),
+* 2 (muc, a(13)), (p40, q(13)),
+* 3 (mu1, a(14)), (deli, q(14)),
+* 4 (mu2, a(15)), (alh1, q(15)),
+* 5 (mu3, a(16)), (beti, q(16)),
+* 6 (s00, a(17)), (gami, q(17)),
+* 7 (s20, a(18)), (p01, q(18)),
+* 8 (s01, a(19)), (p21, q(19)),
+* 9 (s21, a(20)), (p41, q(20))
+
+ END
+
+*******************************************************************
+
+*** unit - Sets internal parameters according to user-specified
+* choice of units. Internal program units are degK(T),
+* and gm/cm**3(D); all other properties are computed in
+* dimensionless form and dimensioned at output time.
+* NOTE: conversion factors for j/g ---> cal/(g,mole)
+* (ffh (4 & 5)) are consistent with those given in
+* Table 1, Helgeson & Kirkham (1974a) for thermal calories,
+* and differ slightly with those given by Haar et al (1984)
+* for international calories.
+
+ SUBROUTINE unit(it,id,ip,ih,itripl)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ DOUBLE PRECISION fft(4), ffd(4), ffvd(4), ffvk(4),
+ 1 ffs(4), ffp(5), ffh(6),
+ 2 ffst(4), ffcd(4), ffch(6)
+
+ COMMON /units/ ft, fd, fvd, fvk, fs, fp, fh, fst, fc
+
+ SAVE
+
+ DATA fft /1.0d0, 1.0d0, 0.555555556d0, 0.555555556d0 /
+ DATA ffd /1.0d-3, 1.0d0, 1.80152d-2, 1.6018d-2/
+ DATA ffvd /1.0d0, 1.0d1, 0.555086816d0, 0.671968969d0 /
+ DATA ffvk /1.0d0, 1.0d4, 1.0d4, 1.076391042d1 /
+ DATA ffs /1.0d0, 1.0d2, 1.0d2, 3.280833d0 /
+ DATA ffp /1.0d0, 1.0d1, 9.869232667d0, 1.45038d2, 1.01971d1/
+ DATA ffh /1.0d0, 1.0d0, 1.80152d1, 2.3901d-1,
+ 1 4.305816d0, 4.299226d-1/
+ DATA ffst /1.0d0, 1.0d3, 0.555086816d2, 0.2205061d1 /
+ DATA ffcd /1.0d0, 1.0d-2, 1.0d-2, 0.3048d0 /
+ DATA ffch /1.0d-3, 1.0d0, 1.0d0, 0.23901d0,
+ 1 0.23901d0, 0.947244d-3 /
+
+
+ ft = fft(it)
+ fd = ffd(id)
+ fvd = ffvd(id)
+ fvk = ffvk(id)
+ fs = ffs(id)
+ fp = ffp(ip)
+ fh = ffh(ih)
+ fst = ffst(id)
+ fc = ffcd(id) * ffch(ih)
+
+ IF (itripl .EQ. 1) CALL tpset
+
+ RETURN
+
+ END
+
+***********************************************************************
+
+*** crtreg - Returns "TRUE" if input state conditions fall within
+* the critical region of H2O; otherwise returns "FALSE".
+* T, P, D, input in user-specified units, are returned in
+* degK, MPa, kg/m3.
+
+ LOGICAL FUNCTION crtreg(isat,iopt,it,T,P,D)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ LOGICAL llim, ulim
+
+ COMMON /crits/ Tc, rhoc, Pc, Pcon, Ucon, Scon, dPcon
+ COMMON /coefs/ a(20), q(20), x(11)
+ COMMON /units/ ft, fd, fvd, fvk, fs, fp, fh, fst, fc
+
+ SAVE
+
+ EQUIVALENCE (Tmin1, x(1)), (Tmin2, x(2)), (Tmax, x(3)),
+ 1 (Dmin, x(4)), (Dmax, x(5)),
+ 2 (Pbase1, x(6)), (Pbase2,x(7)),
+ 3 (PTmins, x(10)), (PTmaxs,x(11))
+
+
+ T = TdegK(it,T)
+ IF (isat .EQ. 0) THEN
+ IF (iopt .EQ. 1) THEN
+ D = D * fd * 1.0d3
+ crtreg = ((T .GE. Tmin1) .AND. (T .LE. Tmax) .AND.
+ 1 (D .GE. Dmin) .AND. (D .LE. Dmax))
+ ELSE
+ P = P / fp
+ IF ((T .LT. Tmin1) .OR. (T .GT. Tmax)) THEN
+ crtreg = .FALSE.
+ ELSE
+ Pmin = Pbase1 + PTmins * (T - Tmin1)
+ Pmax = Pbase2 + PTmaxs * (T - Tmin2)
+ llim = (P .GE. Pmin)
+ ulim = (P .LE. Pmax)
+ 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
+ END IF
+ END IF
+ END IF
+ ELSE
+ IF (iopt .EQ. 1) THEN
+ crtreg = (T .GE. Tmin1)
+ ELSE
+ P = P / fp
+ crtreg = (P .GE. Pbase1)
+ END IF
+ END IF
+
+ RETURN
+ END
+
+*********************************************************************
+
+*** HGKeqn - Computes thermodynamic and transport properties of
+* of H2O from the equation of state given by
+* Haar, Gallagher, & Kell (1984).
+
+ SUBROUTINE HGKeqn(isat,iopt,itripl,Temp,Pres,Dens,epseqn)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ PARAMETER (NPROP = 23)
+
+ INTEGER epseqn
+ DOUBLE PRECISION Dens(2), wprops(NPROP), wpliq(NPROP)
+
+ COMMON /aconst/ wm, gascon, tz, aa, zb, dzb, yb, uref, sref
+ COMMON /wpvals/ wprops, wpliq
+ COMMON /RTcurr/ rt
+
+ SAVE
+
+ rt = gascon * Temp
+
+ CALL HGKsat(isat,iopt,itripl,Temp,Pres,Dens,epseqn)
+
+ IF (isat .EQ. 0) THEN
+ CALL bb(Temp)
+ CALL calcv3(iopt,itripl,Temp,Pres,Dens(1),epseqn)
+ CALL thmHGK(Dens(1),Temp)
+ CALL dimHGK(isat,itripl,Temp,Pres,Dens(1),epseqn)
+ ELSE
+ DO 10 i=1,NPROP
+ 10 wpliq(i) = wprops(i)
+ CALL dimHGK(2,itripl,Temp,Pres,Dens(2),epseqn)
+ END IF
+
+ RETURN
+ END
+
+*****************************************************************
+
+*** HGKsat - If isat=1, computes Psat(T) or Tsat(P) (iopt=1,2),
+* liquid and vapor densities, and associated
+* thermodynamic and transport properties.
+* If isat=0, checks whether T-D or T-P (iopt=1,2)
+* falls on or within TOL of the liquid-vapor
+* surface; if so, sets isat <- 1 and computes
+* properties.
+
+ SUBROUTINE HGKsat(isat,iopt,itripl,Temp,Pres,Dens,epseqn)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ DOUBLE PRECISION Dens(2)
+ INTEGER epseqn
+
+ COMMON /tolers/ TTOL, PTOL, DTOL, XTOL, EXPTOL, FPTOL
+ COMMON /HGKcrt/ tcHGK, dcHGK, pcHGK
+ COMMON /aconst/ wm, gascon, tz, aa, zb, dzb, yb, uref, sref
+ COMMON /units/ ft, fd, fvd, fvk, fs, fp, fh, fst, fc
+ COMMON /tpoint/ Utr, Str, Htr, Atr, Gtr,
+ 1 Ttripl, Ptripl, Dltrip, Dvtrip
+ COMMON /crits/ Tc, rhoC, Pc, Pcon, Ucon, Scon, dPcon
+
+ SAVE
+
+
+ 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.
+ 1 ((iopt .EQ. 2) .AND. (Pres .GT. Pc))) THEN
+ RETURN
+ ELSE
+ CALL pcorr(itripl,Temp,Ptemp,dltemp,dvtemp,epseqn)
+ IF (((iopt .EQ. 2) .AND.
+ 1 (DABS(Pres-Ptemp) .LE. PTOL)) .OR.
+ 2 ((iopt .EQ. 1) .AND.
+ 3 ((DABS(Dens(1)-dltemp) .LE. DTOL) .OR.
+ 4 (DABS(Dens(1)-dvtemp) .LE. DTOL)))) THEN
+ isat = 1
+ Pres = Ptemp
+ Dens(1) = dltemp
+ Dens(2) = dvtemp
+ END IF
+ END IF
+ END IF
+
+ RETURN
+ END
+
+************************************************************************
+
+*** calcv3 - Compute the dependent state variable.
+
+ SUBROUTINE calcv3(iopt,itripl,Temp,Pres,Dens,epseqn)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ INTEGER epseqn
+
+ COMMON /units/ ft, fd, fvd, fvk, fs, fp, fh, fst, fc
+ COMMON /qqqq/ q0, q5
+ COMMON /aconst/ wm, gascon, tz, aa, z, dz, y, uref, sref
+ COMMON /fcts/ ad, gd, sd, ud, hd, cvd, cpd, dpdt, dvdt, dpdd,
+ 1 cjtt, cjth
+ COMMON /RTcurr/ rt
+
+ SAVE
+
+
+ IF (iopt .EQ. 1) THEN
+ CALL resid(Temp,Dens)
+ CALL base(Dens,Temp)
+ CALL ideal(Temp)
+ Pres = rt * Dens * z + q0
+ ELSE
+ IF (Temp .LT. tz) THEN
+ CALL pcorr(itripl,Temp,ps,dll,dvv,epseqn)
+ ELSE
+ ps = 2.0d4
+ dll = 0.0d0
+ END IF
+ IF (Pres .GT. ps) THEN
+ dguess = dll
+ ELSE
+ dguess = Pres / Temp / 0.4d0
+ END IF
+
+ CALL denHGK(Dens,Pres,dguess,Temp,dpdd)
+ CALL ideal(Temp)
+ END IF
+
+ RETURN
+ END
+
+******************************************************************************
+
+*** thmHGK - Computes thermodynamic functions in dimensionless
+* units from the HGK equation of state: Helmholtz, Gibbs,
+* internal energy, and enthalpy functions (ad, gd, ud, hd) are
+* per RT; entropy and heat capacities (sd, cvd, cpd) are per R.
+
+ SUBROUTINE thmHGK(d,t)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ COMMON /aconst/ wm, gascon, tz, aa, zb, dzb, y, uref, sref
+ COMMON /qqqq/ qp, qdp
+ COMMON /basef/ ab, gb, sb, ub, hb, cvb, pb, dpdtb
+ COMMON /resf/ ar, gr, sr, ur, hr, cvr, dpdtr
+ COMMON /idf/ ai, gi, si, ui, hi, cvi, cpi
+ COMMON /fcts/ ad, gd, sd, ud, hd, cvd, cpd, dpdt, dvdt, dpdd,
+ 1 cjtt, cjth
+ COMMON /RTcurr/ rt
+
+ SAVE
+
+
+ z = zb + qp/rt/d
+ dpdd = rt * (zb + y * dzb) + qdp
+ ad = ab + ar + ai - uref/t + sref
+ gd = ad + z
+ ud = ub + ur + ui - uref/t
+ dpdt = rt * d * dpdtb + dpdtr
+ cvd = cvb + cvr + cvi
+ cpd = cvd + t*dpdt*dpdt/(d*d*dpdd*gascon)
+ hd = ud + z
+ sd = sb + sr + si - sref
+ dvdt = dpdt / dpdd / d / d
+ cjtt = 1.0d0 / d - t * dvdt
+ cjth = -cjtt / cpd / gascon
+
+ RETURN
+
+ END
+
+*************************************************************************
+
+*** bb - Computes molecular parameters b, the "excluded volume"
+* (eq A.3), and B, the second virial coefficient (eq A.4),
+* in cm3/g (b1,b2) and their first and second derivatives
+* with respect to temperature (b1t,b1tt,b2t,b2tt).
+
+ SUBROUTINE bb(t)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ DOUBLE PRECISION v(10)
+
+ COMMON /ellcon/ g1, g2, gf, b1, b2, b1t, b2t, b1tt, b2tt
+ COMMON /aconst/ wm, gascon, tz, aa, z, dz, y, uref, sref
+ COMMON /bconst/ bp(10), bq(10)
+
+ SAVE
+
+
+ v(1) = 1.0d0
+
+ DO 2 i=2,10
+ 2 v(i) = v(i-1) * tz / t
+
+ b1 = bp(1) + bp(2) * DLOG(1.0 / v(2))
+ b2 = bq(1)
+ b1t = bp(2) * v(2) / tz
+ b2t = 0.0d0
+ b1tt = 0.0d0
+ b2tt = 0.0d0
+
+ DO 4 i=3,10
+ b1 = b1 + bp(i) * v(i-1)
+ b2 = b2 + bq(i) * v(i-1)
+ b1t = b1t - (i-2) * bp(i) * v(i-1) / t
+ b2t = b2t - (i-2) * bq(i) * v(i-1) / t
+ b1tt = b1tt + bp(i) * (i-2)*(i-2) * v(i-1) / t / t
+ 4 b2tt = b2tt + bq(i) * (i-2)*(i-2) * v(i-1) / t / t
+
+ b1tt = b1tt - b1t / t
+ b2tt = b2tt - b2t / t
+
+ RETURN
+
+ END
+
+***********************************************************************
+
+*** base - Computes Abase, Gbase, Sbase, Ubase, Hbase, Cvbase
+* -- all per RT (dimensionless) -- as well as Pbase & dP/dT
+* -- both per (DRT) -- for the base function (ab, gb, sb, ub,
+* hb, cvb, pb, dpdtb). See Haar, Gallagher & Kell (1979), eq(1).
+
+
+ SUBROUTINE base(d,t)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ COMMON /ellcon/ g1, g2, gf, b1, b2, b1t, b2t, b1tt, b2tt
+ COMMON /basef/ ab, gb, sb, ub, hb, cvb, pb, dpdtb
+ COMMON /aconst/ wm, gascon, tz, a, z, dz, y, uref, sref
+
+ SAVE
+
+
+ y = .25d0 * b1 * d
+ x = 1.0d0 - y
+ z0 = (1.0d0 + g1*y + g2*y*y) / (x*x*x)
+ z = z0 + 4.0d0*y*(b2/b1 - gf)
+ dz0 = (g1 + 2.0d0*g2*y)/(x*x*x) +
+ 1 3.0d0*(1.0d0 + g1*y + g2*y*y)/(x*x*x*x)
+ dz = dz0 + 4.0d0*(b2/b1 - gf)
+
+ pb = z
+
+ ab = -DLOG(x) - (g2 - 1.0d0)/x + 28.16666667d0/x/x +
+ 1 4.0d0*y*(b2/b1 - gf) + 15.166666667d0 +
+ 2 DLOG(d*t*gascon/.101325d0)
+ gb = ab + z
+ ub = -t*b1t*(z - 1.0d0 - d*b2)/b1 - d*t*b2t
+ sb = ub - ab
+ hb = z + ub
+
+ bb2tt = t * t * b2tt
+ cvb = 2.0d0*ub + (z0 - 1.0d0)*(((t*b1t/b1)*(t*b1t/b1)) -
+ 1 t*t*b1tt/b1) - d*(bb2tt - gf*b1tt*t*t) -
+ 2 (t*b1t/b1)*(t*b1t/b1)*y*dz0
+
+ dpdtb = pb/t + d*(dz*b1t/4.0d0 + b2t - b2/b1*b1t)
+
+ RETURN
+
+ END
+
+***********************************************************************
+
+*** resid - Computes residual contributions to pressure (q), the
+* Helmloltz function (ar) , dP/dD (q5), the Gibbs function
+* (gr), entropy (sr), internal energy (ur), enthalpy (hr),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 708
More information about the CHNOSZ-commits
mailing list