subroutine haarH2O(t,pMPa,cp,cv,gH2O,hH2O,sH2O,rhof,vh2o,blkgas, * dvdth2o,sndspeed) C THIS PROGRAM CALCULATES THERMODYNAMIC PROPERTIES OF WATER AND C STEAM USING THE EQUATIONS OF HAAR ET AL. (1984). C PROGRAM AUTHORS:(1) MARK GHIORSO, DEPT. OF GEOLOGY & GEOPHYSICS C UNIVERSITY OF WASHINGTON C Seattle, WA 98195-1310 C ghiorso@u.washington.edu C (2) LARRY G. MASTIN, U.S. GEOLOGICAL SURVEY C CASCADES VOLCANO OBSERVATORY C lgmastin@usgs.gov C tel. 360-993-8925 (USA) C http://vulcan.wr.usgs.gov/Projects/Mastin C DATE: WRITTEN, 1998, REVISED DECEMBER 2002 C LANGUAGE: FORTRAN 90 (98% fortran 77) C THIS PROGRAM WAS ORIGINALLY WRITTEN BY MARK GHIORSO IN C. C IT WAS TRANSLATED INTO FORTRAN BY LARRY MASTIN, JUNE 1998, AND C DONATED TO THE UNIVERSITY OF NEW HAMPSHIRE CONDUIT MODELING C WORKSHOP WEB PAGE DECEMBER, 2002. C INPUT PARAMETERS: C t temperature, Kelvin C pMPa pressure, megapascals C OUTPUT PARAMETERS: C cp specific heat at constant pressure (J/kg K) C cv specific heat at constant volume (J/kg K) C gH2O Gibbs free energy (J/kg) C hH2O enthalpy relative to the elements (H and O) at STP (J/kg) C sH2O entropy relative to the elements at STP (J/kg K) C rhof density (kg/m3) C vh2o molar volume (cm3/mole) C blkgas bulk modulus (Pascals) C dvdth2o change in molar volume with temperature (m3/kg K) C sndspeed sound speed (m/s) implicit none double precision alpi(40), beti(40), bi(0:5), bbi(0:5) double precision r,gref,href,sref,t0,ps,rr,alpha,beta,gamma,p0 double precision taui(0:6), ermi(0:9), vol, rhn, dp, dr, rh, pr, * dpr, q10, qm, x, tr, dtrdt, dpdrh, dpdt, drhdt, d2pdt2, * d2pdtdrh, d2pdrh2, temp double precision ark, brk, oft, buk, cuk, duk, x1, x2, x2i, x3 integer i, icount C TERMS USED IN BASE FUNCTION double precision b, dbdt, d2bdt2, d3bdt3 double precision bb, dbbdt, d2bbdt2, d3bbdt3 double precision y, dydt, dydrh, d2ydt2, d2ydtdrh, d2ydrh2, * d3ydt3, d3ydt2drh, d3ydtdrh2, d3ydrh3 double precision gi(40),ci(18),rhoi(40),ttti(40) integer ki(40), li(40) C BASE FUNCTION AND ITS DERIVATIVES double precision Z, dZdt, dZdrh, d2Zdt2, d2Zdtdrh, d2Zdrh2, * d3Zdt3, d3Zdt2drh,d3Zdtdrh2, d3Zdrh3 double precision Ab, dAbdt, dAbdrh, d2Abdt2, d2Abdtdrh, * d2Abdrh2, d3Abdt3, d3Abdt2drh, d3Abdtdrh2, d3Abdrh3 C TERMS USED IN RESIDUAL FUNCTION double precision del, tau, abc C RESIDUAL FUNCTION AND ITS DERIVATIVES double precision Ar, dArdt, dArdrh, d2Ardt2, d2Ardtdrh, * d2Ardrh2, d3Ardt3, d3Ardt2drh, d3Ardtdrh2, d3Ardrh3 C IDEAL GAS FUNCTION AND ITS DERIVATIVES double precision Ai, dAidt, d2Aidt2, d3Aidt3 C HELMHOLTZ FUNCTION AND ITS DERIVATIVES double precision A, dAdt, dAdrh, d2Adt2, d2Adtdrh, d2Adrh2, * d3Adt3, d3Adt2drh, d3Adtdrh2, d3Adrh3 C FINAL THERMODYNAMIC VALUES double precision blkgas,cp,cpH2O,cv,cvH2O,dcpdtH2O, * gH2O,hH2O,vH2O,dvdtH2O,dvdpH2O,d2vdt2H2O, * d2vdtdpH2O,d2vdp2H2O,rhof,sH2O C OTHER TERMS double precision eps1,expQ,Q,dQdt,dQdrh,d2Qdtdrh, * d2Qdrh2,d2Qdt2,d3Qdt3, * d3Qdt2drh,d3Qdtdrh2,d3Qdrh3,p,pMPa,psat2, * sndspeed,t double precision square, cubic, quartic, quintic data gi/ * -.53062968529023e4, .22744901424408e5, .78779333020687e4, * -.69830527374994e3, .17863832875422e6, -.39514731563338e6, * .33803884280753e6, -.13855050202703e6, -.25637436613260e7, * .48212575981415e7, -.34183016969660e7, .12223156417448e7, * .11797433655832e8, -.21734810110373e8, .10829952168620e8, * -.25441998064049e7, -.31377774947767e8, .52911910757704e8, * -.13802577177877e8, -.25109914369001e7, .46561826115608e8, * -.72752773275387e8, .41774246148294e7, .14016358244614e8, * -.31555231392127e8, .47929666384584e8, .40912664781209e7, * -.13626369388386e8, .69625220862664e7, -.10834900096447e8, * -.22722827401688e7, .38365486000660e7, .68833257944332e5, * .21757245522644e6, -.26627944829770e5, -.70730418082074e6, * -.225e1, -1.68e1, .055e1, * -93.0e1/ data ki/ * 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, * 6, 6, 6, 6, 7, 7, 7, 7, 9, 9, 9, 9, 3, 3, 1, 5, 2, 2, 2, 4/ data li/ * 1, 2, 4, 6, 1, 2, 4, 6, 1, 2, 4, 6, 1, 2, 4, 6, 1, 2, 4, 6, * 1, 2, 4, 6, 1, 2, 4, 6, 1, 2, 4, 6, 0, 3, 3, 3, 0, 2, 0, 0/ data ci/ * .19730271018e2, .209662681977e2, -.483429455355, * .605743189245e1, 22.56023885, -9.87532442, *-.43135538513e1, .458155781, -.47754901883e-1, * .41238460633e-2, -.27929052852e-3, .14481695261e-4, *-.56473658748e-6, .16200446e-7, -.3303822796e-9, * .451916067368e-11, -.370734122708e-13, .137546068238e-15/ data rhoi/ * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.319, 0.310, 0.310, 1.55/ data ttti/ * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 640.0, 640.0, 641.6, 270.0/ data alpi/ * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 34.0, 40.0, 30.0, 1050.0/ data beti/ * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0e4, 2.0e4, 4.0e4, 25.0/ data bi/0.7478629, -0.3540782, 0.0, 0.007159876, * 0.0, -0.003528426/ data bbi/1.1278334, -0.5944001, -5.010996, 0.0, * 0.63684256, 0.0/ C COMMON BLOCKS common/xs/x1,x2,x2i,x3 C SET REFERENCE VALUES: C R SPECIFIC GAS CONSTANT C GREF GIBBS FREE ENERGY C HREF ENTHALPY C SREF ENTROPY C TO TEMPERATURE AT CRITICAL POINT C P PRESSURE AT CRITICAL POINT C RR UNIVERSAL GAS CONSTANT (J/MOLE K) C ALPHA CONSTANT USED TO CALCULATE BASE FUNCTION C BETA CONSTANT USED TO CALCULATE BASE FUNCTION C GAMMA CONSTANT USED TO CALCULATE BASE FUNCTION C PO PRESSURE AT 1 ATM (IN BARS) r = 4.6152 C -54955.2356146121147 25 c and 1 bar gref = -54955.23970014 href = -34099.89230644 sref = 69.94917790942 t0 = 647.25 ps = 220.55 rr = 8.31441 alpha = 11.0 beta = 133.0/3.0 gamma = 7.0/2.0 P0 = 1.01325 p = pMPa*10. !convert from MPa to bars C SET PRESSURE, CYCLE THROUGH INCREASING TEMPERATURES c write (6,12) c12 format (5x,'enter pressure in bars:'$) c read (5,*) p c write (6,13) c13 format (5x,'enter temperature, Celsius:'$) c read (5,*) t c t = t + 273.15 C The values of reduced temperature (T/T0)**i are stored in the C array TAUI(i) taui(0) = 1.0 taui(1) = t/t0 do 1015 i=2,6 1015 taui(i) = taui(i-1)*taui(1) b = bi(1)*log(taui(1)) + bi(0) + bi(3)/taui(3) + bi(5)/taui(5) bb = bbi(0) + bbi(1)/taui(1) + bbi(2)/taui(2) + bbi(4)/taui(4) if (t.le.647.25) ps = psat2(t) C************************************************************************** C set initial guess for rho using thb-fit to redlich-kwong C ark = redlich-kwong constant a C brk = redlich-kwong constant b C the redlich-kwong equation of state is: C p = r*t/(v-brk) - ark/(v*(v+brk)*sqrt(t)) C where t is absolute temp, and v is molar volume (cm3/mole) C******************************************************************** ark = 1.279186e8 - 2.241415e4 * t brk = 1.428062e1 + 6.092237e-4 * t oft = ark/(p*dsqrt(t)) buk = -10.0*rr*t/p cuk = oft - brk*brk + brk*buk duk = - brk*oft call kubik(buk,cuk,duk) C CALCULATE MOLAR VOLUME if (x2i.ne.0.0) then vol = x1 else if (p.lt.ps) then vol = MAX(x1, MAX(x2, x3)) else vol = MIN(x1, MIN(x2,x3)) end if end if C CALCULATE SPECIFIC VOLUME if (vol.le.0.0) then rhn = 1.9 else rhn = (1.0/vol)*18.0152 end if C************************************************************************** C CALCULATE PARAMETERS THAT GO INTO BASE FUNCTION: dp = 9.99e+10 dr = 9.99e+10 icount=1 eps1 = 10.0*epsilon(real(8)) c do while ((icount.le.100).and. c * ((dp.gt.eps1).or. c * (dr.gt.eps1))) do while (icount.le.100) if (dp.lt.1.0e-06) exit if (dr.lt.1.0e-06) exit rh = rhn if (rh.le.0.0) rh = 1.e-8 if (rh.gt.1.9) rh = 1.9 C CALCULATE Y, WHICH IS USED IN BASE FUNCTION (HAAR ET AL., P. 272) y = rh*b/4.0 C CALCULATE TERMS USED IN RESIDUAL FUNCTION (EQ. A.5 OF HAAR ET AL.) ermi(0) = 1.0 ermi(1) = 1.0-exp(-rh) do 1260 i=2,9 1260 ermi(i) = ermi(i-1)*ermi(1) pr = 0.0 dpr = 0.0 C CALCULATE RESIDUAL FUNCTION do 1270 i=1,36 pr = pr + gi(i)/taui(li(i))*ermi(ki(i)-1) dpr = dpr + (2.0+rh*(ki(i)*exp(-rh)-1.0)/ermi(1)) * *gi(i)/taui(li(i))*ermi(ki(i)-1) 1270 continue do 1280 i=37,40 del = rh/rhoi(i) - 1.0 tau = t/ttti(i) - 1.0 abc = -alpi(i)*del**ki(i)-beti(i)*tau*tau if (abc.gt.-100.00) then q10 = gi(i)*del**li(i) * exp(abc) else q10 = 0.0 end if qm = li(i)/del - ki(i)*alpi(i)*del**(ki(i)-1) pr = pr + q10*qm*rh*rh/rhoi(i) dpr = dpr + (q10*qm*rh*rh/rhoi(i)) * (2.0/rh+qm/rhoi(i)) * - rh*rh/(rhoi(i)*rhoi(i))*q10* * (li(i)/del/del + ki(i)*(ki(i)-1)*alpi(i) * *del**(ki(i)-2)) 1280 continue pr = rh*(rh*exp(-rh)*pr * + r*t*((1.0 + alpha*y + beta*y*y)/CUBIC(1.0-y) * + 4.0*y*(bb/b - gamma))) dpr = rh*exp(-rh)*dpr * + r*t*( (1.0 + 2.0*alpha*y + 3.0*beta*y*y)/CUBIC(1.0-y) * + 3.0*y*(1.0 + alpha*y + beta*y*y)/QUARTIC(1.0-y) * + 2.0*4.0*y*(bb/b - gamma)) if (dpr.le.0.0) then if (p.le.ps) then rhn = rhn*0.95 else rhn = rhn*1.05 end if else if (dpr.lt.0.01) dpr = 0.01 x = (p - pr)/dpr if (ABS(x).gt.0.1) x = 0.1*x/ABS(x) rhn = rh + x end if dp = ABS(1.0 - pr/p) dr = ABS(1.0 - rhn/rh) icount=icount+1 end do c write (6,1) t, p, ps, rhn c1 format (5x,'t=',f8.2,' p=',f8.2,' ps=',f8.2,' rhn=',f8.4) C CALCULATE DERIVATIVES W/R TO TEMPERATURE rh = rhn dbdt = bi(1)/t - 3.0*bi(3)/(taui(3)*t) - 5.0*bi(5)/(taui(5)*t) d2bdt2 = -bi(1)/(t*t) + 4.0*3.0*bi(3)/(taui(3)*t*t) * + 6.0*5.0*bi(5)/(taui(5)*t*t) d3bdt3 = 2.0*bi(1)/(t*t*t) - 5.0*4.0*3.0*bi(3)/(taui(3)*t*t*t) * - 7.0*6.0*5.0*bi(5)/(taui(5)*t*t*t) dbbdt = - bbi(1)/(taui(1)*t) - 2.0*bbi(2)/(taui(2)*t) * - 4.0*bbi(4)/(taui(4)*t) d2bbdt2 = 2.0*bbi(1)/(taui(1)*t*t) + 3.0*2.0*bbi(2)/(taui(2)*t*t) * + 5.0*4.0*bbi(4)/(taui(4)*t*t) d3bbdt3 = - 3.0*2.0*bbi(1)/(taui(1)*t*t*t) * - 4.0*3.0*2.0*bbi(2)/(taui(2)*t*t*t) * - 6.0*5.0*4.0*bbi(4)/(taui(4)*t*t*t) y = rh*b/4.0 dydt = rh*dbdt/4.0 dydrh = b/4.0 d2ydt2 = rh*d2bdt2/4.0 d2ydtdrh = dbdt/4.0 d2ydrh2 = 0.0 d3ydt3 = rh*d3bdt3/4.0 d3ydt2drh = d2bdt2/4.0 d3ydtdrh2 = 0.0 d3ydrh3 = 0.0 ermi(0) = 1.0 ermi(1) = 1.0-exp(-rh) do 1300 i=2,9 1300 ermi(i) = ermi(i-1)*ermi(1) C CALCULATE BASE FUNCTION AND ITS DERIVATIVES Z = - log(1.0-y) - (beta-1.0)/(1.0-y) * + (alpha+beta+1.0)/(2.0*(1.0-y)*(1.0-y)) + 4*y*(bb/b-gamma) * - (alpha-beta+3.0)/2.0 + log(rh*r*t/P0) dZdt = 1.0/t + 4.0*(b*dbbdt-bb*dbdt)*y/(b*b) * + 4.0*(bb/b-gamma)*dydt * + dydt/(1.0-y) + (alpha+beta+1.0)*dydt/CUBIC(1.0-y) * - (beta-1.0)*dydt/SQUARE(1.0-y) dZdrh = 1.0/rh + 4.0*(bb/b-gamma)*dydrh + dydrh/(1.0-y) * + (alpha+beta+1.0)*dydrh/CUBIC(1.0-y) * - (beta-1.0)*dydrh/SQUARE(1.0-y) d2Zdt2 = -1.0/(t*t) + 4.0*(((b*d2bbdt2-bb*d2bdt2)*(b*b) * -2.0*(b*dbbdt-bb*dbdt)*b*dbdt)*y)/QUARTIC(b) * + SQUARE(dydt/(1.0-y)) * + 3.0*(alpha+beta+1.0)*SQUARE(dydt)/QUARTIC(1.0-y) * - 2.0*(beta-1.0)*SQUARE(dydt)/CUBIC(1.0-y) * + 8.0*(b*dbbdt-bb*dbdt)*dydt/(b*b) * + 4.0*(bb/b-gamma)*d2ydt2 * + d2ydt2/(1.0-y) + (alpha+beta+1.0)*d2ydt2/CUBIC(1.0-y) * - (beta-1.0)*d2ydt2/SQUARE(1.0-y); d2Zdtdrh = 4.0*(bb/b-gamma)*d2ydtdrh + d2ydtdrh/(1.0-y) * + (alpha+beta+1.0)*d2ydtdrh/CUBIC(1.0-y) * - (beta-1.0)*d2ydtdrh/SQUARE(1.0-y) * + 4.0*(b*dbbdt-bb*dbdt)*dydrh/(b*b) * + dydt*dydrh/SQUARE(1.0-y) * + 3.0*(alpha+beta+1.0)*dydt*dydrh/QUARTIC(1.0-y) * - 2.0*(beta-1.0)*dydt*dydrh/CUBIC(1.0-y) d2Zdrh2 = -1.0/(rh*rh) + SQUARE(dydrh/(1.0-y)) * + 3.0*(alpha+beta+1.0)*SQUARE(dydrh)/QUARTIC(1.0-y) * - 2.0*(beta-1.0)*SQUARE(dydrh)/CUBIC(1.0-y) * + 4.0*(bb/b-gamma)*d2ydrh2 + d2ydrh2/(1.0-y) * + (alpha+beta+1.0)*d2ydrh2/CUBIC(1.0-y) * - (beta-1.0)*d2ydrh2/SQUARE(1.0-y) d3Zdt3 = 2.0/CUBIC(t) + 4.0*(((-2.0*(b*dbbdt-bb*dbdt)*SQUARE(dbdt) * + (b*d3bbdt3+d2bbdt2*dbdt-dbbdt*d2bdt2-bb*d3bdt3)*SQUARE(b) * - 2.0*(b*dbbdt-bb*dbdt)*b*d2bdt2)*y * + ((b*d2bbdt2-bb*d2bdt2)*SQUARE(b) * - 2.0*(b*dbbdt-bb*dbdt)*b*dbdt)*dydt)*QUARTIC(b) * - 4.0*((b*d2bbdt2-bb*d2bdt2)*SQUARE(b) * -2.0*(b*dbbdt-bb*dbdt)*b*dbdt)*CUBIC(b)*y*dbdt) * /(QUARTIC(b)*QUARTIC(b)) * + 2.0*CUBIC(dydt/(1.0-y)) * + 12.0*(alpha+beta+1.0)*CUBIC(dydt)/QUINTIC(1.0-y) * - 6.0*(beta-1.0)*CUBIC(dydt)/QUARTIC(1.0-y) * + 8.0*((b*d2bbdt2-bb*d2bdt2)*SQUARE(b) * -2.0*(b*dbbdt-bb*dbdt)*b*dbdt)*dydt/QUARTIC(b) * + 12.0*(b*dbbdt-bb*dbdt)*d2ydt2/SQUARE(b) * + 3.0*dydt*d2ydt2/SQUARE(1.0-y) * + 9.0*(alpha+beta+1.0)*dydt*d2ydt2/QUARTIC(1.0-y) * - 6.0*(beta-1.0)*dydt*d2ydt2/CUBIC(1.0-y) * + 4.0*(bb/b-gamma)*d3ydt3 + d3ydt3/(1.0-y) * + (alpha+beta+1.0)*d3ydt3/CUBIC(1.0-y) * - (beta-1.0)*d3ydt3/SQUARE(1.0-y); d3Zdt2drh = 4.0*(((b*d2bbdt2-bb*d2bdt2)*SQUARE(b) * -2.0*(b*dbbdt-bb*dbdt)*b*dbdt)*dydrh)/QUARTIC(b) * + 4.0*(bb/b-gamma)*d3ydt2drh + d3ydt2drh/(1.0-y) * + (alpha+beta+1.0)*d3ydt2drh/CUBIC(1.0-y) * - (beta-1.0)*d3ydt2drh/SQUARE(1.0-y) * + 8.0*(b*dbbdt-bb*dbdt)*d2ydtdrh/SQUARE(b) * + 2.0*dydt*d2ydtdrh/SQUARE(1.0-y) * + 6.0*(alpha+beta+1.0)*dydt*d2ydtdrh/QUARTIC(1.0-y) * - 4.0*(beta-1.0)*dydt*d2ydtdrh/CUBIC(1.0-y) * + 2.0*SQUARE(dydt)*dydrh/CUBIC(1.0-y) * + 12.0*(alpha+beta+1.0)*SQUARE(dydt)*dydrh/QUINTIC(1.0-y) * - 6.0*(beta-1.0)*SQUARE(dydt)*dydrh/QUARTIC(1.0-y) * + d2ydt2*dydrh/SQUARE(1.0-y) * + 3.0*(alpha+beta+1.0)*d2ydt2*dydrh/QUARTIC(1.0-y) * - 2.0*(beta-1.0)*d2ydt2*dydrh/CUBIC(1.0-y); d3Zdtdrh2 = 2.0*dydt*SQUARE(dydrh)/CUBIC(1.0-y) * + 12.0*(alpha+beta+1.0)*dydt*SQUARE(dydrh)/QUINTIC(1.0-y) * - 6.0*(beta-1.0)*dydt*SQUARE(dydrh)/QUARTIC(1.0-y) * + 4.0*(bb/b-gamma)*d3ydtdrh2 + d3ydtdrh2/(1.0-y) * + (alpha+beta+1.0)*d3ydtdrh2/CUBIC(1.0-y) * - (beta-1.0)*d3ydtdrh2/SQUARE(1.0-y) * + 2.0*d2ydtdrh*dydrh/SQUARE(1.0-y) * + 6.0*(alpha+beta+1.0)*d2ydtdrh*dydrh/QUARTIC(1.0-y) * - 4.0*(beta-1.0)*d2ydtdrh*dydrh/CUBIC(1.0-y) * + 4.0*(b*dbbdt-bb*dbdt)*d2ydrh2/SQUARE(b) * + dydt*d2ydrh2/SQUARE(1.0-y) * + 3.0*(alpha+beta+1.0)*dydt*d2ydrh2/QUARTIC(1.0-y) * - 2.0*(beta-1.0)*dydt*d2ydrh2/CUBIC(1.0-y) d3Zdrh3 = 2.0/CUBIC(rh) + 2.0*CUBIC(dydrh/(1.0-y)) * + 12.0*(alpha+beta+1.0)*CUBIC(dydrh)/QUINTIC(1.0-y) * - 6.0*(beta-1.0)*CUBIC(dydrh)/QUARTIC(1.0-y) * + 3.0*dydrh*d2ydrh2/SQUARE(1.0-y) * + 9.0*(alpha+beta+1.0)*dydrh*d2ydrh2/QUARTIC(1.0-y) * - 6.0*(beta-1.0)*dydrh*d2ydrh2/CUBIC(1.0-y) * + 4.0*(bb/b-gamma)*d3ydrh3 + d3ydrh3/(1.0-y) * + (alpha+beta+1.0)*d3ydrh3/CUBIC(1.0-y) * - (beta-1.0)*d3ydrh3/SQUARE(1.0-y) Ab = r*t*Z dAbdt = r*Z + r*t*dZdt dAbdrh = r*t*dZdrh d2Abdt2 = 2.0*r*dZdt + r*t*d2Zdt2 d2Abdtdrh = r*dZdrh + r*t*d2Zdtdrh d2Abdrh2 = r*t*d2Zdrh2 d3Abdt3 = 3.0*r*d2Zdt2 + r*t*d3Zdt3 d3Abdt2drh = 2.0*r*d2Zdtdrh + r*t*d3Zdt2drh d3Abdtdrh2 = r*d2Zdrh2 + r*t*d3Zdtdrh2 d3Abdrh3 = r*t*d3Zdrh3 C CALCULATE RESIDUAL FUNCTION AND ITS DERIVATIVES Ar = 0.0 dArdt = 0.0 dArdrh = 0.0 d2Ardt2 = 0.0 d2Ardtdrh = 0.0 d2Ardrh2 = 0.0 d3Ardt3 = 0.0 d3Ardt2drh = 0.0 d3Ardtdrh2 = 0.0 d3Ardrh3 = 0.0 do 1400 i=1,36 Ar = Ar + gi(i)/ki(i)/taui(li(i))*ermi(ki(i)) dArdt = dArdt + (-li(i)*gi(i)/ki(i)/(taui(li(i))*t) * * ermi(ki(i))) dArdrh = dArdrh + gi(i)/taui(li(i))*ermi(ki(i)-1)*exp(-rh) d2Ardt2 = d2Ardt2 + (li(i)+1.0)*li(i)*gi(i)/ki(i) * /(taui(li(i))*t*t)*ermi(ki(i)) d2Ardtdrh = d2Ardtdrh + (-li(i)*gi(i)/ * (taui(li(i))*t)*ermi(ki(i)-1)*exp(-rh)) if (ki(i).gt.1) then d2Ardrh2 = d2Ardrh2 +gi(i)/taui(li(i)) * *((ki(i)-1.0)*ermi(ki(i)-2) * *exp(-rh)-ermi(ki(i)-1))*exp(-rh) else d2Ardrh2 = d2Ardrh2 + (-gi(i)/taui(li(i))*exp(-rh)) end if d3Ardt3 = d3Ardt3 - (li(i)+2.0)*(li(i)+1.0)*li(i)*gi(i)/ki(i) * /(taui(li(i))*t*t*t)*ermi(ki(i)) d3Ardt2drh = d3Ardt2drh +(li(i)+1.0)*li(i)*gi(i) * /(taui(li(i))*t*t)*ermi(ki(i)-1)*exp(-rh) if (ki(i).gt.1) then d3Ardtdrh2 = d3Ardtdrh2 - li(i)*gi(i)/(taui(li(i))*t) * *((ki(i)-1.0)*ermi(ki(i)-2) * *exp(-rh)-ermi(ki(i)-1))*exp(-rh) else d3Ardtdrh2 = d3Ardtdrh2 + * li(i)*gi(i)/(taui(li(i))*t)*exp(-rh) end if if (ki(i).gt.2) then d3Ardrh3 = d3Ardrh3 + gi(i)/taui(li(i)) * *(((ki(i)-2.0)*ermi(ki(i)-3)*exp(-rh) * - 3.0*ermi(ki(i)-2)) * *(ki(i)-1.0)*exp(-rh)+ermi(ki(i)-1))*exp(-rh) else if (ki(i).gt.1) then d3Ardrh3 = d3Ardrh3 - gi(i)/taui(li(i)) * *(4.0*exp(-rh)-1.0)*exp(-rh) else d3Ardrh3 = d3Ardrh3 + gi(i)/taui(li(i))*exp(-rh) end if 1400 continue do 1500 i=37,40 del = rh/rhoi(i) - 1.0 tau = t/ttti(i) - 1.0 Q = -alpi(i)*del**ki(i) - beti(i)*tau*tau dQdt = -beti(i)*2.0*tau/ttti(i) if (ki(i).eq.0) then dQdrh = 0.0 else dQdrh = -alpi(i)*ki(i)*del**(ki(i)-1)/rhoi(i) end if d2Qdt2 = -beti(i)*2.0/SQUARE(ttti(i)) d2Qdtdrh = 0.0 if ((ki(i).eq.0).or.(ki(i).eq.1)) then d2Qdrh2 = 0.0 else d2Qdrh2 = -alpi(i)*ki(i)*(ki(i)-1.0)*del**(ki(i)-2) * /SQUARE(rhoi(i)) end if d3Qdt3 = 0.0 d3Qdt2drh = 0.0 d3Qdtdrh2 = 0.0 if ((ki(i).eq.0).or.(ki(i).eq.1).or.(ki(i).eq.2)) then d3Qdrh3 = 0.0 else d3Qdrh3 = -alpi(i)*ki(i)*(ki(i)-1.0)*(ki(i)-2.0) * *del**(ki(i)-3)/CUBIC(rhoi(i)) end if if (Q.gt.-100.0) then expQ = dexp(Q) else expQ = 0.0 end if Ar = Ar + gi(i)*del**li(i)*expQ dArdt = dArdt + gi(i)*del**li(i)*expQ*dQdt if (li(i).eq.0) then dArdrh = dArdrh + gi(i)*expQ*dQdrh else dArdrh = dArdrh + gi(i)*li(i)*del**(li(i)-1)*expQ/rhoi(i) * + gi(i)*del**li(i)*expQ*dQdrh end if d2Ardt2 = d2Ardt2 + gi(i)*del**li(i)*expQ*(SQUARE(dQdt)+d2Qdt2) if (li(i).eq.0) then d2Ardtdrh = d2Ardtdrh + gi(i)*expQ*dQdt*dQdrh * + gi(i)*expQ*d2Qdtdrh else d2Ardtdrh = d2Ardtdrh + gi(i)*li(i) * *del**(li(i)-1)*expQ*dQdt/rhoi(i) * + gi(i)*del**li(i)*expQ*(dQdt*dQdrh+d2Qdtdrh) end if if (li(i).eq.0) then d2Ardrh2 = d2Ardrh2 + gi(i)*expQ*SQUARE(dQdrh) * + gi(i)*expQ*d2Qdrh2 else if (li(i).eq.1) then d2Ardrh2 = d2Ardrh2 + 2.0*gi(i)*expQ*dQdrh/rhoi(i) * + gi(i)*del*expQ*SQUARE(dQdrh) + gi(i)*del*expQ*d2Qdrh2 else d2Ardrh2 = d2Ardrh2 + gi(i)*li(i)*(li(i)-1.0)*del**(li(i)-2) * *expQ/SQUARE(rhoi(i)) * + 2.0*gi(i)*li(i)*(del**(li(i)-1))*expQ*dQdrh/rhoi(i) * + gi(i)*del**li(i)*(expQ*SQUARE(dQdrh)+expQ*d2Qdrh2) end if d3Ardt3 = d3Ardt3 + gi(i)*del**li(i) * *expQ*(CUBIC(dQdt)+3.0*dQdt*d2Qdt2+d3Qdt3) if (li(i).eq.0) then d3Ardt2drh = d3Ardt2drh + gi(i)*(expQ*SQUARE(dQdt)*dQdrh * + expQ*(d2Qdt2*dQdrh + dQdt*d2Qdtdrh)) * + gi(i)*(expQ*dQdt*d2Qdtdrh + expQ*d3Qdt2drh) else d3Ardt2drh = d3Ardt2drh + gi(i)*li(i)*del**(li(i)-1) * *(expQ*SQUARE(dQdt)+expQ*d2Qdt2)/rhoi(i) * + gi(i)*del**li(i)*(expQ*dQdt*(dQdt*dQdrh+d2Qdtdrh) * + expQ*(d2Qdt2*dQdrh+dQdt*d2Qdtdrh+d3Qdt2drh)) end if if (li(i).eq.1) then d3Ardtdrh2 = d3Ardtdrh2 + 2.0*gi(i)*li(i)*(expQ*dQdt*dQdrh * + expQ*d2Qdtdrh)/rhoi(i) * + gi(i)*del*(expQ*dQdt*SQUARE(dQdrh) * + expQ*2.0*dQdrh*d2Qdtdrh * + expQ*dQdt*d2Qdrh2 + expQ*d3Qdtdrh2) else d3Ardtdrh2 = d3Ardtdrh2 + gi(i)*li(i)*(li(i)-1.0) * *del**(li(i)-2)*expQ*dQdt/SQUARE(rhoi(i)) * + 2.0*gi(i)*li(i)*del**(li(i)-1) * *(expQ*dQdt*dQdrh + expQ*d2Qdtdrh)/rhoi(i) * + gi(i)*del**li(i)*(expQ*dQdt*SQUARE(dQdrh) * + expQ*2.0*dQdrh*d2Qdtdrh + expQ*dQdt*d2Qdrh2 * + expQ*d3Qdtdrh2) end if if (li(i).eq.0) then d3Ardrh3 = d3Ardrh3 + gi(i)*(expQ*CUBIC(dQdrh) * + expQ*2.0*dQdrh*d2Qdrh2 * + expQ*dQdrh*d2Qdrh2 + expQ*d3Qdrh3) else if (li(i).eq.1) then d3Ardrh3 = d3Ardrh3 + 3.0*gi(i)*(expQ*SQUARE(dQdrh) * + expQ*d2Qdrh2)/rhoi(i) * + gi(i)*del*(expQ*CUBIC(dQdrh) + expQ*2.0*dQdrh*d2Qdrh2 * + expQ*dQdrh*d2Qdrh2 + expQ*d3Qdrh3) else if (li(i).eq.2) then d3Ardrh3 = d3Ardrh3 + 3.0*gi(i)*2.0*expQ*dQdrh/SQUARE(rhoi(i)) * + 3.0*gi(i)*2.0*del*(expQ*SQUARE(dQdrh) * + expQ*d2Qdrh2)/rhoi(i) * + gi(i)*SQUARE(del)*(expQ*CUBIC(dQdrh) * + expQ*2.0*dQdrh*d2Qdrh2 * + expQ*dQdrh*d2Qdrh2 + expQ*d3Qdrh3) else d3Ardrh3 = d3Ardrh3 + gi(i)*li(i)*(li(i)-1.0)*(li(i)-2.0) * *del**(li(i)-3)*expQ/CUBIC(rhoi(i)) * + 3.0*gi(i)*li(i)*(li(i)-1.0)*del**(li(i)-2) * *expQ*dQdrh/SQUARE(rhoi(i)) * + 3.0*gi(i)*li(i)*del**(li(i)-1) * *(expQ*SQUARE(dQdrh) + expQ*d2Qdrh2)/rhoi(i) * + gi(i)*del**li(i)*(expQ*CUBIC(dQdrh) * + expQ*2.0*dQdrh*d2Qdrh2 * + expQ*dQdrh*d2Qdrh2 + expQ*d3Qdrh3) end if 1500 continue C CALCULATE IDEAL GAS FUNCTION tr = t/1.0e2 dtrdt = 1.0/100.0 Z = 1.0 + (ci(1)/tr + ci(2))*log(tr) dZdt = (-ci(1)*dtrdt/SQUARE(tr))*log(tr) * + (ci(1)/tr + ci(2))*dtrdt/tr d2Zdt2 = (2.0*ci(1)*SQUARE(dtrdt)/CUBIC(tr))*log(tr) * + (-ci(1)*dtrdt/SQUARE(tr))*dtrdt/tr * + (-ci(1)*dtrdt/SQUARE(tr))*dtrdt/tr * - (ci(1)/tr + ci(2))*SQUARE(dtrdt/tr) d3Zdt3 = (-3.0*2.0*ci(1)*CUBIC(dtrdt)/QUARTIC(tr))*log(tr) * + (2.0*ci(1)*SQUARE(dtrdt)/CUBIC(tr))*dtrdt/tr * + (2.0*ci(1)*SQUARE(dtrdt)/CUBIC(tr))*dtrdt/tr * - (-ci(1)*dtrdt/SQUARE(tr))*SQUARE(dtrdt/tr) * + (2.0*ci(1)*SQUARE(dtrdt)/CUBIC(tr))*dtrdt/tr * - (-ci(1)*dtrdt/SQUARE(tr))*SQUARE(dtrdt/tr) * - (-ci(1)*dtrdt/SQUARE(tr))*SQUARE(dtrdt/tr) * + 2.0*(ci(1)/tr + ci(2))*CUBIC(dtrdt/tr) do 1600 i=3,18 Z = Z + ci(i)*tr**(i-6) dZdt = dZdt + (i-6.0)*ci(i)*tr**(i-7)*dtrdt d2Zdt2 = d2Zdt2 + (i-6.0)*(i-7.0)*ci(i)*tr**(i-8)*SQUARE(dtrdt) d3Zdt3 = d3Zdt3 + (i-6.0)*(i-7.0)*(i-8.0)*ci(i)*tr**(i-9) * *CUBIC(dtrdt) 1600 continue Ai = -r*t*Z dAidt = - r*Z - r*t*dZdt d2Aidt2 = - 2.0*r*dZdt - r*t*d2Zdt2 d3Aidt3 = - 3.0*r*d2Zdt2 - r*t*d3Zdt3 C CALCULATE THE TOTAL HELMHOLTZ FUNCTION AND ITS DERIVIATIVES A = Ab + Ar + Ai dAdt = dAbdt + dArdt + dAidt dAdrh = dAbdrh + dArdrh d2Adt2 = d2Abdt2 + d2Ardt2 + d2Aidt2 d2Adtdrh = d2Abdtdrh + d2Ardtdrh d2Adrh2 = d2Abdrh2 + d2Ardrh2 d3Adt3 = d3Abdt3 + d3Ardt3 + d3Aidt3 d3Adt2drh = d3Abdt2drh + d3Ardt2drh d3Adtdrh2 = d3Abdtdrh2 + d3Ardtdrh2 d3Adrh3 = d3Abdrh3 + d3Ardrh3 C CALCULATE G = A + P/RH AND V = 1/RH c p = rh*rh*dAdrh dpdrh = 2.0*rh*dAdrh + rh*rh*d2Adrh2 dpdt = rh*rh*d2Adtdrh drhdt = -dpdt/dpdrh d2pdt2 = rh*rh*d3Adt2drh d2pdtdrh = 2.0*rh*d2Adtdrh + rh*rh*d3Adtdrh2 d2pdrh2 = 2.0*dAdrh + 4.0*rh*d2Adrh2 + rh*rh*d3Adrh3 C CALCULATE OTHER THERMODYNAMIC PROPERTIES C ("TEMP" IS A TEMPORARY VARIABLE, NOT TEMPERATURE) gH2O = A + p/rh hH2O = A + p/rh - t*dAdt sH2O = - dAdt cpH2O = -t*d2Adt2 + (t/(rh*rh))*SQUARE(dpdt)/(dpdrh) cvH2O = -t * d2Adt2 vH2O = 1.0/rh dvdtH2O = -(1.0/SQUARE(rh))*drhdt dvdpH2O = -(1.0/SQUARE(rh))/dpdrh temp = (-d2pdt2 - 2.0*d2pdtdrh*drhdt * - d2pdrh2*SQUARE(drhdt))/dpdrh d2vdt2H2O = 2.0*SQUARE(drhdt)/CUBIC(rh) - temp/SQUARE(rh) temp = (-d2pdtdrh/dpdrh - d2pdrh2*drhdt/dpdrh)/dpdrh d2vdtdpH2O = -2.0*dpdt/(CUBIC(rh)*SQUARE(dpdrh)) - temp/SQUARE(rh) d2vdp2H2O = (1.0/SQUARE(rh))*d2pdrh2/CUBIC(dpdrh) * + (2.0/CUBIC(rh))/SQUARE(dpdrh) temp = - d2Adt2 - t*d3Adt3 + SQUARE(dpdt/rh)/dpdrh * + t*(2.0*dpdrh*dpdt*d2pdt2 - SQUARE(dpdt)*d2pdtdrh) * /SQUARE(rh*dpdrh); dcpdtH2O = temp + t*(d2vdt2H2O)*dpdt C CONVERT FROM SPECIFIC (PER G) VALUES TO MOLAR VALUES gH2O = gH2O * 1.80152 hH2O = hH2O * 1.80152 sH2O = sH2O * 1.80152 cpH2O = cpH2O * 1.80152 cvH2O = cvH2O * 1.80152 dcpdtH2O = dcpdtH2O * 1.80152 vH2O = vH2O * 18.0152 dvdtH2O = dvdtH2O * 18.0152 dvdpH2O = dvdpH2O * 18.0152 d2vdt2H2O = d2vdt2H2O * 18.0152 d2vdtdpH2O = d2vdtdpH2O * 18.0152 d2vdp2H2O = d2vdp2H2O * 18.0152 C SUBSTRACT OUT BERMAN'S VALUE (CHECK THIS!) C Berman 1988 Haar 1977 gH2O = gH2O - 285829.96 - (298.15*69.9146) - gref hH2O = hH2O - 285829.96 - href C CONVERT TO SI UNITS cp = cpH2O/.0180152 cv = cvH2O/.0180152 rhof = (1.0D+06*.0180152)/vH2O gamma = cpH2O/cvH2O blkgas = rhof*dpdrh*100. sndspeed = sqrt((cp/cv)*blkgas/rhof) C CONVERT DVDTH2O FROM CM3/MOLE TO M3/(KG K) dvdth2o = dvdth2o/(0.0180152*1.0d+06) return end C*************************************************************************** C FUNCTION PSAT C FUNCTION THAT CALCULATES THE SATURATION PRESSURE AT A GIVEN TEMPERATURE C**************************************************************************** function psat2(t) implicit double precision (a-h,o-z) double precision a(8), ff, t, v, w, wsq integer i data a/-7.8889166, 2.5514255, -6.716169, 33.239495, * -105.38479, 174.35319, -148.39348, 48.631602/ if (t.le.314.0) then psat2 = exp(6.3573118-8858.843/t + 607.56335/(t**0.6)) return else v = t/647.25 w = dabs(1.0-v*.9997) wsq = sqrt(w) ff = 0.0 do 2340 i=1,8 z = i ff = ff + a(i)*w**((z+1.)/2.) 2340 continue psat2 = 1.001*220.93*exp(ff/v) end if return end C************************************************************************* C SUBROUTINE KUBIK C************************************************************************* subroutine kubik(b,c,d) implicit double precision (a-h,o-z) double precision b,c,d,ff,p,phi3,pi,q,r,x1,x2,x2i,x3 common/xs/x1,x2,x2i,x3 C SET CONSTANTS, INITIALIZE VARIABLES pi = 3.14159263538979 x2 = 0.0 x2i = 0.0 x3 = 0.0 if ((c.eq.0.0).and.(d.eq.0.0)) then x1 = -b return end if q = (2.0*cubic(b)/27.0 - b*c/3.0 + d)/2.0 p = (3.0*c - b*b)/9.0 ff = dabs(p) r = sqrt(ff) ff = r*q if (ff.lt.0.0) r = -r ff = q/cubic(r); if (p.gt.0.0) then phi3 = log(ff + sqrt(ff*ff+1.0))/3.0 x1 = -r*(exp(phi3) - exp(-phi3)) - b/3.0 x2i = 1.0 else if (q*q + p*p*p > 0.0) then phi3 = log(ff + sqrt(ff*ff-1.0))/3.0 x1 = -r*(exp(phi3) + exp(-phi3)) - b/3.0 x2i = 1.0 else phi3 = atan(sqrt(1.0-ff*ff)/ff)/3.0 x1 = -2.0*r*cos(phi3) - b/3.0 x2 = 2.0*r*cos(pi/3.0-phi3) - b/3.0 x2i = 0.0 x3 = 2.0*r*cos(pi/3.0+phi3) - b/3.0 end if end if return end C**************************************************************************** C SQUARE, CUBIC, QUARTIC, QUINTIC FUNCTIONS C*************************************************************************** function SQUARE(x) double precision x, SQUARE square = x*x return end function CUBIC(x) double precision x, CUBIC cubic = x*x*x return end function QUARTIC(x) double precision x, QUARTIC quartic = x*x*x*x return end function QUINTIC(x) double precision x, QUINTIC quintic = x*x*x*x*x return end