subroutine melt(tabs,pMPa,SiO2,Al2O3,Fe2O3,FeO,MgO,CaO,TiO2, * Na2O,K2O,dH2O,xcomp,blkmag,cpm,enthalpy,entropy, * mwttot,rhom,hwm,swm,gwm,dgdxw,dvmdt,sndspeed) C SUBROUTINE THAT CALCULATES THERMODYNAMIC PROPERTIES OF MELTS USING C THE METHODOLOGY OF GHIORSO & SACK (1995; CONTRIB. MIN. PET., 199: C 197-212). C PROGRAM AUTHOR: 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 SUBROUTINE CALLS THE ADDITIONAL SUBROUTINE: C haarH2O C INPUT VARIABLES C tabs temperature (Kelvin) C pMPa pressure (MPa) C SiO2 wt% SiO2 in melt (anhydrous) C Al2O3 wt% Al2O3 in melt (anhydrous) C Fe2O3 wt% Fe2O3 in melt (anhydrous) C FeO wt% FeO in melt (anhydrous) C MgO wt% MgO in melt (anhydrous) C CaO wt% CaO in melt (anhydrous) C TiO2 wt% TiO2 in melt (anhydrous) C Na2O wt% Na2O in melt (anhydrous) C K2O wt% K2O in melt (anhydrous) C dH2O wt% dissolved water in melt C OUTPUT VARIABLES C blkmag bulk modulus of melt (Pascals) C cpm specific heat of melt (J/kg K) C enthalpy enthalpy of melt, relative to that of constituent elements at STP (J/kg) C entropy entropy of melt relative to that of constituent elements at STP (J/kg K) C mwttot total molecular weight of melt (kg/mole), calculating moles as described below C rhom melt density (kg/m3) C xcomp(10) mole fractions of the various end-member components C hwm partial molar enthalpy of dissolved water (J/kg) C swm partial molar entropy of dissolved water (J/kg K) C dgdxw partial of Gibbbs free energy with mole fraction water (J/mole) C dvmdt change in specific volume with temperature (m3/kg K) C sndspeed speed of sound in melt (m/s) C THE METHOD OF GHIORSO AND SACK CALCULATES THERMODYNAMIC PROPERTIES OF THE MELT BY CONSIDERING C IT TO BE A SOLUTION OF END-MEMBER COMPONENTS WHOSE PROPERTIES ARE WELL-CHARACTHERIZED. C THOSE END-MEMBER COMPONENTS ARE THE FOLLOWING MINERALS: C COMPOUND NAME INDEX NUMBER C SIO2 SILICA (AMORPHOUS) 1 C TIO2 RUTILE 2 C AL2O3 CORUNDUM 3 C FE2O3 HEMATITE 4 C FE2SIO4 FAYALITE 5 C MG2SIO4 FORSTERITE 6 C CASIO3 PSEUDOWOLLASTONITE 7 C NA2SIO3 SODIUM METASILICATE 8 C KALSIO4 KALSILITE 9 C H2O WATER (LIQUID) 10 C C VARIABLE NAMES (ALL VALUES ARE PER MOLE) C COMP COMPRESSIBILITY OF WATER (1/MPA) C CPLIQ(10) CP OF LIQUID C DH2O MASS FRACTION WATER DISSOLVED IN THE TOTAL SYSTEM (MELT+XTLS+GAS) C DPDT CHANGE IN PRESSURE W/TEMP FOR WATER (MPA/K) C DVDT0(10) CHANGE IN SPECIFIC VOLUME WITH TEMP AT STANDARD T AND P C DVDT(10) CHANGE IN SPECIFIC VOLUME WITH TEMP AT SPECIFIED T AND P C DVDP0(10) CHANGE IN SPECIFIC VOLUME WITH PRESSURE AT STANDARD T AND P C DVDP(10) CHANGE IN VOLUME WITH PRESSURE AT SPECIFIED T AND P C D2VDTDP(10) SECOND DERIVATIVE OF SP. VOLUME WITH T AND P C D2VDP2(10) SECOND DERIVATIVE OF SP. VOLUME WITH P C GAPP(10) APPARENT GIBBS FREE ENERGY OF MELT (J/MOLE) C GFUSION(10) GIBBS FREE ENERGY OF THE LIQUID AT THE MELTING TEMP C HAPP(10) APPARENT ENTHALPY OF MELT (J/MOLE) C HFUSION(10) ENTHALPY OF LIQUID AT TFUSION (J/MOLE) C HWM PARTIAL MOLAR ENTHALPY OF DISSOLVED WATER C LABEL(10) LABELS OF END-MEMBER COMPONENTS C LABELOX(10) LABELS OF OXIDE COMPONENTS C MASSOX(10) MOLECULAR WEIGHTS OF OXIDE COMPONENTS C MCOMP(10) MOLES OF END-MEMBER COMPONENTS C MWTCOMP(10) MOLECULAR WEIGHT OF EACH COMPONENT (KG/MOLE) C MWTTOT AVERAGE MOLECULAR WEIGHT OF MELT (FROM END-MEMBER COMPONENTS) C PR STANDARD PRESSURE (1 ATM) C PBARS PRESSURE (BARS) C R UNIVERSAL GAS CONSTANT (8.314 J/MOLE K) C RHO1ATM DENSITY AT 1 ATM, T (KG/M3) C SAPP(10) APPARENT ENTROPY OF MELT C SFUSION(10) THIRD-LAW ENTROPY OF LIQUID AT THE MELTING TEMPERATURE C SUMWXX(10) SUM OF W(I,J)*X(J) C SWM PARTIAL MOLAR ENTROPY OF DISSOLVED WATER IN MELT C TFUSION(10) MELTING TEMPERATURE (K) C TR STANDARD TEMPERATURE (298.15 K) C VO(10) SPECIFIC VOLUME AT P=1 ATM, T=1673 K) C V(10) SPECIFIC VOLUME OF EACH COMPONENT AT SPECIFIED T, P (J/BAR) C VTOT SPECIFIC VOLUME OF TOTAL MELT C WTOXAN(9) WEGHT PERCENT OF OXIDE COMPONENTS (ANHYDROUS) C WTOXIDE(10) WEIGHT PERCENT OF OXIDE COMPONENTS (HYDROUS) C I USE UNITS OF TEMPERATURE (K) AND PRESSURE (BARS) TO REMAIN C CONSISTENT WITH OTHER THERMODYNAMIC MODELS OF MINERALS AND C MELTS (GHIORSO & SACK, 1995; BERMAN, 1988). implicit none real*8 cpliq(10), dvdp(10), dvdp0(9), dvdt(10), dvdt0(9), * d2vdpdt(9), d2vdp2(9), * gapp(10), gwm, happ(10), * hfusion(9), hwm, massox(10), * mcomp(10), mwtcomp(10), * sapp(10), sfusion(9), sumwxx(10), swm, * tfusion(9), vo(9), v(10), w(10,10), wtoxan(9), * wtoxide(10), xcomp(10), xoxide(10) real*8 SiO2,Al2O3,Fe2O3,FeO,MgO,CaO,TiO2,Na2O,K2O,totox real*8 ah2o,bh2o,blkmag,cpm,dgdt,d2gdt2,d3gdt3,dgdxw, dh2o, * dtmp,dvmdp,dvmdt,enthalpy,entropy,moles, mwttot, * pbars,pMPa,pr,r,rhom,sndspeed,sumgibbs, * sums,sumxcomp,sumh,sumw, * sumwh,sumwx,tabs,tr,v1atm,vdp,vtot real*8 cpg,cvg,gH2O,hH2O,sH2O,vH2O,dvdth2o,rhof,blkgas,ch2o real*8 grobie,dgrobiedt,d2grobiedt2,d3grobiedt3 real*8 phip,dphipdt,d2phipdt2,d3phipdt3,dphipdp,d2phipdpdt, * d2phipdp2 integer icomp,icp,imin,ioxide,itemp character answer*1 C MOLAR HEAT OF END-MEMBER COMPONENTS data cpliq/70.823,109.2,170.3,240.9,240.2,271.,172.4,180.2, * 217.,0./ C PARTIAL OF MOLAR VOLUME WITH P AT STP FOR END-MEMBER COMPONENTS data dvdp0/-1.89e-05,-2.31e-05,-2.26e-05,-2.53e-05,-2.79e-05, * -1.35e-05,-1.55e-05,-4.29e-05,-6.4e-05/ C SECOND PARTIAL OF MOLAR VOLUME WITH P FOR END-MEMBER COMPONENTS data d2vdp2/3.6e-10,5.0e-10,4.0e-10,4.4e-10,14.6e-10,4.14e-10, * 3.89e-10,8.49e-10,12.1e-10/ C SECOND PARTIAL OF MOLAR VOLUME WITH T & P FOR END-MEMBER COMPONENTS data d2vdpdt/1.3e-08,0.,2.7e-08,3.1e-08,-2.3e-08,-1.3e-08, * -1.6e-08,-5.3e-08,-4.6e-08/ C PARTIAL OF MOLAR VOLUME WITH TEMPERATURE FOR END-MEMBER COMPONENTS data dvdt0/0.,7.25e-04,2.62e-04,9.09e-04,5.84e-04,5.24e-04, * 2.92e-04,7.41e-04,7.27e-04/ c data gfusion/-.10431E+07,-.11662E+07,-.20978E+07,-.12459E+07, c * -.19011E+07,-.27549E+07,-.19749E+07,-.18546E+07, c * -.27205E+07/ C ENTHALPY OF FUSION OF END-MEMBER COMPONENTS data hfusion/-.82346E+06,-.762692E+06,-.130972E+07, * -.47792E+06,-.11763E+07,-.172335E+07,-.138822E+07, * -.13461E+07,-.176486E+07/ C MASS PER MOLE OF OXIDES data massox/.060084,.101961,.159691,.071846,.040304, * .056077,.079878,.061979,.094195,.018015/ C MASS PER MOLE OF END-MEMBER COMPONENTS data mwtcomp/.060088,.079898,.101957,.159697,.203786,.14071, * .116167,.122067,.142169,.01802/ C MOLAR ENTROPY OF FUSION OF END-MEMBER COMPONENTS data sfusion/.14842E+03,.215792E+03,.33971E+03,.40528E+03, * .48645E+03,.47691E+03,.322874E+03,.37362E+03,.47240E+03/ C TEMPERATURE OF FUSION OF END-MEMBER COMPONENTS data tfusion/1480.,1870.,2320.,1895.,1490.,2163.,1817., * 1361.,2023./ C MOLAR VOLUME OF END-MEMBER COMPONENTS data vo/2.690,2.3161,3.711,4.2131,5.4201,4.9801, * 4.347,5.5681,6.8376/ C BINARY INTERACTION COEFFICIENTS data w/ 0.0000D+00, 2.6267D+07, -3.9120D+07, 8.1100D+06, * 2.3661D+07, 3.4210D+06, -8.6400D+05, -9.9039D+07, * -3.3922D+07, 3.0976D+07, * 2.6267D+07, 0.0000D+00, -2.9450D+07, -8.7457D+07, * 5.2090D+06, -4.1780D+06, -3.5373D+07, -1.5416D+07, * -4.8095D+07, 8.1879D+07, * -3.9120D+07, -2.9450D+07, 0.0000D+00, -1.7089D+07, * -3.0509D+07, -3.2880D+07, -5.7918D+07, -1.3079D+08, * -2.3859D+07, -1.6098D+07, * 8.1100D+06, -8.7457D+07, -1.7089D+07, 0.0000D+00, * -1.7907D+08, -7.1519D+07, 1.2077D+07, -1.4966D+08, * 5.7556D+07, 3.1406D+07, * 2.3661D+07, 5.2090D+06, -3.0509D+07, -1.7907D+08, * 0.0000D+00, -3.7257D+07, -1.2971D+07, -9.0534D+07, * 2.3649D+07, 2.8874D+07, * 3.4210D+06, -4.1780D+06, -3.2880D+07, -7.1519D+07, * -3.7257D+07, 0.0000D+00, -3.1732D+07, -4.1877D+07, * 2.2323D+07, 3.5634D+07, * -8.6400D+05, -3.5373D+07, -5.7918D+07, 1.2077D+07, * -1.2971D+07, -3.1732D+07, 0.0000D+00, -1.3247D+07, * 1.7111D+07, 2.0375D+07, * -9.9039D+07, -1.5416D+07, -1.3079D+08, -1.4966D+08, * -9.0534D+07, -4.1877D+07, -1.3247D+07, 0.0000D+00, * 6.5230D+06, -9.6938D+07, * -3.3922D+07, -4.8095D+07, -2.3859D+07, 5.7556D+07, * 2.3649D+07, 2.2323D+07, 1.7111D+07, 6.5230D+06, * 0.0000D+00, 1.0374D+07, * 3.0976D+07, 8.1879D+07, -1.6098D+07, 3.1406D+07, * 2.8874D+07, 3.5634D+07, 2.0375D+07, -9.6938D+07, * 1.0374D+07, 0.0000D+00/ C MAKE SURE OXIDES ADD UP TO 100%: QUERY USER IF THEY DON'T totox = SiO2+Al2O3+Fe2O3+FeO+MgO+CaO+TiO2+Na2O+K2O if(abs(totox-100.).gt.0.005) then write (6,584) 584 format (5x,'total of component oxides does not equal 100%.', * ' Adjust automatically? (y/n)'$) read (5,585) answer 585 format (a1) if (answer.ne.'y') stop SiO2 = SiO2*(100./totox) Al2O3 = Al2O3*(100./totox) Fe2O3 = Fe2O3*(100./totox) FeO = FeO*(100./totox) MgO = MgO*(100./totox) CaO = CaO*(100./totox) TiO2 = TiO2*(100./totox) Na2O = Na2O*(100./totox) K2O = K2O*(100./totox) end if C DEFINE CONSTANTS r = 8.314 tr=298.15 pr=1.013 dtmp = 25. pbars = pMPa*10. C ASSIGN OXIDES TO ARRAY wtoxan(1) = SiO2/100. wtoxan(2) = Al2O3/100. wtoxan(3) = Fe2O3/100. wtoxan(4) = FeO/100. wtoxan(5) = MgO/100. wtoxan(6) = CaO/100. wtoxan(7) = TiO2/100. wtoxan(8) = Na2O/100. wtoxan(9) = K2O/100. C CALCULATE APPARENT ENTROPY, ENTHALPY, GIBBS FREE ENERGY FOR MELT do 100 imin=1,9 sapp(imin) = sfusion(imin) + cpliq(imin)*dlog(tabs/tfusion(imin)) happ(imin) = hfusion(imin) + cpliq(imin)*(tabs-tfusion(imin)) gapp(imin) = happ(imin) - tabs*sapp(imin) C CALCULATE SPECIFIC VOLUME AT P=1 ATM, T=tabs v1atm = vo(imin)+dvdt0(imin)*(tabs-1673.) C CALCULATE THERMODYNAMIC PROPERTIES AT pbars, tabs vdp = v1atm*(pbars-pr) + * (dvdp0(imin) + d2vdpdt(imin)*(tabs-1673.))* * (pbars**2/2. - pbars*pr + pr**2/2.) * + d2vdp2(imin)* * (pbars**3/6. - pbars*pr**2 * + pbars**2*pr - pr**3/6.) v(imin) = v1atm + (dvdp0(imin) * + d2vdpdt(imin)*(tabs-1673.))*(pbars-pr) * + d2vdp2(imin)*(pbars**2/2. - pbars*pr + pr**2/2.) dvdt(imin) = dvdt0(imin) dvdp(imin) = dvdp0(imin) + d2vdpdt(imin)*(tabs-1673.) * + d2vdp2(imin)*(pbars-pr) sapp(imin) = sapp(imin) - dvdt0(imin)*(pbars-pr) happ(imin) = happ(imin) + vdp - tabs*dvdt0(imin)*(pbars-pr) gapp(imin) = happ(imin) - tabs*sapp(imin) 100 continue C DETERMINE THERMODYNAMIC PROPERTIES OF DISSOLVED WATER USING C HAAR ET AL. EQUATION OF STATE, FOR WATER OF THE GIVEN C TEMPERATURE AND AT 1 BAR PRESSURE call haarH2O(tabs,1.0d-01,cpg,cvg,gH2O,hH2O,sH2O,rhof,vH2O,blkgas, * dvdth2o,ch2o) C DETERMINE THE MOLAR GIBBS FREE ENERGY OF WATER AT THE SPECIFIED C TEMPERATURE AND 1 BAR PRESSURE USING THE METHOD OF ROBIE. THE C EQUATIONS BELOW WERE PROVIDED BY MARK GHIORSO. gRobie = r*tabs*(2.9147*dlog(tabs) - 9.6863d-04*tabs * + 6.8593d-08*tabs**2 + 77.8899/dsqrt(tabs) * - 28954.8/tabs - 2263.27/tabs**2 - 15.8997) dgRobiedt = r*(2.9147*dlog(tabs) - 9.6863d-04*tabs * + 6.8593d-08*tabs**2 + 77.8899/dsqrt(tabs) * - 28954.8/tabs - 2263.27/tabs**2 - 15.8997) * + r*tabs*(2.9147/tabs - 9.6863d-04 * + 2.0*6.8593d-08*tabs - 0.5*77.8899/tabs**1.5 * + 28954.8/tabs**2 + 2.0*2263.27/tabs**3) d2gRobiedt2 = 2.0*r*(2.9147/tabs - 9.6863d-04 * + 2.0*6.8593d-08*tabs - 0.5*77.8899/tabs**1.5 * + 28954.8/tabs**2 + 2.0*2263.27/tabs**3) * + r*tabs*(-2.9147/tabs**2 * + 2.0*6.8593d-08 + 0.75*77.8899/tabs**2.5 * - 2.0*28954.8/tabs**3 - 6.0*2263.27/tabs**4) d3gRobiedt3 = 3.0*r*(-2.9147/tabs**2 * + 2.0*6.8593d-08 + 0.75*77.8899/tabs**2.5 * - 2.0*28954.8/tabs**3 - 6.0*2263.27/tabs**4) * + r*tabs*(2.0*2.9147/tabs**3 * - 2.25*77.8899/tabs**3.5 * + 6.0*28954.8/tabs**4 + 24.0*2263.27/tabs**5) C CALCULATE PHIP AND ITS DERIVATIVES USING AN EQUATION LIFTED FROM c A PROGRAM BY MARK GHIORSO. phiP = (0.110/tabs + 4.432d-05 + 1.405d-07*tabs * - 2.394d-11*tabs**2)*pbars * + (7.337d-08/tabs - 1.170d-08 - 9.502d-13*tabs)*pbars**2 * + (1.876d-10/tabs + 4.586d-13)*pbars**3 * - 1.191d-14*pbars**4/tabs dphiPdt = (-0.110/tabs**2 + 1.405d-07 - 2.0*2.394d-11*tabs)*pbars * + (-7.337d-08/tabs**2 - 9.502d-13)*pbars**2 * - 1.876d-10*pbars**3/tabs**2 * + 1.191d-14*pbars**4/tabs**2 d2phiPdt2 = (2.0*0.110/(tabs**3) - 2.0*2.394d-11)*pbars * + 2.0*(7.337d-08*pbars**2 + 1.876d-10*pbars**3 * - 1.191d-14*pbars**4)/tabs**3 d3phiPdt3 = -6.0*(0.110*pbars + 7.337d-08*pbars**2 * + 1.876d-10*pbars**3 - 1.191d-14*pbars**4)/tabs**4 dphiPdp = (0.110/tabs + 4.432d-05 + 1.405d-07*tabs * - 2.394d-11*tabs**2) * + 2.*(7.337d-08/tabs - 1.170d-08 - 9.502d-13*tabs)*pbars * + 3.*(1.876d-10/tabs + 4.586d-13)*pbars**2 * - 4.*1.191d-14*pbars**3/tabs d2phiPdpdt = -0.110/tabs**2 + 1.405d-07 - 2.*2.394d-11*tabs * + 2.*(-7.337d-08/tabs**2 - 9.502d-13)*pbars * + 3.*(-1.876d-10/tabs**2)*pbars**2 * + 4.*1.191d-14*pbars**3/tabs**2 d2phiPdp2 = 2.*(7.337d-08/tabs - 1.170d-08 - 9.502d-13*tabs) * + 6.*(1.876d-10/tabs + 4.586d-13)*pbars * - 12.*1.191d-14*pbars**2/tabs v(10) = dphiPdp*r*tabs dvdp(10) = d2phiPdp2*r*tabs dvdt(10) = d2phiPdpdt*r*tabs + r*dphiPdp C COMPUTE GIBBS FREE ENERGY OF WATER IN MELT USING METHODS OF GHIORSO C ET AL. (1983, CONTRIB. MIN. PET., 84:107-145). GHIORSO ET AL. C CALCULATE THE GIBBS FREE ENERGY AT TEMP=T, P=1 ATM USING THE C RELATION GAPP=(A/T + B)*RT, WHERE A AND B ARE CONSTANTS WITH VALUES C A=-33676.0 K**-1, AND B=18.3527 (FROM MARK GHIORSO). C THE GIBBS FREE ENERGY AT ELEVATED PRESSURE IS CALCULATED USING C A FORMULA FROM MARK GHIORSO ah2o = -33676.0 bh2o = 18.3527 gapp(10) = (ah2o/tabs + bh2o + phiP)*r*tabs + gh2o - grobie dgdt = r*(ah2o/tabs + bh2o + phiP) * + r*tabs*(-ah2o/tabs**2 + dphiPdt) - dgRobiedt d2gdt2 = 2.0*r*(-ah2o/tabs**2 + dphiPdt) * + r*tabs*(2.0*ah2o/tabs**3 + d2phiPdt2) - d2gRobiedt2 d3gdt3 = 3.0*r*(2.0*ah2o/tabs**3 + d2phiPdt2) * + r*tabs*(-6.0*ah2o/tabs**4 + d3phiPdt3) - d3gRobiedt3 C THE FOLLOWING FORMULAS FOR ENTROPY AND ENTHALPY OF WATER WERE TAKEN C FROM A PROGRAM BY MARK GHIORSO. sapp(10) = sh2o - dgdt happ(10) = gapp(10) + tabs*sapp(10) cpliq(10) = cpg*0.0180152 - tabs*d2gdt2 C*************************************************************************** C CALCULATE THERMODYNAMIC PROPERTIES OF THE MELT, AND PARTIAL C MOLAR PROPERTIES OF EACH COMPONENT C*************************************************************************** wtoxide(10) = dh2o/100. C CALCULATE WT% OXIDES IN HYDROUS MELT do 55 ioxide=1,9 55 wtoxide(ioxide) = 100.*wtoxan(ioxide)/(100.*(1+wtoxide(10))) C CALCULATE MOLE FRACTIONS OF OXIDES moles=0. do 200 ioxide=1,10 200 moles = moles + wtoxide(ioxide)/massox(ioxide) do 300 ioxide=1,10 xoxide(ioxide) = wtoxide(ioxide)/(massox(ioxide)* * moles) 300 continue C CALCULATE MOLES OF END-MEMBER COMPONENTS mcomp(1) = xoxide(1)-0.5*(xoxide(4)+xoxide(5)) - xoxide(6) * -xoxide(8) - 2.*xoxide(9) mcomp(2) = xoxide(7) mcomp(3) = xoxide(2)-xoxide(9) mcomp(4) = xoxide(3) mcomp(5) = 0.5*xoxide(4) mcomp(6) = 0.5*xoxide(5) mcomp(7) = xoxide(6) mcomp(8) = xoxide(8) mcomp(9) = 2.*xoxide(9) mcomp(10) = xoxide(10) C CALCULATE MOLE FRACTIONS (XCOMP) OF END-MEMBER COMPONENTS, C AVERAGE MOLECULAR WEIGHT (MWTTOT) moles=0. do 400 icomp=1,10 400 moles = moles + mcomp(icomp) mwttot = 0. do 500 icomp=1,10 xcomp(icomp) = mcomp(icomp)/moles 500 mwttot = mwttot + xcomp(icomp)*mwtcomp(icomp) C CALCULATE THERMODYNAMIC PROPERTIES OF OVERALL MELT cpm = 0.d00 sumgibbs = 0.d00 sums = 0.d00 sumxcomp = 0.d00 sumw = 0.d00 sumwx = 0.d00 sumwh = 0.d00 sumh = 0.d00 vtot = 0.d00 dvmdt = 0.d00 dvmdp = 0.d00 do 550 icomp=1,10 550 sumwxx(icomp)=0. do 600 icomp=1,10 cpm = cpm + xcomp(icomp)*cpliq(icomp) sumgibbs = sumgibbs + xcomp(icomp)*gapp(icomp) sums = sums + xcomp(icomp)*sapp(icomp) sumwx = sumwx + .001*w(10,icomp)*xcomp(10)*xcomp(icomp) sumwh = sumwh + .001*w(10,icomp)*xcomp(icomp) sumh = sumh + xcomp(icomp)*happ(icomp) vtot = vtot + xcomp(icomp)*v(icomp) dvmdt = dvmdt + xcomp(icomp)*dvdt(icomp) dvmdp = dvmdp + xcomp(icomp)*dvdp(icomp) if (xcomp(icomp).gt.0.) then sumxcomp = sumxcomp + xcomp(icomp)*log(xcomp(icomp)) end if do 650 icp = 1,10 sumw = sumw + .001*w(icomp,icp)*xcomp(icomp)*xcomp(icp) sumwxx(icomp) = sumwxx(icomp) + .001*w(icomp,icp)*xcomp(icp) 650 continue 600 continue enthalpy = sumh + 0.5*sumw if (xcomp(10).gt.0.) then entropy = sums - r*sumxcomp * - r*(xcomp(10) *log(xcomp(10)) * + (1.-xcomp(10))*log(1.-xcomp(10))) swm = sapp(10) - 2.*r*log(xcomp(10)) hwm = happ(10) + sumwxx(10) - 0.5*sumw gwm = gapp(10) + 2.*r*tabs*log(xcomp(10)) * + sumwxx(10) - 0.5*sumw dgdxw = (-2.*sumwxx(10) + sumw * + 2.*r*tabs*(1/xcomp(10) - 1.)) * / (1.-xcomp(10)) else entropy = sums - r*sumxcomp end if C CONVERT AVERAGE MOLECULAR WEIGHT FROM G/MOLE TO KG/MOLE mwttot = mwttot C CONVERT VTOT, DVMDT FROM J/(BAR MOLE K) TO M3/(KG K) vtot = vtot/(1.0D+05*mwttot) rhom = 1./vtot dvmdt = dvmdt/(1.0d+05*mwttot) dvmdp = dvmdp/(1.0d+11*mwttot) cpm = cpm/mwttot blkmag = -1/(rhom*dvmdp) sndspeed = sqrt(blkmag/rhom) C CONVERT ENTHALPY FROM J/MOLE TO J/KG; ENTROPY FROM J/MOLE K TO J/KG K enthalpy = enthalpy/mwttot entropy = entropy/mwttot hwm = hwm/0.0180152 swm = swm/0.0180152 gwm = gwm/0.0180152 dgdxw = dgdxw/0.0180152 return end