	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

