SUBROUTINE NI_OLIVINE_CALC(TK,P,Q,R,MU_NIOL,A_NIOL,
     1    MU_NIMG,MU_NIFE)

C	*****************************************************************
C       SUBROUTINE TO CALCULATE CHEMICAL POTENTIALS FOR NI-MG-FE OLIVINES
C       M. HIRSCHMANN  AMERICAN MINERALOGIST 1991 76:1232-1249
C       *****************************************************************
C       INPUTS:
C       TK       (REAL*8)   T IN KELVINS
C       P        (REAL*8)   P IN BARS
C       Q        (REAL*8)   2* X NI2SIO4 - 1
C       R        (REAL*8)   2* X FE2SIO4 - 1
C       OUTPUTS:  
C       MU_NIOL  (REAL*8)   CHEMICAL POTENTIAL OF NI2SIO4  (2 ATOM BASIS)
C       A_NIOL   (REAL*8)   ACTIVITY OF NiSi(0.5)O(2) (1 ATOM BASIS)
C	MU_NIMG  (REAL*8)   NI-MG EXCHANGE POTENTIAL
C	MU_NIFE  (REAL*8)   NI-FE EXCHANGE POTENTIAL

C	CHEMICAL POTENTIAL IS IN JOULES/GFW
C	####### NOTE THAT CHEMICAL POTENTIAL IS ON 2 ATOM BASIS #########
C	######          AND ACTIVITY IS ON 1 ATOM BASIS         #########
C       ***************************************************************
	IMPLICIT REAL*8 (A-H,O-Z)
	IMPLICIT INTEGER(I,N)
	DIMENSION XM(2,3)
	REAL*8 MU_NIOL,MU_NIMG,MU_NIFE,MD(12)
	REAL*8 K0,K1,K2,K3,L1,L2
	EXTERNAL OL_ORDER

C	RR IS GAS CONSTANT
	RR=8.31437
C	*******************************************************************
C	ERROR TRAP FOR OUT OF BOUNDS COMPOSITIONS OR FOR NI-FREE COMPOSITIONS
	IF (DABS(Q).GT.1.0.OR.DABS(R).GT.1.0) THEN
          WRITE(*,100)
	  RETURN
	END IF
	IF (Q.EQ.-1.0) THEN
	  WRITE(*,200)
	  RETURN
	END IF
  100  FORMAT(1X,'INPUT COMPOSITION IS OUT OF RANGE.',/,
     1  1X,'Q AND R MUST BOTH BE BETWEEN -1 AND 1.')
  200  FORMAT(1X,'OLIVINE MUST CONTAIN SOME NI2SIO4.',/,
     1  1X,'Q MUST BE > -1.0.')	
C	*******************************************************************
C	SOLUTION PARAMETERS FOR NI-MG-FE OLIVINES ARE HERE PRESENTED IN
C	MODULAR FORM:

C	MD(1) = 0.25*(GX+WM1+WM2  NIFE) 
C	MD(2) = 0.25*(GX+WM1+WM2  NIMG)
C	MD(3) = 0.25*(GX+WM1+WM2  FEMG)
C	MD(4) = 0.25*(GEX-WM1+WM2 NIFE)
C	MD(5) = 0.25*(GEX-WM1+WM2 NIMG)
C	MD(6) = 0.25*(GEX-WM1+WM2 FEMG)
C	MD(7) = 0.25*(GEX+WM1-WM2 NIFE)
C	MD(8) = 0.25*(GEX+WM1-WM2 NIMG)
C	MD(9) = 0.25*(GEX+WM1-WM2 FEMG)
C	MD(10)= 0.25*(GX-WM1-WM2  NIFE)
C	MD(11)= 0.25*(GX-WM1-WM2  NIMG)
C	MD(12)= 0.25*(GX-WM1-WM2  FEMG)


	MD(1)  = 0.005500E+05
	MD(2)  = 0.005500E+05
	MD(3)  = 0.050836E+05
	MD(4)  = -.031250E+05
	MD(5)  = -.031250E+05
	MD(6)  = 0.000000E+05
	MD(7)  = -.031250E+05
	MD(8)  = -.031250E+05
	MD(9)  = 0.000000E+05
	MD(10) = -.032500E+05
	MD(11) = -.032500E+05
	MD(12) = -.000000E+05
