C  ********************************************************
C    "NIGGLI.FOR".   WRITTEN IN FORTRAN 77.
C    This user friendly program calculates the NIGGLI values
C    from the composition of a rock, expressed in weight
C    percent of constituent oxides.
C    Then, it performs a search through a database called
C    "MAGMA.DAT" in order to determine the corresponding
C    Magma Type of the rock.
C    Dr. G. Tsambourakis & S.F. Youl
C  ********************************************************
      INTEGER   IN,OUT,REF,RESULT,OXOUT(27), NOXOUT
      PARAMETER ( IN=5,OUT=6,REF=9,RESULT=10 )
      REAL MWV,MG,K
      CHARACTER*1 ANSWER
      CHARACTER*5 OXIDE(27),TITLEX
      CHARACTER*30 NAME
      DIMENSION MWV(27),WTP(27),V(20),R1(27),R2(27)
C   OXIDE - Oxide, MWV - MOLecular Weight Values.
      DATA OXIDE(1),MWV(1)/'SiO2 ',60.084 /
      DATA OXIDE(2),MWV(2)/'TiO2 ',79.86  /
      DATA OXIDE(3),MWV(3)/'Al2O3',101.96 /
      DATA OXIDE(4),MWV(4)/'Fe2O3',159.69 /
      DATA OXIDE(5),MWV(5)/'FeO  ',71.846 /
      DATA OXIDE(6),MWV(6)/'MnO  ',70.9374/
      DATA OXIDE(7),MWV(7)/'NiO  ',74.69  /
      DATA OXIDE(8),MWV(8)/'MgO  ',40.304 /
      DATA OXIDE(9),MWV(9)/'CaO  ',56.077 /
      DATA OXIDE(10),MWV(10)/'BaO  ',153.326/
      DATA OXIDE(11),MWV(11)/'SrO  ',103.6  /
      DATA OXIDE(12),MWV(12)/'Na2O ',61.9798/
      DATA OXIDE(13),MWV(13)/'K2O  ',94.196 /
      DATA OXIDE(14),MWV(14)/'P2O5 ',141.944/
      DATA OXIDE(15),MWV(15)/'F-   ',18.9984/
      DATA OXIDE(16),MWV(16)/'S    ',32.066 /
      DATA OXIDE(17),MWV(17)/'CO2  ',44.009 /
      DATA OXIDE(18),MWV(18)/'Cr2O3',151.99 /
      DATA OXIDE(19),MWV(19)/'ZnO  ',81.38  /
      DATA OXIDE(20),MWV(20)/'Li2O ',29.881 /
      DATA OXIDE(21),MWV(21)/'SO3  ',80.063 /
      DATA OXIDE(22),MWV(22)/'H2O+ ',18.0153/
      DATA OXIDE(23),MWV(23)/'H2O- ',18.0153/
      DATA OXIDE(24),MWV(24)/'V2O3 ',149.88 /
      DATA OXIDE(25),MWV(25)/'CoO  ',74.932 /
      DATA OXIDE(26),MWV(26)/'Rb2O ',186.935/
      DATA OXIDE(27),MWV(27)/'Cs2O ',281.8  /
      SUM1=0.0
      SUM2=0.0
      WRITE(OUT,20)
20    FORMAT(/' Sample No: ? ')
      READ(IN,30)NAME
30    FORMAT(A)
      WRITE(OUT,40)
40    FORMAT(/' How many elements have been analysed: ?')
      READ(IN,50)NUM
50    FORMAT(I2)
      WRITE(OUT,60)
60    FORMAT(/' Type oxides and Wt-%(A6,F8.3)'/)
      DO 4 IJ=1,NUM
1     READ(IN,70)TITLEX,WTPER
70    FORMAT(A,F8.4)
      DO 2 I=1,27
      NUMCOM=I
      IF(TITLEX.EQ.OXIDE(I)) GOTO 3
2     CONTINUE
      WRITE(OUT,80)TITLEX
80    FORMAT(10X,20H No compound called  ,A,8H in data)
      WRITE(OUT,90)
90    FORMAT(/' Do you want to try again (Y/y))?'/)
      READ(IN,100) ANSWER