C	*******************************************************************
C	VOLUME CORRECTIONS TO SOLUTION PARAMETERS
C	VEX   = VOLUME CONTRIBUTION TO ENERGY OF GEX NIMG AND GEX NIFE
C	VGW   = VOLUME CONTRIBUTION TO WM1, WM2 FOR NIFE AND NIMG
C	VMGFE = VOLUME CONTRIBUTION TO WM1, WM2 FOR FEMG

	VEX   = 0.349000E-01  
	VGW   = -.010000E+00
	VMGFE = 0.030000E+00

C	ADJUST PARAMETERS FOR PRESSURE:
	MD(1)=MD(1)+0.25*VGW*(P-1.0)
	MD(2)=MD(2)+0.25*VGW*(P-1.0)
	MD(3)=MD(3)+0.25*VMGFE*(P-1.0)
	MD(4)=MD(4)+0.25*VEX*(P-1.0)
	MD(5)=MD(5)+0.25*VEX*(P-1.0)
	MD(7)=MD(7)+0.25*VEX*(P-1.0)
	MD(8)=MD(8)+0.25*VEX*(P-1.0)

C	***********************************************************************
C       CALCULATE STANDARD STATE FREE ENERGY FOR NI2SIO4 OLIVINE
C	AT THIS T AND P FROM HIRSCHMANN, 1991	

        H=-1395.3 *1000.0
        S=   128.1
        V=   4.260
        CP_T= 0.0
        CP_H= 0.0
        K0= 0.214997 D+03
        K1= -.103075 D+04
        K2= -.494453 D+07
        K3= 0.623705 D+09
        V1=  -0.61  D-06   
        V2=   0.00  D-12    
        V3=   28.422D-06    
        V4=   34.355D-10    

        CALL G_INTEGRAL(G_NIOL,P,TK,H,S,V,K0,K1,K2,K3,CP_H,CP_T,
     1                  L1,L2,V1,V2,V3,V4)

C	***********************************************************************
C       CALCULATE STANDARD STATE FREE ENERGY FOR MG2SIO4 OLIVINE
C	AT THIS T AND P FROM BERMAN, 1988

	H=   -2174.420*1000.0
        S=   94.010
	V=   4.366

	CP_T= 0.0
	CP_H= 0.0           

	K0=  238.64          
	K1=  -20.013 D+02  
	K2=   0.0    D+05   
	K3=  -11.624 D+07   
	L1=   0.0    D-02   
	L2=   0.0    D-05   

	V1=   -.791  D-06
	V2=   1.351  D-12
	V3=  29.464  D-06
	V4=  88.633  D-10

        CALL G_INTEGRAL(G_FO,P,TK,H,S,V,K0,K1,K2,K3,CP_H,CP_T,
     1                  L1,L2,V1,V2,V3,V4)

C	***********************************************************************
C       CALCULATE STANDARD STATE FREE ENERGY FOR FE2SIO4 OLIVINE
C	AT THIS T AND P FROM BERMAN, 1988

	H= -1479360. 
	S=  150.930                  
	V=  4.630

 	K0=     248.93           
	K1=     -19.239      D+02  
	K2=     0.0
	K3=   -13.910        D+07  
	CP_T=   0.0
	CP_H=   0.0
	L1  =   0.0
	L2  =   0.0

	V1  =  -0.730        D-06
	V2  =   0.0
	V3  =   26.546       D-06
	V4  =   79.482       D-10

        CALL G_INTEGRAL(G_FA,P,TK,H,S,V,K0,K1,K2,K3,CP_H,CP_T,
     1                  L1,L2,V1,V2,V3,V4)

C	********************************************************************
C       CALCULATE EQUILIBRIUM ORDERING STATE OF OLIVINE:

	CALL OL_ORDER(TK,Q,R,MD,S,T,XM)

C       OUTPUT 0F THIS:
C       S,T    = ORDERING PARAMETERS
C       XM(1,1)=X M1 NI
C       XM(1,2)=X M1 MG
C       XM(1,3)=X M1 FE
C       XM(2,1)=X M2 NI
C       XM(2,2)=X M2 MG
C       XM(2,3)=X M2 FE            
            
      MU_NIOL=G_NIOL+RR*TK*DLOG(XM(1,1)*XM(2,1))+MD(1)
     1 +Q*(MD(3)-MD(2)-MD(1))       
     2 +(Q**2)*MD(2)+R*(MD(1)-MD(2)+MD(3) +Q*(MD(3)+MD(2)-MD(1))        
     3 )+(R**2)*MD(3)+S*(MD(4)-MD(5)-MD(9) +Q*(MD(6)+MD(5)-MD(4)       
     4                                     )+R*(MD(6)-MD(9)))
     5 -(S**2)*MD(12)+T*(MD(9)-MD(7)+MD(5)+Q*(MD(8)-MD(5))
     6                                     +R*(MD(9)+MD(8)-MD(7))
     7                                    +S*(MD(12)+MD(11)-MD(10)))       
     8  +(T**2)*MD(11)

C	TWO ATOM BASIS ACTIVITY:	
      A_NIOL=DEXP((MU_NIOL-G_NIOL)/(RR*TK))
C	ONE ATOM BASIS ACTIVITY:
      A_NIOL=DSQRT(A_NIOL)


C 	## PLEASE NOTE, THERE IS AN ERROR IN EQUATION A2 OF HIRSCHMANN, 1991
C	## THE TERM (1+R-2Q)/4 SHOULD BE (R-1-2Q)/4


	XFO=(-Q-R)/2.0
	IF (XFO.NE.0.0) THEN
      MU_NIMG=0.5*(G_NIOL-G_FO+
     1             RR*TK*DLOG(XM(1,1)*XM(2,1)/XM(1,2)/XM(2,2)))
     2  -MD(3)-MD(2)+MD(1)-2.0*MD(2)*Q-(MD(3)-MD(2)+MD(1))*R
     3  -(MD(6)-MD(5)+MD(4))*S-(MD(8)+MD(5))*T
	ELSE
	MU_NIMG=0.0
	ENDIF

C 	## PLEASE NOTE, THERE IS AN ERROR IN EQUATION A3 OF HIRSCHMANN, 1991
C	## THE TERM (R-Q)/2 SHOULD BE (R-Q)/4

	XFA=(R+1.0)/2.0
	IF (XFA.NE.0.0) THEN
      MU_NIFE=0.5*(G_NIOL-G_FA)+
     1             RR*TK*DLOG(XM(1,1)*XM(2,1)/XM(1,3)/XM(2,3))
     2   +Q*(MD(3)-MD(2)-MD(1))
     3   +R*(MD(3)-MD(2)+MD(1))
     4   +S*(MD(4)-MD(5)-MD(9))
     5   +T*(MD(5)-MD(7)+MD(9))
	ELSE
	MU_NIFE=0.0
	END IF
	RETURN
	END

       SUBROUTINE G_INTEGRAL(G_0,P,T,HS,SS,VS,K0,K1,K2,K3,
     1                        CP_H,CP_T,L1,L2,V1,V2,V3,V4)

C       CALCULATES STD_STATE G, FREE ENERGY OF FORMATION FROM
C	THE ELEMENTS AT THE TEMPERATURE T (KELVINS) AND PRESSURE (BARS)
C	OF INTEREST USING THE BERMAN(1988) FORM.

        IMPLICIT REAL*8(A-H,O-Z)
	IMPLICIT INTEGER(I-N)
        REAL*8 K0,K1,K2,K3,L1,L2
        TR=298.15
        PR=1.0

         HS = HS + K0*(T-TR) + 2.0*K1*(SQRT(T)-SQRT(TR))
     1        - K2*(1.0/T-1.0/TR) - 0.5*K3*(1.0/T**2-1.0/TR**2)
         SS = SS + K0*LOG(T/TR)
     1        - 2.0*K1*(1.0/SQRT(T)-1.0/SQRT(TR))
     2        - 0.5*K2*(1.0/T**2-1.0/TR**2)
     3        - (1.0/3.0)*K3*(1.0/T**3-1.0/TR**3)
         IF(CP_T.NE.0.0) THEN
            IF(T.GT.CP_T) THEN
               HS = HS + CP_H
     1              + 0.5*L1**2*(CP_T**2-TR**2)
     2              + (2.0/3.0)*L1*L2*(CP_T**3-TR**3)
     3              + 0.25*L2**2*(CP_T**4-TR**4)
              SS = SS + CP_H/CP_T
     1              + L1**2*(CP_T-TR)
     2              + L1*L2*(CP_T**2-TR**2)
     3              + (1.0/3.0)*L2**2*(CP_T**3-TR**3)
            ELSE
              HS = HS + 0.5*L1**2*(T**2-TR**2)
     1              + (2.0/3.0)*L1*L2*(T**3-TR**3)
     2              + 0.25*L2**2*(T**4-TR**4)
               SS = SS + L1**2*(T-TR)
     1              + L1*L2*(T**2-TR**2)
     2              + (1.0/3.0)*L2**2*(T**3-TR**3)
            END IF
         END IF
            VS = VS*((V1/2.0-V2)*(P**2-PR**2)+V2*(P**3-PR**3)/3.0
     1              +(1.0-V1+V2+V3*(T-TR)+V4*(T-TR)**2)*(P-PR))

      G_0 = HS - T*SS + VS
        RETURN
        END
	SUBROUTINE OL_ORDER(TTK,Q,R,MD1,S,T,XM)
	