100   FORMAT(A)
      IF(ANSWER.EQ.'Y'.OR.ANSWER.EQ.'y') GOTO 1
      STOP
3     CONTINUE
      SUM1=SUM1+WTPER
      XMW=MWV(NUMCOM)
      RES1=WTPER/XMW
      SUMMOLS=SUMMOLS+0.5*RES1
      WTP(NUMCOM)=WTPER
4     CONTINUE
      OPEN (UNIT=RESULT, FILE='NIGGLI.OUT')
      WRITE(RESULT,120)NAME
120   FORMAT(//,15X,8H SAMPLE:,A)
      WRITE(RESULT,130)
130   FORMAT(16X,'-------------------------------'/)
      WRITE(RESULT,140)
140   FORMAT(2X,'OXIDES',4X,'wt-%',4X,'OXIDES',4X,'wt-%',4X,
     &'OXIDES',4X,'wt-%',4X,'OXIDES',4X,'wt-%')
      WRITE(RESULT,150)
150   FORMAT(1X,73('*')/)
      NOXOUT=0
      DO 5 IJ=1,27
	IF (WTP(IJ).NE.0.0) THEN
	  NOXOUT=NOXOUT+1
	  OXOUT(NOXOUT)=IJ
	  R1(IJ)=WTP(IJ)/MWV(IJ)
	  R2(IJ)=R1(IJ)*1000
	END IF
5     CONTINUE
      V(1)=R2(1)
      V(2)=R2(2)
      V(3)=R2(3)
      V(4)=2*R2(4)+R2(5)+R2(6)+R2(7)+R2(8)+2*R2(18)+2*R2(24)+R2(25)
      V(5)=2*R2(4)
      V(6)=R2(8)+0.5*R2(20)
      V(7)=2*R2(10)+2*R2(11)+R2(12)+R2(13)+R2(26)+R2(27)
      V(8)=R2(8)
      V(9)=R2(9)
      V(10)=R2(13)
      V(11)=R2(14)
      V(12)=R2(15)
      V(13)=R2(17)
      V(14)=R2(19)
      V(15)=R2(16)
      V(16)=R2(22)
      V(17)=2*R2(4)+R2(5)
      WRITE(RESULT,160)(OXIDE(OXOUT(I)),WTP(OXOUT(I)), I=1,NOXOUT)
160   FORMAT(3X,A,F9.3,4X,A,F9.3,4X,A,F9.3,4X,A,F9.3)
      WRITE(RESULT,180)SUM1
180   FORMAT(2X,6H TOTAL,F9.3,//)
      WRITE(RESULT,190)
190   FORMAT(1X,14H NIGGLI VALUES)
      WRITE(RESULT,200 )
200   FORMAT(1X,14H *************/)
      SUM2=0.0
      SUM2= V(3)+V(4)+V(7)+V(9)
      AL=V(3)/SUM2*100
      FM=V(4)/SUM2*100
      C =V(9)/SUM2*100
      ALK=V(7)/SUM2*100
      SI=V(1)/SUM2*100
      TI=V(2)/SUM2*100
      P =V(11)/SUM2*100
      H =V(16)/SUM2*100
      K =V(10)/V(7)
      MG=V(8)/V(4)
      W =V(5)/V(17)
      WRITE(RESULT,210 )
210   FORMAT(4X,' al      fm      c      alk      ')
      WRITE(RESULT,220 )AL,FM,C,ALK
220   FORMAT(4F8.2)
      WRITE(RESULT,230 )
230   FORMAT(/4X,' si      k       ti      p       h       mg      w
     &	   ')
      WRITE(RESULT,240 )SI,K,TI,P,H,MG,W
240   FORMAT(7F8.2/)
      IF (AL.LT.ALK) THEN
	  SI1=100+3*AL+ALK
	  QZ=SI-SI1
	WRITE(RESULT,250)QZ
	ELSE IF (AL.LT.(ALK+C)) THEN
	  SI1=100+4*ALK
	  QZ=SI-SI1
	WRITE(RESULT,250)QZ
	ELSE
	  SI1=100+5*ALK-AL
	  QZ=SI-SI1
	WRITE(RESULT,250)QZ
	WRITE(RESULT,255)
255     FORMAT('Excess Alumina calculated as Corundum')
      END IF
      CALL GET MAG TYPE (REF,RESULT,AL,FM,C,ALK,SI,K,MG)
      CLOSE (UNIT=RESULT)
      STOP
250   FORMAT('  "Quartz zahl" qz = ',F6.2)
      END
      SUBROUTINE GET MAG TYPE (LUN REF, LUN RES,
     &  AL,FM,C,ALK,SI,K,MG)
      INTEGER   LUN REF,LUN RES
      REAL      AL,ALK,C,FM,K,MG,SI
c
      INTEGER   MAX MNRLS, N COMP
      PARAMETER (MAX MNRLS=200, N COMP=7)
      INTEGER   MAX MAG, MAX SER
      PARAMETER (MAX MAG=50, MAX SER=5)
c
c       Variable Declarations
c
      CHARACTER MAG NAME(MAX MAG,2)*78, NAME(MAX MNRLS)*25
      CHARACTER SER NAME(MAX SER)*78
      INTEGER   I,INDX,J,MAG TYPE(MAX MNRLS), NUM MINRLS
      INTEGER   S LIST(MAX MNRLS), S NUM, SER TYPE(MAX MNRLS)
      REAL      REF COMP(MAX MNRLS,N COMP), IN COMP(N COMP)
c
      REAL      TOL1(N COMP),TOL2(N COMP),TOL3(N COMP),TOL4(N COMP),
     &          TOL5(N COMP),TOL6(N COMP), TOL7(N COMP)
      DATA TOL1 / 0.08, 0.08, 0.10, 0.12, 0.12, 0.35, 0.35 /
      DATA TOL2 / 0.10, 0.15, 0.15, 0.15, 0.15, 0.44, 0.44 /
      DATA TOL3 / 0.10, 0.18, 0.18, 0.18, 0.22, 7.50, 0.60 /
      DATA TOL4 / 0.20, 0.20, 0.20, 0.20, 0.20, 0.60, 0.70 /
      DATA TOL5 / 0.25, 0.22, 0.22, 0.30, 0.22, 0.60, 0.60 /
      DATA TOL6 / 0.27, 0.27, 0.27, 0.35, 0.35, 1.00, 1.00 /
      DATA TOL7 / 0.30, 0.30, 0.30, 0.35, 0.35, 2.00, 2.00 /
c
c       Read the mineralogical data
c
      CALL READ IN (LUN REF,NAME,MAX MNRLS,REF COMP,N COMP,
     &               SER NAME,MAX SER,MAG NAME,MAX MAG,
     &               SER TYPE,MAG TYPE,NUM MINRLS)
c
c       Set the input composition
c
      IN COMP(1) = SI
      IN COMP(2) = AL
      IN COMP(3) = FM
      IN COMP(4) = C
      IN COMP(5) = ALK
      IN COMP(6) = K
      IN COMP(7) = MG
      CALL LOCATE (IN COMP,N COMP,REF COMP,MAX MNRLS,NUM MINRLS,
     & TOL1,S LIST,S NUM)
      IF ( SNUM.EQ.0 ) THEN
      WRITE (*,*) 'SECOND SET OF DETECTION LIMITS TRIED'
      CALL LOCATE (IN COMP,N COMP,REF COMP,MAX MNRLS,NUM MINRLS,
     &  TOL2,S LIST,S NUM)
      ENDIF
      IF ( SNUM.EQ.0 ) THEN
	WRITE (*,*) 'THIRD SET OF DETECTION LIMITS TRIED'
	CALL LOCATE (IN COMP,N COMP,REF COMP,MAX MINRLS,NUM MINRLS,
     &  TOL3, S LIST, S NUM)
      ENDIF
      IF (SNUM.EQ.0) THEN
	WRITE (*,*) 'FOURTH SET OF DETECTION LIMITS TRIED'
	CALL LOCATE (IN COMP, N COMP, REF COMP, MAX MNRLS, NUM MINRLS,
     &  TOL4, S LIST, S NUM)
      ENDIF
      IF (SNUM.EQ.0) THEN
	WRITE (*,*) 'FIFTH SET OF DETECTION LIMITS TRIED'
	CALL LOCATE (IN COMP, N COMP, REF COMP, MAX MNRLS, NUM MINRLS,
     &  TOL5, S LIST, S NUM)
      ENDIF
      IF (SNUM.EQ.0) THEN
	WRITE (*,*) 'SIXTH SET OF DETECTION LIMITS TRIED'
	CALL LOCATE (IN COMP, N COMP, REF COMP, MAX MNRLS, NUM MINRLS,
     &  TOL6, S LIST, S NUM)
      ENDIF
      IF (SNUM.EQ.0) THEN
	WRITE (*,*) 'SEVENTH SET OF DETECTION LIMITS TRIED'
	CALL LOCATE (IN COMP, N COMP, REF COMP, MAX MNRLS, NUM MINRLS,
     &  TOL7, S LIST, S NUM)
      ENDIF
c     Output match/es, if any.
c
      IF (SNUM.GT.0) THEN
	WRITE(LUN RES,260)
	IF (SNUM.GT.1) WRITE(LUN RES,270)SNUM
	OLD SER=0
	OLD MAG=0
	DO I=1,SNUM
	  INDX=SLIST(I)
	  IF (OLD SER.NE.SER TYPE(INDX)) THEN
	     IF ( OLD SER .NE. 0 ) WRITE (LUN RES,310)
	     WRITE (LUN RES,280) SER NAME(SER TYPE(INDX))
	     OLD SER=SER TYPE(INDX)
	     OLD MAG = 0
	  END IF
	  IF (OLD MAG.NE.MAG TYPE(INDX)) THEN
	     IF ( OLD MAG .NE. 0 ) WRITE (LUN RES,310)
	     WRITE (LUN RES,280) MAG NAME(MAG TYPE(INDX),1)
	     IF ( MAG NAME(MAG TYPE(INDX),2) .NE. ' ' )
     &         WRITE (LUN RES,285) MAG NAME(MAG TYPE(INDX),2)
	     WRITE (LUN RES,290)
	     OLD MAG=MAG TYPE (INDX)
	  ENDIF
	  WRITE (LUN RES,300) NAME(INDX),
     &                  (REF COMP(INDX,J), J = 1,N COMP)
	END DO
	  WRITE (LUN RES,310)
      ELSE
	WRITE (LUN RES,320 )
      END IF
      STOP
260   FORMAT(//1X,' MAGMA TYPE',/,1X' **********',/)
270   FORMAT(1X,' The following ',I1,' magma types match closely',/)
280   FORMAT(/,A)
285   FORMAT(A)
290   FORMAT(/1X' Name',22X,'si  al   fm    c   alk   k    mg',/,
     &          2X,59('_'))
300   FORMAT(2X,A25,F4.0,4F5.1,2F5.2)
310   FORMAT(2X,59('-'),//)
320   FORMAT(//,'No matches found in the data base')
      END
      SUBROUTINE READ IN (LUN, NAME, MAX MNRLS, REF COMP, N COMP,
     &                     SER NAME, MAX SER, MAG NAME, MAX MAG,
     &                     SER TYPE, MAG TYPE, NUM MINRLS)
      INTEGER   LUN, MAX MAG, MAX MNRLS, MAX SER, N COMP, NUM MINRLS
      INTEGER   MAG TYPE(MAX MNRLS), SER TYPE(MAX MNRLS)
      CHARACTER MAG NAME(MAX MAG,2)*78, NAME(MAX MNRLS)*25
      CHARACTER SER NAME(MAX SER)*78
      REAL      REF COMP(MAX MNRLS,N COMP)
c
      CHARACTER REF NAME*40
      PARAMETER (REF NAME='MAGMA.DAT')
c
      CHARACTER INLINE*78
      LOGICAL   LEX
      INTEGER   I,J,MAG KT, SER KT
c
      INQUIRE (FILE=REF NAME, EXIST=LEX)
      IF (.NOT.LEX) THEN
	WRITE (*,*) 'Data file does not exist: ',REF NAME
	STOP
      END IF
      OPEN (UNIT=LUN, FILE=REF NAME, STATUS='OLD')
      I = 1
      MAG KT = 0
      SER KT = 0
10    CONTINUE
	CALL GET NXT LINE (LUN,INLINE,*20)
	IF (INLINE(1:2).EQ.'**') THEN
	  CALL GET NXT LINE (LUN,INLINE,*20)
	  SER KT = SER KT + 1
	  SER NAME(SER KT) = INLINE
	ELSE IF (INLINE(1:2).EQ.'  ') THEN
	  MAG KT = MAG KT + 1
	  MAG NAME(MAG KT,1) = INLINE
	  READ (LUN,330 ) MAG NAME(MAG KT,2)
	ELSE
	  NAME(I) = INLINE(1:25)
	  INLINE  = INLINE(26:)
	  READ (INLINE,340 ,ERR=20) (REF COMP(I,J), J = 1,N COMP)
	  MAG TYPE(I) = MAG KT
	  SER TYPE(I) = SER KT
	  I = I + 1
	END IF
      IF (I.LE.MAX MNRLS) GOTO 10
20    CONTINUE
      NUM MINRLS = I - 1
      CLOSE (UNIT=LUN)
      IF (NUM MINRLS.EQ.0) THEN
	WRITE (*,*) 'Data file is empty or corrupt: ',REF NAME
	STOP
      END IF
      RETURN
330   FORMAT (A)
340   FORMAT (F3.0,1X,4(F4.0,1X),2(F4.2,1X))
      END

      SUBROUTINE GET NXT LINE (LUN,INLINE,*)
c
c       Get the next line to be processed from LUN and return
c       it in INLINE.
c
c       Ignore blank lines and lines beginning with #
c
c       Take the first return on an error or EOF
c
      INTEGER   LUN
      CHARACTER INLINE*78
c
10    CONTINUE
	READ (LUN,350,ERR=20,END=20) INLINE
      IF (INLINE.EQ.' '.OR.INLINE(1:1).EQ.'#') GOTO 10
      RETURN
20    CONTINUE
      RETURN 1
350   FORMAT(A)
      END

      SUBROUTINE LOCATE (IN COMP, N COMP, REF COMP, MAX REF,
     &                    NUM MINRLS, TOLRNCE, S LIST, S NUM)
      INTEGER   MAX REF, N COMP, NUM MINRLS, S NUM, S LIST(S NUM)
      REAL      IN COMP(N COMP), REF COMP(MAX REF,N COMP)
      REAL      TOLRNCE(N COMP)
c
      INTEGER   I,J,K
      REAL      TOL
c
c       Function Declarations
c
      LOGICAL   COMPARE
c
c       Start searching
c
      S NUM = NUM MINRLS
      J = 0
c
c       The strategy is to compare the first component of the input
c       composition IN COMP(1) with the first component of the reference
c       data REF COMP(I,1). The index of successful comparisons is placed
c       in S LIST.
c
      DO I = 1,S NUM
	TOL = TOLRNCE(1)
	IF (COMPARE (IN COMP(1),REF COMP(I,1),TOL)) THEN
	  J = J + 1
	  S LIST(J) = I
	END IF
      END DO
c
c       S LIST is now further pruned by doing sequential comparisons on
c       the 2nd through N COMPth component of the input composition.
c
      DO K = 2, N COMP
	S NUM = J
	J = 0
	TOL = TOLRNCE(K)
	DO I = 1,SNUM
	  IF (COMPARE (IN COMP(K),REF COMP(S LIST(I),K),TOL)) THEN
	    J = J + 1
	    S LIST(J) = S LIST(I)
	  END IF
	END DO
      END DO
      S NUM = J
      RETURN
      END

      LOGICAL FUNCTION COMPARE (A,B,TOL)
      REAL      A,B,TOL
c
c       Compare A to B with a relative tolerance of TOL.
c       Zero or negative values of B indicate that the comparison
c       should be considered true.
c
      COMPARE = .TRUE.
      IF (B.LE.0.0) RETURN
      IF (ABS (1.0 - A/B).LE.TOL) RETURN
      COMPARE = .FALSE.
      RETURN
      END
     