C	SUBROUTINE TO CALCULATE THE EQUILIBRIUM ORDERING STATE
C	OF A MINERAL USING PARAMETERS OF SACK-THOMPSON SOLUTION
C	MODEL AT A GIVEN T AND COMPOSITION.

C	Marc Hirschmann March, 1991.
C	THIS SUBROUTINE WAS WRITTEN TO CALCULATE EQUILIBRIUM ORDERING
C	STATES OF NI-MG-FE OLIVINES (HIRSCHMANN, 1991 AM. MIN. 76:1232-1249)

C	INPUT:
C	TTK  = DEGREES KELVIN
C	   SACK-THOMPSON CENTROSYMMETRIC COMPOSITIONAL VARIABLES:
C	Q= 2 * XNI OL -1.0
C	R= 2 * XFE OL -1.0	

C	OUTPUT:
C	SACK-THOMPSON CENTROSYMMETRIC ORDERING VARIABLES
C	S=X(FE)M2-X(FE)M1
C	T=X(NI)M1-X(NI)M2
C
C	XM(2,3)= SITE OCCUPANCIES:

C	XM(1,1)=X M1 NI
C 	XM(1,2)=X M1 MG
C	XM(1,3)=X M1 FE
C	XM(2,1)=X M2 NI
C	XM(2,2)=X M2 MG
C	XM(2,3)=X M2 FE

	IMPLICIT REAL*8(A-H,O-Z)
	IMPLICIT INTEGER(I-N)
    
	REAL*8 JAC
	DIMENSION HESS(2,2),HESSI(2,2),JAC(2)
	DIMENSION SS(2),S0(2),SHI(2),SLO(2)
	DIMENSION XM(2,3)
	REAL*8 MD1(12),MD
	REAL*8 TEST
	INTEGER BINARY
	COMMON/PARAM/MD(12),Q1,R1,TK
	LOGICAL SINGMAT,CONVERGE
	Q1=Q
	R1=R
	TK=TTK

C	*********************************************************************
C	** CONVERT SOLUTION PARAMETERS PASSED FROM MAIN PROGRAM TO
C	** PARAMETERS PASSABLE VIA COMMON BLOCK param TO SUBROUTINE 
C	** DERIV_CALC

	DO K=1,12
	MD(K)=MD1(K)
	END DO
C	*********************************************************************
C	ERROR TRAP FOR CASE OF PURE PHASES

	XFE=(R+1.0)/2.0
	XNI=(Q+1.0)/2.0
	XMG=1.0-XFE-XNI

	IF(XFE.EQ.1.0.OR.XNI.EQ.1.0.OR.XMG.EQ.1.0) THEN
	  S=0.0
	  T=0.0
	  RETURN
	END IF
C	*********************************************************************
C	DETERMINE MAXIMUM AND MINIMUM PERMISSIBLE VALUES FOR ORDERING
C	VARIABLE FOR THIS COMPOSITION. 

	SLO(1)=-1.0
	SLO(2)=-1.0
	SHI(1)=1.0
	SHI(2)=1.0
	IF(R.GT.0.0) THEN
	  SLO(1)=R-1.0
	  SHI(1)=1.0-R
	ELSE
	  SLO(1)=-1.0-R
	  SHI(1)=1.0+R
	END IF

	IF(Q.GT.0.0) THEN
	  SLO(2)=Q-1.0
	  SHI(2)=1.0-Q
	ELSE
	  SLO(2)=-Q-1.0
	  SHI(2)=Q+1.0
	END IF

	IF (XFE.NE.0.0.AND.XNI.NE.0.0.AND.XMG.NE.0.0) THEN
C	** THIS BRANCH IS FOR A TRUE TERNARY.	
	
C	CALCULATE ORDERING FOR TERNARY SOLUTION:

C	DETERMINATION OF EQUILIBRIUM ORDERING STATE IS ITERATIVE
C	BEGIN BY GUESSING INITIAL VALUES FOR ORDERING VARIABLES
C	WHICH ARE IN VECTOR S0
C	CALCULATED ORDERING VARIABLES WILL BE IN VECTOR SS.

C	FIRST GUESS FOR ORDERING
	S0(1)=SHI(1)/2.0
	S0(2)=SHI(2)/2.0


	TOL=1.0E-06	
	CONVERGE=.FALSE.
	SINGMAT=.FALSE.
C	TOL IS TOLERANCE FOR TEST OF CONVERGENCE
	ICOUNT=0
	DO ICOUNT=1,100
	IF (.NOT.CONVERGE) THEN

C	*********************************************************************	
C	CALCULATE NEW ORDERING VARIABLES SS

C	CALCULATE JACOBIAN AND HESSIAN OF GIBBS FREE ENERGY IN
C	SUBROUTINE DERIV_CALC

	CALL DERIV_CALC(S0,JAC,HESS)

C	INVERT HESSIAN
	CALL INVERT(2,HESS,HESSI,SINGMAT)
	IF (SINGMAT) THEN
C	(ERROR MESSAGE)
	WRITE(*,200)
  200   FORMAT(1X,'SINGULARITY IN MATRIX INVERSION IN SUBROUTINE OL
     1 ORDER. FATAL ERROR.')
	STOP
	END IF

C	SOLVE FOR NEW (GUESS AT) ORDERING VARIABLES:
 
	DO I=1,2
	SS(I)=S0(I)
	DO J=1,2
	SS(I)=SS(I)-HESSI(I,J)*JAC(J)
	END DO
	END DO

C	COMPARE 

	TEST=0.0
	DO I=1,2
	TEST=TEST+DABS(SS(I)-S0(I))
	END DO
	IF (TEST.LT.TOL) THEN
	  CONVERGE=.TRUE.
	ELSE
	  S0(1)=SS(1)
	  S0(2)=SS(2)
	END IF

C	CHECK TO MAKE SURE THAT NEW GUESSES FOR S0 DO NOT VIOLATE
C	PERMISSIBLE COMPOSITIONS.

	DO I=1,2
	IF (S0(I).GT.SHI(I)) S0(I)=SHI(I)-1.0E-08
	IF (S0(I).LT.SLO(I)) S0(I)=SLO(I)+1.0E-08
	ENDDO

	END IF	
C	END OF "CONVERGE" ITERATIVE LOOP
	END DO   
C	END OF ICOUNT LOOP

	IF(.NOT.CONVERGE) THEN
	WRITE(*,100)ICOUNT
  100   FORMAT(/,1X,'ORDERING ROUTINE FAILED TO CONVERGE AFTER ',
     1 I3,' LOOPS.')
	RETURN
	END IF  	

	S=SS(1)
	T=SS(2)
	XM(1,1)=(Q+T+1.0)/2.0      
 	XM(2,1)=(Q-T+1.0)/2.0      
	XM(1,2)=(S-Q-T-R)/2.0      
 	XM(2,2)=(T-Q-R-S)/2.0     
	XM(1,3)=(R-S+1.0)/2.0      
 	XM(2,3)=(R+S+1.0)/2.0     

	ELSE  
C	***********************************************************************
C	CALCULATE ORDERING FOR A BINARY SOLUTION

	IF (XFE.EQ.0.0) THEN     
C	**NI-MG BINARY

	BINARY=1

	CALL BINARY_ORDER(BINARY,T)
	XM(1,1)=(Q+T+1.0)/2.0      
 	XM(2,1)=(Q-T+1.0)/2.0      
	XM(1,2)=1.0-XM(1,1)        
 	XM(2,2)=1.0-XM(2,1)        	
	XM(1,3)=0.0		     
	XM(2,3)=0.0		   
	S=0.0
	END IF

	IF (XMG.EQ.0.0) THEN	
C	** NI-FE BINARY
	BINARY=2
	CALL BINARY_ORDER(BINARY,S)
	XM(1,3)=(R-S+1.0)/2.0     
 	XM(2,3)=(R+S+1.0)/2.0     
	XM(1,1)=1.0-XM(1,3)       
 	XM(2,1)=1.0-XM(2,3)       
	XM(1,2)=0.0               
 	XM(2,2)=0.0               
	T=S
	END IF

	IF (XNI.EQ.0.0) THEN	
C	** FE-MG BINARY
	BINARY=3
	CALL BINARY_ORDER(BINARY,S)
	XM(1,1)=0.0                
 	XM(2,1)=0.0                
	XM(1,3)=(R-S+1.0)/2.0      
 	XM(2,3)=(R+S+1.0)/2.0      
	XM(1,2)=1.0-XM(1,3)        
 	XM(2,2)=1.0-XM(2,3)        
	T=0.0
	END IF

	ENDIF
C	END OF BINARY VS. TERNARY IF-THEN

	RETURN
	END

	SUBROUTINE DERIV_CALC(SS,JAC,HESS)
	
C	** CALCULATES JACOBIAN AND HESSIAN OF GIBBS FREE ENERGY
C	WITH RESPECT TO ORDERING VARIABLES (SS).
C       FOR OLIVINE TERNARY ORDERING STATE.
	
	IMPLICIT REAL*8 (A-H,O-Z)
	IMPLICIT INTEGER(I-N)
      
	REAL*8 MD,JAC
	DIMENSION HESS(2,2),JAC(2),SS(2)
	COMMON/PARAM/MD(12),Q1,R1,TK
C	** RR = GAS CONSTANT

	RR=8.31437	
	Q=Q1
	R=R1
	S=SS(1)
	T=SS(2)

	X2MG=  (T-S-R-Q)/2.0
	X1MG= (-T+S-R-Q)/2.0
	X1NI= ( T+Q+1.0)/2.0
	X2NI= (-T+Q+1.0)/2.0
	X2FE=  (S+R+1.0)/2.0
	X1FE= (-S+R+1.0)/2.0

	IF (X2MG.LT.1.0E-12) X2MG = 1.0E-12
	IF (X1MG.LT.1.0E-12) X1MG = 1.0E-12
	IF (X2NI.LT.1.0E-12) X2NI = 1.0E-12
	IF (X1NI.LT.1.0E-12) X1NI = 1.0E-12
	IF (X2FE.LT.1.0E-12) X2FE = 1.0E-12
	IF (X1FE.LT.1.0E-12) X1FE = 1.0E-12

      JAC(1)=0.5*RR*TK*DLOG(X1MG*X2FE/X2MG/X1FE)+(-1.0*MD(12)-MD(11)
     1 +MD(10))*T+2.0*MD(12)*S+(MD(9)-MD(6))*R+(MD(4)-MD(5)-MD(6))
     2 *Q+MD(9)+MD(4)-MD(5)
      JAC(2)=0.5*RR*TK*DLOG(X2MG*X1NI/X1MG/X2NI)+2.0*MD(11)*T+(MD(10)-
     1 MD(11)-MD(12))*S+(MD(7)-MD(8)-MD(9))*R+(MD(5)-MD(8))*Q-MD(9)+
     2 MD(7)+MD(5)

      HESS(1,1) = 0.5*RR*TK*(1.0/X1MG+1.0/X2MG+1.0/X2FE+1.0/X1FE)
     1  +2.0*MD(12) 
      HESS(1,2) = -0.5*RR*TK*(1.0/X2MG+1.0/X1MG)-MD(12)-MD(11)+MD(10)
      HESS(2,1) = -0.5*RR*TK*(1.0/X2MG+1.0/X1MG)-MD(12)-MD(11)+MD(10)
      HESS(2,2) = 0.5*RR*TK*(1.0/X1MG+1.0/X2MG+1.0/X2NI+1.0/X1NI)
     1  +2.0*MD(11)

	RETURN
	END

	SUBROUTINE BINARY_ORDER(BINARY,S)

C	CALCULATES EQUILIBRIUM ORDERING STATE FOR BINARY OLIVINE SOLUTIONS
C	MMH MARCH, 1991

C 	INPUTS: 
C	BINARY: INDICATES WHICH OLIVINE BINARY
C	  BINARY=1  NI-MG OLIVINE
C	  BINARY=2  FE-NI OLIVINE
C	  BINARY=3  FE-MG OLIVINE

C	OUTPUT
C	S:  VALUE OF ORDERING VARIABLE

	IMPLICIT REAL*8(A-H,O-Z)
	IMPLICIT INTEGER(I-N)

	REAL*8 MD,K1,K2
	REAL*8 QR,VAR1,VAR2,VAR3,S,SMAX,SMIN,S1A,ROOT,A,B,C,
     1  TEST1,TEST2,DIFF
	INTEGER BINARY,NCOUNT
	COMMON/PARAM/MD(12),Q1,R1,TK
	LOGICAL BCONVERGE,PASS,SA,SB
	TOLERANCE=1.0E-8
	RR=8.31437

C	*********************************************************************
	IF(BINARY.EQ.1) THEN
	QR=Q1
	VAR1=MD(11)
	VAR2=MD(5)-MD(8)
	VAR3=MD(5)+MD(8)
	ELSE IF (BINARY.EQ.2) THEN
	QR=R1
	VAR1=MD(10)
	VAR2=MD(7)-MD(4)
	VAR3=MD(7)+MD(4)
	ELSE IF (BINARY.EQ.3) THEN
	QR=R1
	VAR1=MD(12)
	VAR2=MD(9)-MD(6)
	VAR3=MD(9)+MD(6)
	END IF

C	*********************************************************************
C	MAKE SURE THAT THERE IS SOME ENERGY OF ORDERING ON THIS BINARY
	IF (VAR1.EQ.0.0.AND.VAR2.EQ.0.0.AND.VAR3.EQ.0.0) THEN
	S=0.0
	RETURN
	END IF
C	*********************************************************************
C	ESTABLISH PERMISSIBLE LIMITS ON VALUE OF ORDERING VARIABLE FOR
C	THIS COMPOSITION:
	SMAX=1.0
	IF (SMAX.GT.(1.0-QR)) SMAX = (1.0-QR)
	IF (SMAX.GT.QR+1.0) SMAX=QR+1.0
	SMIN=-1.0
	IF (SMIN.LT.(-QR-1.0)) SMIN=-QR-1.0
	IF (SMIN.LT.(QR-1.0)) SMIN=QR-1.0
	IF (SMIN.GE.SMAX) THEN 
	WRITE(*,20)Q1,R1,SMIN,SMAX
   20   FORMAT(1X,'ERROR IN BINARY OLIVINE ORDERING ROUTINE',/,
     1  1X,'SMIN IS NOT LESS THAN SMAX! ',4(E12.6,1X))
	STOP
	END IF   !SMIN.GE.SMAX
C	*********************************************************************
C	FIRST GUESS FOR S:
	S=SMAX/2.0
C	*********************************************************************
C	ITERATIVE LOOP:
	BCONVERGE=.FALSE.
	NCOUNT=0
	DO NCOUNT=1,200
	IF (.NOT.BCONVERGE) THEN	
	K1=2.0*VAR1*S+VAR2*QR+VAR3
	K2=DEXP(-2.0*K1/(RR*TK))

	A=1.0-K2
	B=2.0+2.0*K2
	C=(K2-1.0)*(QR**2-1.0)

	ROOT=B**2-4.0*A*C
	IF (ROOT.LT.0.0) THEN
	WRITE(*,1001)
 1001   FORMAT(1X,'IMAGINARY ROOT IN OLIVINE ORDERING ROUTINE.  STOP.')
        STOP
	END IF
	S1A=(-B+DSQRT(ROOT))/(2.0*A)
	S1B=(-B-DSQRT(ROOT))/(2.0*A)

C	THIS ALGORITHM SELECTS WHICH QUADRATIC ROOT TO USE:
	PASS=.FALSE.
	DO WHILE (.NOT.PASS)
	  SA=.TRUE.
	  SB=.TRUE.
	  IF (S1A.GT.SMAX)  SA=.FALSE.
	  IF (S1B.GT.SMAX)  SB=.FALSE.
	  IF (S1A.LT.SMIN) SA=.FALSE.
	  IF (S1B.LT.SMIN) SB=.FALSE.
	  IF (SA.OR.SB) THEN
	    PASS=.TRUE.
	  ELSE	
	    IF (.NOT.SA) S1A=(S+S1A)/2.0
	    IF (.NOT.SB) S1B=(S+S1B)/2.0
	  END IF
	END DO
	
	IF(SA.AND.SB) THEN	
	  TEST1=DABS(S1A-S)
	  TEST2=DABS(S1B-S)
	   IF (TEST1.LT.TEST2) THEN
	     S1=S1A
	   ELSE
	     S1=S1B
	   ENDIF
	ELSE 
	   IF(SA) S1=S1A
	   IF(SB) S1=S1B
	END IF

C	** TEST FOR CONVERGENCE:
	DIFF=DABS(S1-S)
	IF(DIFF.LT.TOLERANCE) THEN
	  BCONVERGE=.TRUE.
	ELSE
	  S=S1
	END IF

	END IF  
C	(END OF CONVERGE IF-THEN)

	END DO 
C	(END OF NCOUNT DO LOOP)
	IF (.NOT.BCONVERGE) THEN
	WRITE(*,30)DIFF
   30   FORMAT(1X,'FAILURE TO CONVERGE IN OLIVINE ORDERING ROUTINE.',
     1  /,1X,'DIFFERENCE= ',E12.6)
	END IF
C	*********************************************************************
	RETURN
	END


	SUBROUTINE INVERT(N,M,INV_M,SINGMAT)
C	******************************************************************
C	THIS SUBROUTINE INVERTS A SYMMETRIC POSITIVE DEFINITE MATRIX
C	IT MUST BE 10X10 OR SMALLER
C	##  IT WILL NOT TEST MATRIX FOR POSITIVE DEFINITENESS     ##
C	##  IT WILL NOT TEST MATRIX FOR SYMMETRY!!!              ##
C	
C	SUBROUTINE USES BAUER-REINSCH INVERSION
C	ALGORITHM TRANSLATED FROM TURBO-PASCAL ALOGORITHM A9
C	NASH, J.C. 1990 COMPACT METHODS FOR COMPUTERS P. 99-100.
C	2ND EDITION.
C
C	TRANSLATED BY MARC HIRSCHMANN, MARCH 1991
C	***************************************************************
C	INPUTS:
C	N  - ORDER OF NXN SQUARE MATRIX
C	M  - AN NXN SYMMETRIC POSITIVE-DEFINITE MATRIX N<=10
C
C	OUTPUTS:
C	INV_M -INVERSE OF M
C	SINGMAT- SINGULAR MATRIX ERROR MESSAGE

	INTEGER I,J,KK,MM,N,Q,QQ,NVECTOR
	REAL*8 S,T,INV_M,M
	DIMENSION M(N,N),INV_M(N,N)

	REAL*8 AVECTOR(100),X(100)
	LOGICAL SINGMAT

C	READ M INTO VECTOR FORMAT LOWER TRIANGULAR MATRIX AVECTOR
C
C	E.G.	1000  SEQUENCE IS ORDER IN AVECTOR
C		2300
C		4560
C		7890

	NVECTOR=0
	DO I=1,N
	DO J=1,I
	NVECTOR=NVECTOR+1
	AVECTOR(NVECTOR)=M(J,I)
	END DO
	END DO

C	BEGIN MATRIX INVERSION
C	AT COMPLETION, AVECTOR WILL CONTAIN ELEMENTS OF INV_M

	SINGMAT=.FALSE.

	DO KK=1,N
	K=N-KK+1	 
	IF (.NOT.SINGMAT) THEN
	S=AVECTOR(1)
	IF (S.GT.0.0) THEN
	
	MM=1
	DO I=2,N
	
	Q=MM
	MM=MM+I
	T=AVECTOR(Q+1)
	X(I)=-T/S
	IF (I.GT.K) X(I)=-X(I)
	QQ=Q+2

	DO J=QQ,MM
	AVECTOR(J-I)=AVECTOR(J)+T*X(J-Q)
	END DO

	END DO

	Q=Q-1
	AVECTOR(MM)=1.0/S
	
	DO I=2,N
	AVECTOR(Q+I)=X(I)
	END DO
	END IF
	
	ELSE
	SINGMAT=.TRUE.
	END IF
	END DO

C	END OF MATRIX INVERSION
C 	NOW READ INVERTED MATRIX INTO INV_M

	NVECTOR=0
	DO I=1,N
	DO J=1,I
	NVECTOR=NVECTOR+1
	INV_M(J,I)=AVECTOR(NVECTOR)
	IF (J.NE.I) INV_M(I,J)=AVECTOR(NVECTOR)
	END DO
	END DO
	RETURN
	END