text:
| The Computer Language Shootout
| http://shootout.alioth.debian.org/
| Contributed by Ian Osgood
| modified by Anton Ertl
| Joshua Grams removed most of the double cell manipulation
|  and added island detection.
| Modified for Reva 6.1.5 (http://ronware.org/reva/wiki/6.1.5)
|  by Ron Aaron

: d2/ inline{ D1 F8 D1 3E } ;
: d2* inline{ D1 26 D1 D0 } ;

: enum ( n -- )  0 do i constant loop ;
: table create does> ( i -- t[i] )  swap cells + @ ;

7 enum         E   SE   SW    W   NW   NE   STOP
table offset   1 ,  7 ,  6 , -1 , -7 , -6 ,    0 ,
table rotate  SE , SW ,  W , NW , NE ,  E , STOP ,
table reflect  E , NE , NW ,  W , SW , SE , STOP ,

| paths are more easily transformable than bit masks

create path    5 cells allot
create offsets 6 cells allot

1024 allot | padding for Pentium 4 and bigforth shortcomings

: init-path ( 4*dirs -- )
	E path 5 0do tuck ! cell+ loop drop ;
: rotate-path
	path 5 0do dup @ rotate  over ! cell+ loop drop ;
: reflect-path
	path 5 0do dup @ reflect over ! cell+ loop drop ;

: path-offsets
	0 offsets !
	path offsets
	5 0do
		over @ offset
		over @ +
		over cell+ !
		swap cell+ swap cell+
	loop 2drop ;
: minimum-offset ( -- n )
	offsets @
	6 1 do offsets i cells + @ min loop ;
: normalize-offsets
	minimum-offset negate
	6 0do dup offsets i cells + +! loop drop ;
: offsets-mask ( -- mask )
	0
	6 0do
		offsets i cells + @
		1 swap << or
	loop ;

| make and store the twelve transformations of the path

: path-mask ( -- mask )  path-offsets normalize-offsets offsets-mask ;
: path-masks ( 4*dirs -- )
	                       false , | used flag
	         init-path path-mask ,
	5 0 do rotate-path path-mask , loop
	      reflect-path path-mask ,
	5 0 do rotate-path path-mask , loop ;

13 cells constant /piece

|  all paths start with an implicit E

create pieces
 STOP SE  E  E path-masks
 STOP NE  E SE path-masks
 STOP SW SE  E path-masks
 STOP SE SW  E path-masks
   SW  W  E SE path-masks | one backtrack, since this shape branches
 STOP SE NE SE path-masks
 STOP NE SE SE path-masks
 STOP  E SW SE path-masks
 STOP  E SE  E path-masks
 STOP NE SE  E path-masks

variable #solutions
create smallest 64 allot
create largest  64 allot

variable board    | high word of board; low word on stack

1024 allot | padding for Pentium 4 and bigforth shortcomings

: put-piece ( piece shift -- )
	over pieces - /piece / '0 + >r ( R: piece-char )
	here + swap @ ( buf mask )
	repeat
		dup 1 and if over r@ swap c! then
		swap 1+ swap 2/
	dup
	| 0= until
	while
	2drop r> drop ;

| extract solution from stack and store at HERE
| (non-destructive because we need the data for backtracking).
: store-solution ( pieceN shiftN ... piece0 board )
	| here 64 [char] * fill
	0 ( absolute-shift )
	depth 1- 2 swap do
		i pick over put-piece
		i 1- pick +
	-3 skip | -2 +loop
	loop
	drop ;

: .line ( line -- line+6 )
	5 0do dup c@ emit space 1+ loop cr 1+ ;
: .solution ( buffer -- )
	5 0do .line 1+  space .line loop drop cr ;

: record ( [st] -- [st] )
	store-solution  | here .solution
	here 64 smallest 64 cmp 0 <if here smallest 64 move then
	largest 64  here 64 cmp 0 <if here  largest 64 move then
	1 #solutions +! ;  | throw if #solutions == NUM

| initial board, with edges filled in
| 2 base !
| 0000011.000001.0000011.000001.0000011.000001.0000011.000001.0000011.00000
| %00000110000010000011000001000001 %10000010000011000001000001100000

| board mask for a hexagon with an empty center
%110000101000011 constant empty-hex

$80000000 constant hi-bit

| is it a single-cell island?
	| the center (empty) cell is 7 bits in.
: island? ( board bit -- flag )  empty-hex * 7 >> tuck and = ;

| fun with bit manipulation :)
: fill-leading ( u -- u' )  dup 1- or ;
: first-empty ( u -- bit )  dup dup 1+ or xor ;

| return a bit-mask for the second empty cell on the board.
: second ( board -- bit )  fill-leading first-empty ;

| check two spots for single-cell islands
: prune? ( board -- flag )
	dup 1 island? if drop true else dup second island? then ;

| remove filled cells at the beginning of the board
: (shift-board)
	0 swap board @
: (shift-again)
	over 1 and 0if ;; then
	| while true:
	d2/ hi-bit or  rot 1+ -rot
	(shift-again)
	;

: shift-board ( board -- shift board' ) (shift-board) board ! ;

| restore filled cells at the beginning of the board
: unshift-board ( shift board -- board' )
	board @ rot 0 ?do d2* swap 1+ swap loop board ! ;

: invert -1 xor ;
: solve ( board -- board )
	dup prune? if ;; then
	pieces  10 0do
		dup @ if
			/piece +
		else
			true over ! cell+        | mark used
			12 0do
				2dup @ and 0if
					tuck @ xor       | add piece
					dup invert if
						shift-board solve unshift-board
					else record then
					over @ xor swap  | remove piece
				then
			cell+ loop
			false over /piece - !    | mark unused
		then
	loop drop ;

: main
	.s cr
	#solutions off
	smallest 64 '9 fill
	largest  64 '0 fill
	.s cr
		%10000010000011000001000001100000
		%00000110000010000011000001000001
		board ! solve drop
	#solutions ? ." solutions found" cr cr
	smallest .solution
	largest  .solution ;

main bye

|||

 build & benchmark results

RUNNING SCRIPT: meteor.bigforth-4.bigforth

Wed Jan 24 20:54:26 PST 2007

=================================================================
COMMAND LINE (%A is single numeric argument):

 /opt/bigforth/bigforth $BIGFORTH_FLAGS ../meteor.bigforth-4.bigforth %A

PROGRAM OUTPUT
==============
2098 solutions found

0 0 0 0 1
 2 2 2 0 1
2 6 6 1 1
 2 6 1 5 5
8 6 5 5 5
 8 6 3 3 3
4 8 8 9 3
 4 4 8 9 3
4 7 4 7 9
 7 7 7 9 9

9 9 9 9 8
 9 6 6 8 5
6 6 8 8 5
 6 8 2 5 5
7 7 7 2 5
 7 4 7 2 0
1 4 2 2 0
 1 4 4 0 3
1 4 0 0 3
 1 1 3 3 3

max_stars_repo_name: ronaaron/reva

max_issues_repo_name: ronaaron/reva

max_forks_repo_name: ronaaron/reva


text:
      subroutine NumberModel(cmodel,model)
      character cmodel*21,cmodel2*21
      cmodel2=cmodel
      model2=model
      stop'   ***** This program can only run EPOS *****      '
      end

      subroutine IniModel(model)
      model2=model
      end

      subroutine IniEvtModel
      end

      subroutine emsaaaModel(model,id,iret)
      model2=model
      id2=id
      iret2=iret
      end

      function crseModel(model,ekin,maproj,matarg,idtarg)
      crseModel=0.
      model2=model
      ekin2=ekin
      maproj2=maproj
      matarg2=matarg
      idtarg2=idtarg
      end

      subroutine crseaaModel(sigt,sigi,sigc,sigel)
      sigt=0.
      sigi=0.
      sigc=0.
      sigel=0.
      end

      subroutine m2XXFZ( a,b)
      double precision a, b(2),c,d(2)
      c=a
      d(1)=b(1)
      d(2)=b(2)

      end

      subroutine m3SIGMA(ek,idpro,idtar,latar,matar,sigi,sige)
      sige=0.
      sigi=0.
      ek2=ek
      idpro2=idpro
      matar2=matar
      idtar2=idtar
      latar2=latar
      end

      subroutine m6SIGMA(icl,engy,stot,sela,sine,sdifr,slela,Rho)
      icl2=icl
      engy2=engy
      stot=0.
      sela=0.
      sine=0.
      sdifr=0.
      slela=0.
      Rho=0.
      end

      subroutine m7SIGMA(stot,scut,sine,slela)
      stot=0.
      scut=0.
      sine=0.
      slela=0.
      end

      subroutine m8SIGMA(stot,scut,sine,sela,slela,ssd)
      stot=0.
      scut=0.
      sine=0.
      sela=0.
      slela=0.
      ssd=0.
      end

      subroutine m9SIGMA(stot,sine,sela)
      stot=0.
      sine=0.
      sela=0.
      end

      subroutine decaymod(ip)
      idum=ip
      end

ext: f

max_stars_repo_name: kotania/impy

max_issues_repo_name: kotania/impy

max_forks_repo_name: kotania/impy


text:
c gmm01f
c  THIS FILE IS PART OF THE ORIGINAL GMM01F DISTRIBUTION
      subroutine cofd0(nmax)
      implicit double precision (a-h,o-z)
      include 'gmm01f.par'
      parameter (nmp=np*(np+2))
      parameter (ni0=np*(np+1)*(2*np+1)/3+np*np)
      integer v
      double precision lnfacd
      common/cofmnv0/cof0(ni0)
      common/crot/cofsr(nmp)
          common/fnr/fnr(0:2*(np+2))
      i=0
      sm=-0.5d0*dble((-1)**nmax)
      do m=-nmax,nmax
         ns=max(1,iabs(m))
         sm=-sm
         do n=ns,nmax
            inm=n*(n+1)-m
            do v=ns,nmax
               i=i+1
               ivm=v*(v+1)+m
               c=cofsr(inm)+cofsr(ivm)
               c=sm*dexp(c)
               c0=fnr(2*n+1)*fnr(2*v+1)
               c1=fnr(n)*fnr(v)*fnr(n+1)*fnr(v+1)
               c0=c0/c1
                   cof0(i)=c*c0
                enddo
             enddo
          enddo
      return
      end

max_stars_repo_name: m-ringler/GMM-DIP

max_issues_repo_name: m-ringler/GMM-DIP

max_forks_repo_name: m-ringler/GMM-DIP


text:
      SUBROUTINE CHECKINPUT_ST
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/CONST5/ILO,INLO,IONLYLO
      COMMON/CUT1/PTMIN,PTMAX
      COMMON/CUT2/YMIN,YMAX
      COMMON/FLAVOR/IFLAVOR,ITOTAL
      COMMON/CHARCONJ/ICHARCONJ

C***  CHECK MST1, MG, MT, MS AND MST2 ( SHOULD BE ABOVE 1 GEV )
C***  CHECK SIN2T
      IF (MST1.LT.1D0) THEN
         MST1 = 1D0
         PRINT *,'MST1 CHANGED TO 1.'
      END IF
      IF (MG.LT.1D0) THEN
         MG = 1D0
         PRINT *,'MG CHANGED TO 1.'
      END IF
      IF (MT.LT.1D0) THEN
         MT = 1D0
         PRINT *,'MT CHANGED TO 1.'
      END IF
      IF (MS.LT.1D0) THEN
         MS = 1D0
         PRINT *,'MS CHANGED TO 1.'
      END IF
      IF (MST2.LT.1D0) THEN
         MST2 = 1D0
         PRINT *,'MST2 CHANGED TO 1.'
      END IF
      IF (ABS(SIN2T).GT.1D0) THEN
         SIN2T = 0D0
         PRINT *,'SIN2T CHANGED TO 0.'
      END IF

C***  CHECK ICOLL
      IF ((ICOLL.LT.0).OR.(ICOLL.GT.1)) THEN
         PRINT *,'COLLIDER TYPE SPECIFIED INCORRECTLY; ICOLL = ',ICOLL
         PRINT *,'PROGRAM ABORTED'
         call HARD_STOP
      END IF

C***  CHECK ENERGY ( SHOULD BE ABOVE PRODUCTION THRESHOLD )
      IF (ENERGY.LT.2*MST1) THEN
         PRINT *,'ENERGY TOO SMALL TO PRODUCE STOP-ANTISTOP'
         PRINT *,'PROGRAM ABORTED'
         call HARD_STOP
      END IF

C***  CHECK IFLAVOR
      IF ((IFLAVOR.LT.0).OR.(IFLAVOR.GT.5)) THEN
         PRINT *,'IFLAVOR IN THE WRONG RANGE; IFLAVOR = ',IFLAVOR
         PRINT *,'IFLAVOR CHANGED TO 0 (ALL)'
         IFLAVOR = 0
      END IF

C***  CHECK ITOTAL
      IF ((ITOTAL.LT.0).OR.(ITOTAL.GT.1)) THEN
         PRINT *,'ITOTAL IN THE WRONG RANGE; ITOTAL = ',ITOTAL
         PRINT *,'ITOTAL CHANGED TO 1'
         ITOTAL = 1
      END IF

C***  CHECK ISCAPT AND SCAFAC
      IF ((ISCAPT.LT.0).OR.(ISCAPT.GT.1)) THEN
         PRINT *,'ISCAPT IN THE WRONG RANGE; ISCAPT = ',ISCAPT
         PRINT *,'ISCAPT CHANGED TO 0'
         ISCAPT = 0
      END IF
      IF (SCAFAC.LT.1D-2) THEN
         PRINT *, 'SCAFAC TOO SMALL'
         PRINT *, 'SCAFAC CHANGED TO 1.'
         SCAFAC = 1.D0
      END IF
      IF ((ITOTAL.EQ.1).AND.(ISCAPT.EQ.1)) THEN
         PRINT *,'WARNING: TRANSVERSE MASS IS NOT USED AS A SCALE'
      END IF

C***  CHECK PTMIN AND PTMAX
      IF (PTMIN.LT.0D0) THEN
         PTMIN = 0D0
         PRINT *,'PTMIN NEGATIVE'
         PRINT *,'PTMIN CHANGED TO 0.'
      END IF
      IF (PTMAX.LT.0D0) THEN
         PTMAX = 0D0
         PRINT *,'PTMAX NEGATIVE'
         PRINT *,'PTMAX CHANGED TO 0.'
      END IF
      IF (PTMAX.LT.PTMIN) THEN
         PTMAX = PTMIN
         PRINT *,'PTMAX SMALLER THAN PTMIN'
         PRINT *,'PTMAX CHANGED TO PTMIN'
      END IF
      IF ((PTMAX.EQ.PTMIN).AND.(PTMAX.EQ.0D0)) THEN
         PTMAX = 0.1D0
         PTMIN = 0.1D0
         PRINT *,'PTMAX = PTMIN CHANGED TO 0.1'
      END IF

C***  CHECK YMIN AND YMAX
      IF (YMIN.LT.0D0) THEN
         YMIN = 0D0
         PRINT *,'YMIN NEGATIVE'
         PRINT *,'YMIN CHANGED TO 0.'
      END IF
      IF (YMAX.LT.0D0) THEN
         YMAX = 0D0
         PRINT *,'YMAX NEGATIVE'
         PRINT *,'YMAX CHANGED TO 0.'
      END IF
      IF (YMAX.LT.YMIN) THEN
         YMAX = YMIN
         PRINT *,'YMAX SMALLER THAN YMIN'
         PRINT *,'YMAX CHANGED TO YMIN'
      END IF

C***  CHECK ICHACONJ
      IF ((ICHARCONJ.LT.-1).OR.(ICHARCONJ.GT.1)) THEN
         PRINT *,'ICHARCONJ IN THE WRONG RANGE; ICHARCONJ = ',ICHARCONJ
         PRINT *,'ICHARCONJ CHANGED TO 0'
         ICHARCONJ = 0
      END IF

C***  CHECK IONLYLO
      IF ((IONLYLO.LT.0).OR.(IONLYLO.GT.1)) THEN
         PRINT *,'IONLYLO IN THE WRONG RANGE'
         PRINT *,'IONLYLO CHANGED TO 0'
         IONLYLO = 0
      END IF

C***  CHECK ILO, INLO
      IF (ILO.LT.10) THEN
         ILO = 10
         PRINT *,'ILO CHANGED TO 10'
      END IF
      IF (INLO.LT.10) THEN
         INLO = 10
         PRINT *,'INLO CHANGED TO 10'
      END IF

      RETURN
      END

C ======================================================================

      SUBROUTINE INTEGST(RESLO,ERRLO,RESNLO,ERRNLO)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      integer  ivegas(1:4)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST5/ILO,INLO,IONLYLO

      EXTERNAL SIGLOST, SIGNLOST

      CALL CHECKINPUT_ST

      S = ENERGY**2

      CALL INILO(ivegas)
      idim = 3
      acc = 1.d-4
      call INTEG(SIGLOST,idim,ivegas,acc,reslo,errlo)

      IF (IONLYLO.EQ.0) THEN
         CALL ININLO(ivegas)
         idim = 5
         call INTEG(SIGNLOST,idim,ivegas,acc,resnlo,errnlo)
      ELSE
         RESNLO = 0D0
         ERRNLO = 0D0
      END IF

      RETURN
      END

C ======================================================================

C***  THE CROSS-SECTION IN LO

      REAL*8 FUNCTION SIGLOST(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION VAR(3)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/FLAVOR/IFLAVOR,ITOTAL
      COMMON/CHARCONJ/ICHARCONJ
      CALL DVEGLO_ST(VAR)
      IF (ITOTAL.EQ.0) THEN
         CALL DPTYST
         CALL DEFSCAST
         CALL DDILOST
         SIGLOST1 = ESTGGLO()
         SIGLOST2 = ESTQBLO()

C***     SUM OVER POSITIVE AND NEGATIVE RAPIDITY
         IF ((ICOLL.EQ.0).AND.(ICHARCONJ.NE.0)) THEN
            YVAL = -YVAL
            CALL DDILOST
            SIGLOST1 = 2 * SIGLOST1
            SIGLOST2 = SIGLOST2 + ESTQBLO()
         ELSE
            SIGLOST1 = 2 * SIGLOST1
            SIGLOST2 = 2 * SIGLOST2
         END IF
      ELSE
         CALL DEFSCAST
         CALL DTOLOST
         SIGLOST1 = FSTGGLO()
         SIGLOST2 = FSTQBLO()
      END IF
      SIGLOST =  SUMLOST(SIGLOST1,SIGLOST2,IFLAVOR)
      RETURN
      END

C ======================================================================

C***  THE CROSS-SECTION IN NLO

      REAL*8 FUNCTION SIGNLOST(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION VAR(4)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/FLAVOR/IFLAVOR,ITOTAL
      COMMON/CHARCONJ/ICHARCONJ

      CALL DVEGNLO_ST(VAR)
      CALL DMINLST
      IF (ITOTAL.EQ.0) THEN
         CALL DPTYST
         CALL DEFSCAST
         CALL DDINLST
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.1).OR.(IFLAVOR.EQ.4))
     +        SIGNLOST1 = ESTGGNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.2).OR.(IFLAVOR.EQ.4))
     +        SIGNLOST2 = ESTQBNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5))
     +        SIGNLOST3 = ESTGBNLO()

C***     SUM OVER POSITIVE AND NEGATIVE RAPIDITY
         IF ((ICOLL.EQ.0).AND.(ICHARCONJ.NE.0)) THEN
            SIGNLOST1 = 2 * SIGNLOST1
            YVAL = -YVAL
            CALL DDINLST
            IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.2).OR.(IFLAVOR.EQ.4))
     +           SIGNLOST2 = SIGNLOST2 + ESTQBNLO()
            IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5))
     +           SIGNLOST3 = SIGNLOST3 + ESTGBNLO()
         ELSE
            SIGNLOST1 = 2 * SIGNLOST1
            SIGNLOST2 = 2 * SIGNLOST2
            SIGNLOST3 = 2 * SIGNLOST3
         END IF
      ELSE
         CALL DEFSCAST
         CALL DTONLST
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.1).OR.(IFLAVOR.EQ.4))
     +        SIGNLOST1 = FSTGGNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.2).OR.(IFLAVOR.EQ.4))
     +        SIGNLOST2 = FSTQBNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5))
     +        SIGNLOST3 = FSTGBNLO()
      END IF
      SIGNLOST =  SUMNLOST(SIGNLOST1,SIGNLOST2,SIGNLOST3,IFLAVOR)
      RETURN
      END

C ======================================================================

      SUBROUTINE DVEGLO_ST(VAR)
C***  FILLS UP THE COMMON BLOCK VARVEGAS WITH 3 VARIABLES
      IMPLICIT REAL*8 (A-H,M-Z)
      DIMENSION VAR(3)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      Y1 = VAR(1)
      Y2 = VAR(2)
      Y3 = VAR(3)
      RETURN
      END

C ======================================================================

      SUBROUTINE DVEGNLO_ST(VAR)
C***  FILLS UP THE COMMON BLOCK VARVEGAS WITH 5 VARIABLES
      IMPLICIT REAL*8 (A-H,M-Z)
      DIMENSION VAR(4)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      Y1 = VAR(1)
      Y2 = VAR(2)
      Y3 = VAR(3)
      Y4 = VAR(4)
      RETURN
      END

C ======================================================================

      SUBROUTINE SWAP_ST(X,Y)
C***  INTERCHANGES X AND Y
      IMPLICIT REAL*8 (A-H,M-Z)
      Z = X
      X = Y
      Y = Z
      RETURN
      END

C ======================================================================

      SUBROUTINE AVER_ST(X,Y)
C***  AVERAGES X AND Y
      IMPLICIT REAL*8 (A-H,M-Z)
      X = ( X + Y )/2D0
      Y = X
      RETURN
      END

C ======================================================================

      SUBROUTINE DEFSCAST()
C***  DEFINES THE RENORMALIZATION AND FACTORIZATION SCALE
C***  FILLS UP PART OF THE COMMON BLOCKS CONST6 AND CONST2
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/FLAVOR/IFLAVOR,ITOTAL

      IF ((ITOTAL.EQ.0).AND.(ISCAPT.EQ.1)) THEN
         SCALE = MTRANS * SCAFAC
      ELSE
         SCALE = MST1 * SCAFAC
      END IF

C***  THE SCALE SHOULD BE BELOW 1 TEV, BECAUSE SOME PARTON DENSITIES
C***  ARE RESTRICTED TO THIS REGION

      IF (SCALE.GE.1000D0) SCALE = 1000D0

      IF (SCALE.LE.5D0) SCALE = 5D0

      ALPHAS = ALPS(SCALE)

      RETURN
      END

C ======================================================================

      SUBROUTINE DPTYST()
C***  DEFINES PT AND Y (INCLUDING CUTS)
C***  FILLS UP THE COMMON BLOCK PTYST
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CUT1/PTMINC,PTMAXC
      COMMON/CUT2/YMINC,YMAXC
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      PT2MAX = S/4.D0 -MST1**2
      PTMAX = SQRT( DMIN1(PT2MAX,PTMAXC**2) )
      PTMIN = PTMINC
      IF (PTMAX.GT.PTMIN) THEN
         PT = (PTMAX-PTMIN) * Y1 +PTMIN
         PREPT = 2.D0 * PT * (PTMAX -PTMIN)
      ELSE IF (PTMAX.EQ.PTMIN) THEN
         PT = PTMIN
         PREPT = 2.D0 * PT
      ELSE IF (PTMAX.LT.PTMIN) THEN
         PT = 0D0
         PREPT = 0D0
      END IF
      MTRANS = SQRT(PT**2 +MST1**2)
      YMAXK = +DACOSH(SQRT(S)/2.D0/MTRANS)
      YMINK = -DACOSH(SQRT(S)/2.D0/MTRANS)
      YMAX = DMIN1( YMAXK, YMAXC)
      YMIN = DMAX1( YMINK, YMINC, 0D0)

      IF (YMAX.GT.YMIN) THEN
         YVAL = (YMAX -YMIN) * Y2 + YMIN
         PREY = YMAX - YMIN
      ELSE IF (YMAX.EQ.YMIN) THEN
         YVAL = YMIN
         PREY = 1.D0
      ELSE IF (YMAX.LT.YMIN) THEN
         YVAL = 0D0
         PREY = 0D0
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE DDILOST()
C***  DEFINES THE KINEMATICS FOR DIFFERENTIAL CROSS-SECTIONS IN LO
C***  FILLS UP THE COMMON BLOCK DILOST AND PDLOST
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/DILOST/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      T1 = -SQRT(S)*MTRANS*DEXP(+YVAL)
      U1 = -SQRT(S)*MTRANS*DEXP(-YVAL)
      X1MIN = -T1/(S+U1)
      X1 = X1MIN ** Y3
      X2B= -X1*U1/(X1*S +T1)
      SHATB  = X1 *X2B*S
      T1HATB = X2B*T1
      U1HATB = X1 *U1
      PREFB = S * PREPT * PREY *LOG(X1MIN)*(-X1*X2B)/(-X1*U1)
      CALL PDFLOST(X1,X2B)
      RETURN
      END

C ======================================================================

      SUBROUTINE DTOLOST()
C***  DEFINES THE KINEMATICS FOR TOTAL CROSS-SECTIONS IN LO
C***  FILLS UP THE COMMON BLOCK TOLOST AND PDLOST
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/TOLOST/SHAT,T1HAT,PREF,X1,X2
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      X1 = (4*MST1**2/S) ** Y1
      X2 = (4*MST1**2/S/X1) ** Y2
      SHAT = X1*X2*S
      BETA = SQRT(1D0 -4*MST1**2/SHAT)
      T1HAT = -SHAT/2D0 +(Y3 -0.5D0)*SHAT*BETA
      PREF = SHAT*BETA *LOG(4*MST1**2/S)*LOG(4*MST1**2/S/X1)
      CALL PDFLOST(X1,X2)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESTGGLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN LO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/DILOST/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      ESTGGLO =
     +     PDGGB *PREFB * DSTGGB(ALPHAS,SHATB,T1HATB,MST1,MG)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESTQBLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN LO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      INTEGER ICHARCONJ
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/DILOST/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/CHARCONJ/ICHARCONJ

C***  DISTINGUISH BETWEEN STOPS AND ANTISTOPS IN THE FINAL STATE
C***  ICHARCONJ =  1  <-- STOPS
C***  ICHARCONJ = -1  <-- ANTISTOPS
C***  ICHARCONJ =  0  <-- AVERAGE OF STOPS AND ANTISTOPS

      IF (ICHARCONJ.EQ.1) THEN
         CALL SWAP_ST(PDQBB, PDBQB)
      ELSE IF (ICHARCONJ.EQ.0) THEN
         CALL AVER_ST(PDQBB, PDBQB)
      END IF

      ESTQBLO =
     +     +  PDQBB * PREFB*DSTQBB(ALPHAS,SHATB,T1HATB,MST1,MG)
     +     +  PDBQB * PREFB*DSTQBB(ALPHAS,SHATB,U1HATB,MST1,MG)

      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSTGGLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN LO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/TOLOST/SHAT,T1HAT,PREF,X1,X2
      FSTGGLO =
     +     PDGGB * PREF*DSTGGB(ALPHAS,SHAT,T1HAT,MST1,MG)
      RETURN
      END

C ======================================================================
      REAL*8 FUNCTION FSTQBLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN LO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/TOLOST/SHAT,T1HAT,PREF,X1,X2
      FSTQBLO =
     +     + (PDQBB +PDBQB ) * PREF*DSTQBB(ALPHAS,SHAT,T1HAT,MST1,MG)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION SUMLOST(SIGLOST1,SIGLOST2,IFLAVOR)
C***  SUMS UP DIFFERENT FLAVOR STATES FOR CROSS-SECTION IN LO
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I)
      IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.4))
     +                  SUMLOST = SIGLOST1 + SIGLOST2
      IF (IFLAVOR.EQ.1) SUMLOST = SIGLOST1
      IF (IFLAVOR.EQ.2) SUMLOST = SIGLOST2
      IF ((IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5)) SUMLOST = 0D0
      RETURN
      END

C ======================================================================

      SUBROUTINE DMINLST()
C***  DEFINES SOME MISCELLANEOUS VARIABLES FOR CROSS-SECTION IN NLO
C***  FILLS UP THE COMMON BLOCK MINLST
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/MINLST/DEL,SCA
      DEL = 1.D-4 * MST1**2
      SCA = SCALE**2/MST1**2
      RETURN
      END

C ======================================================================

      SUBROUTINE DDINLST()
C***  DEFINES THE KINEMATICS FOR DIFFERENTIAL CROSS-SECTIONS IN NLO
C***  FILLS UP THE COMMON BLOCK DINLST1-2 AND PDLOST, PDNLST1
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/MINLST/DEL,SCA
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/PDNLST1/PDGG,PDQB,PDBQ,PDQG, PDGQ,PDBG, PDGB

      T1 = -SQRT(S)*MTRANS*DEXP(+YVAL)
      U1 = -SQRT(S)*MTRANS*DEXP(-YVAL)
      X1MIN = -T1/(S+U1)
      X1 = X1MIN ** Y3
      S4MAX = X1 * (S+U1) +T1
      S4 = (S4MAX -DEL)* Y4 +DEL
      X2 = (S4 -X1*U1)/(X1*S +T1)
      X2B= -X1*U1/(X1*S +T1)
      SHAT   = X1 *X2 *S
      SHATB  = X1 *X2B*S
      T1HAT  = X2 *T1
      T1HATB = X2B*T1
      U1HAT  = X1 *U1
      U1HATB = X1 *U1
      PREF  = S * PREPT * PREY *LOG(X1MIN)*(-X1*X2 )
     +     *(S4MAX -DEL)/(S4-X1*U1)
      PREFB = S * PREPT * PREY *LOG(X1MIN)*(-X1*X2B)/(-X1*U1)
      CALL PDFNLOST1(X1,X2,X2B)
      RETURN
      END

C ======================================================================

      SUBROUTINE DTONLST()
C***  DEFINES THE KINEMATICS FOR TOTAL CROSS-SECTIONS IN NLO
C***  FILLS UP THE COMMON BLOCKS DINLST1-2 AND PDNLSTT
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/MINLST/DEL,SCA
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/PDNLSTT/PDGGT, PDQBT, PDGQT
      X1 = (4*MST1**2/S) ** Y1
      X2 = (4*MST1**2/S/X1) ** Y2
      SHAT = X1*X2*S
      DEL = 1.D-4 * MST1**2 * (1D0 -4*MST1**2/SHAT)
      BETA = SQRT(1D0 -4*MST1**2/SHAT)
      T1HAT =
     + -(SHAT -DEL)/2D0 +(Y3 -0.5D0)*SQRT((SHAT-DEL)**2-4*MST1**2*SHAT)
      S4MAX = SHAT +T1HAT +SHAT*MST1**2/T1HAT
      S4 = (S4MAX -DEL)*Y4 +DEL
      PREF  = SHAT*BETA *S4MAX*LOG(4*MST1**2/S)*LOG(4*MST1**2/S/X1)
      PREFB = SHAT*BETA *      LOG(4*MST1**2/S)*LOG(4*MST1**2/S/X1)
      CALL PDFNLOSTT(X1,X2)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESTGGNLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/MINLST/DEL,SCA
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/PDNLST1/PDGG,PDQB,PDBQ,PDQG, PDGQ,PDBG, PDGB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      ESTGGNLO =
     +     +PDGGB * PREFB * DSTGGB(ALPHAS,SHATB,T1HATB,MST1,MG)
     +     +PDGGB * PREFB *
     +     DSTGGV(ALPHAS,SHATB,T1HATB,MST1,MG,MT,MS,MST2,sin2t)
     +     +PDGGB * PREFB * DSTGG1(ALPHAS,SHATB,T1HATB,MST1,MG,SCA)
     +     +PDGGB * PREFB *
     +     DSTGGD(ALPHAS,SHATB,T1HATB,S4,MST1,MG,DEL,S4MAX)
     +     +PDGGB * PREFB *
     +     DSTGG2(ALPHAS,SHATB,T1HATB,S4,MST1,MG,DEL,S4MAX,SCA)
     +     +PDGG  * PREF  * DSTGGH(ALPHAS,SHAT,T1HAT,S4,MST1,MG)
     +     +PDGG  * PREF  * DSTGG3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESTQBNLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      INTEGER ICHARCONJ
      COMMON/MINLST/DEL,SCA
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/PDNLST1/PDGG,PDQB,PDBQ,PDQG, PDGQ,PDBG, PDGB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/CHARCONJ/ICHARCONJ

C***  DISTINGUISH BETWEEN STOPS AND ANTISTOPS IN THE FINAL STATE
C***  ICHARCONJ =  1  <-- STOPS
C***  ICHARCONJ = -1  <-- ANTISTOPS
C***  ICHARCONJ =  0  <-- AVERAGE OF STOPS AND ANTISTOPS

      IF (ICHARCONJ.EQ.1) THEN
         CALL SWAP_ST(PDQBB, PDBQB)
         CALL SWAP_ST(PDQB , PDBQ )
      ELSE IF (ICHARCONJ.EQ.0) THEN
         CALL AVER_ST(PDQBB, PDBQB)
         CALL AVER_ST(PDQB , PDBQ )
      END IF

      ESTQBNLO =
     +        +PDQBB * PREFB * DSTQBB(ALPHAS,SHATB,T1HATB,MST1,MG)
     +        +PDQBB * PREFB *
     +        DSTQBV(ALPHAS,SHATB,T1HATB,MST1,MG,MT,MS,MST2,SIN2T)
     +        +PDQBB * PREFB * DSTQB1(ALPHAS,SHATB,T1HATB,MST1,MG,SCA)
     +        +PDQBB * PREFB *
     +        DSTQBD(ALPHAS,SHATB,T1HATB,S4,MST1,MG,DEL,S4MAX)
     +        +PDQBB * PREFB *
     +        DSTQB2(ALPHAS,SHATB,T1HATB,S4,MST1,MG,DEL,S4MAX,SCA)
     +        +PDQB  * PREF  *
     +        DSTQBH(ALPHAS,SHAT,T1HAT,S4,MST1,MG)
     +        +PDQB  * PREF  *
     +        DSTQB3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA)

      ESTQBNLO = ESTQBNLO
     +        +PDBQB * PREFB * DSTQBB(ALPHAS,SHATB,U1HATB,MST1,MG)
     +        +PDBQB * PREFB *
     +        DSTQBV(ALPHAS,SHATB,U1HATB,MST1,MG,MT,MS,MST2,SIN2T)
     +        +PDBQB * PREFB * DSTQB1(ALPHAS,SHATB,U1HATB,MST1,MG,SCA)
     +        +PDBQB * PREFB *
     +        DSTQBD(ALPHAS,SHATB,U1HATB,S4,MST1,MG,DEL,S4MAX)
     +        +PDBQB * PREFB *
     +        DSTQB2(ALPHAS,SHATB,U1HATB,S4,MST1,MG,DEL,S4MAX,SCA)
     +        +PDBQ  * PREF  *
     +        DSTQBH(ALPHAS,SHAT,U1HAT,S4,MST1,MG)
     +        +PDBQ  * PREF  *
     +        DSTQB3(ALPHAS,SHAT,U1HAT,S4,MST1,MG,SCA)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESTGBNLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR G QB AND G Q
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/MINLST/DEL,SCA
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/PDNLST1/PDGG,PDQB,PDBQ,PDQG, PDGQ,PDBG, PDGB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/CHARCONJ/ICHARCONJ

C***  DISTINGUISH BETWEEN STOPS AND ANTISTOPS IN THE FINAL STATE
C***  ICHARCONJ =  1  <-- STOPS
C***  ICHARCONJ = -1  <-- ANTISTOPS
C***  ICHARCONJ =  0  <-- AVERAGE OF STOPS AND ANTISTOPS

      IF (ICHARCONJ.EQ.1) THEN
         CALL SWAP_ST(PDGB, PDGQ)
         CALL SWAP_ST(PDBG, PDQG)
      ELSE IF (ICHARCONJ.EQ.0) THEN
         CALL AVER_ST(PDGB, PDGQ)
         CALL AVER_ST(PDBG, PDQG)
      END IF

      ESTGBNLO =
     +        +PDGB  * PREF  *
     +        DSTGBH(ALPHAS,SHAT,T1HAT,S4,MST1,MG)
     +        +PDGB  * PREF  *
     +        DSTGB3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA)
     +        +PDBG  * PREF  *
     +        DSTGBH(ALPHAS,SHAT,U1HAT,S4,MST1,MG)
     +        +PDBG  * PREF  *
     +        DSTGB3(ALPHAS,SHAT,U1HAT,S4,MST1,MG,SCA)
     +
     +        +PDQG  * PREF  *
     +        DSTQGH(ALPHAS,SHAT,T1HAT,S4,MST1,MG)
     +        +PDQG  * PREF  *
     +        DSTQG3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA)
     +        +PDGQ  * PREF  *
     +        DSTQGH(ALPHAS,SHAT,U1HAT,S4,MST1,MG)
     +        +PDGQ  * PREF  *
     +        DSTQG3(ALPHAS,SHAT,U1HAT,S4,MST1,MG,SCA)

      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSTGGNLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN NLO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/PDNLSTT/PDGGT, PDQBT, PDGQT
      COMMON/MINLST/DEL,SCA
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      FSTGGNLO =
     +        +PDGGT * PREFB * DSTGGB(ALPHAS,SHAT,T1HAT,MST1,MG)
     +        +PDGGT * PREFB *
     +        DSTGGV(ALPHAS,SHAT,T1HAT,MST1,MG,MT,MS,MST2,SIN2T)
     +        +PDGGT * PREFB * DSTGG1(ALPHAS,SHAT,T1HAT,MST1,MG,SCA)
     +        +PDGGT * PREFB *
     +        DSTGGD(ALPHAS,SHAT,T1HAT,S4,MST1,MG,DEL,S4MAX)
     +        +PDGGT * PREFB *
     +        DSTGG2(ALPHAS,SHAT,T1HAT,S4,MST1,MG,DEL,S4MAX,SCA)
     +        +PDGGT * PREF  * DSTGGH(ALPHAS,SHAT,T1HAT,S4,MST1,MG)
     +        +PDGGT * PREF  * DSTGG3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSTQBNLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN NLO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDNLSTT/PDGGT, PDQBT, PDGQT
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/MINLST/DEL,SCA
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      FSTQBNLO =
     +        +PDQBT * PREFB * DSTQBB(ALPHAS,SHAT,T1HAT,MST1,MG)
     +        +PDQBT * PREFB *
     +        DSTQBV(ALPHAS,SHAT,T1HAT,MST1,MG,MT,MS,MST2,SIN2T)
     +        +PDQBT * PREFB * DSTQB1(ALPHAS,SHAT,T1HAT,MST1,MG,SCA)
     +        +PDQBT * PREFB *
     +        DSTQBD(ALPHAS,SHAT,T1HAT,S4,MST1,MG,DEL,S4MAX)
     +        +PDQBT * PREFB *
     +        DSTQB2(ALPHAS,SHAT,T1HAT,S4,MST1,MG,DEL,S4MAX,SCA)
     +        +PDQBT * PREF  *
     +        DSTQBH(ALPHAS,SHAT,T1HAT,S4,MST1,MG)
     +        +PDQBT * PREF  *
     +        DSTQB3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSTGBNLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN NLO FOR G QB AND G Q
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/PDNLSTT/PDGGT, PDQBT, PDGQT
      COMMON/MINLST/DEL,SCA
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      FSTGBNLO =
     +        +PDGQT * PREF  *
     +        DSTGBH(ALPHAS,SHAT,T1HAT,S4,MST1,MG)
     +        +PDGQT * PREF  *
     +        DSTGB3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION SUMNLOST(SIGNLOST1,SIGNLOST2,SIGNLOST3,IFLAVOR)
C***  SUMS UP DIFFERENT FLAVOR STATES FOR CROSS-SECTION IN NLO
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I)
      IF (IFLAVOR.EQ.0) SUMNLOST = SIGNLOST1 +SIGNLOST2 +SIGNLOST3
      IF (IFLAVOR.EQ.1) SUMNLOST = SIGNLOST1
      IF (IFLAVOR.EQ.2) SUMNLOST = SIGNLOST2
      IF (IFLAVOR.EQ.3) SUMNLOST = SIGNLOST3
      IF (IFLAVOR.EQ.4) SUMNLOST = SIGNLOST1 +SIGNLOST2
      IF (IFLAVOR.EQ.5) SUMNLOST = SIGNLOST3
      IF (.NOT.(ABS(SUMNLOST).LT.1.D40)) THEN
         SUMNLOST = 0.D0
         PRINT *,'NAN'
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE PDFLOST(X1,X2)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB

      CALL PARTONDF(X1,SCALE,XPDF)
      CALL PARTONDF(X2,SCALE,YPDF)

      PDGGB = XPDF(0) * YPDF(0)
      PDQBB = 0D0
      PDBQB = 0D0
      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF( I)
            PDBQB = PDBQB + XPDF(-I) * YPDF(-I)
 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF(-I)
            PDBQB = PDBQB + XPDF(-I) * YPDF( I)
 1010    CONTINUE
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE PDFNLOSTT(X1,X2)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDNLSTT/PDGGT, PDQBT, PDGQT

      CALL PARTONDF(X1,SCALE,XPDF)
      CALL PARTONDF(X2,SCALE,YPDF)

      PDGGT = XPDF(0) * YPDF(0)
      PDQBT = 0D0
      PDGQT = 0D0

      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDQBT = PDQBT + XPDF( I) * YPDF( I) + XPDF(-I) * YPDF(-I)
            PDGQT = PDGQT + XPDF(0)  * YPDF(-I) + XPDF( I) * YPDF( 0)
     +                    + XPDF(0)  * YPDF(I)  + XPDF(-I) * YPDF( 0)
 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDQBT = PDQBT + XPDF( I) * YPDF(-I) + XPDF(-I) * YPDF(I)
            PDGQT = PDGQT + XPDF(0)  * YPDF(-I) + XPDF( I) * YPDF(0)
     +                    + XPDF(0)  * YPDF(I)  + XPDF(-I) * YPDF(0)
 1010    CONTINUE
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE PDFNLOST1(X1,X2,X2B)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      DIMENSION ZPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/PDNLST1/PDGG,PDQB,PDBQ,PDQG, PDGQ,PDBG, PDGB

      CALL PARTONDF(X1 ,SCALE,XPDF)
      CALL PARTONDF(X2B,SCALE,YPDF)
      CALL PARTONDF(X2 ,SCALE,ZPDF)

      PDGGB = XPDF(0) * YPDF(0)
      PDGG  = XPDF(0) * ZPDF(0)
      PDQBB = 0D0
      PDQB  = 0D0
      PDBQB = 0D0
      PDBQ  = 0D0
      PDQG  = 0D0
      PDGQ  = 0D0
      PDBG  = 0D0
      PDGB  = 0D0

      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF( I)
            PDQB  = PDQB  + XPDF( I) * ZPDF( I)
            PDBQB = PDBQB + XPDF(-I) * YPDF(-I)
            PDBQ  = PDBQ  + XPDF(-I) * ZPDF(-I)
            PDQG  = PDQG  + XPDF(+I) * ZPDF( 0)
            PDGQ  = PDGQ  + XPDF( 0) * ZPDF(-I)
            PDBG  = PDBG  + XPDF(-I) * ZPDF( 0)
            PDGB  = PDGB  + XPDF( 0) * ZPDF( I)

 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF(-I)
            PDQB  = PDQB  + XPDF( I) * ZPDF(-I)
            PDBQB = PDBQB + XPDF(-I) * YPDF( I)
            PDBQ  = PDBQ  + XPDF(-I) * ZPDF( I)
            PDQG  = PDQG  + XPDF( I) * ZPDF( 0)
            PDGQ  = PDGQ  + XPDF( 0) * ZPDF( I)
            PDBG  = PDBG  + XPDF(-I) * ZPDF( 0)
            PDGB  = PDGB  + XPDF( 0) * ZPDF(-I)
 1010    CONTINUE
      END IF
      RETURN
      END

C ======================================================================

max_stars_repo_name: sliem/docker-hep

max_issues_repo_name: sliem/docker-hep

max_forks_repo_name: sliem/docker-hep


    !-------------------------------------------------------------------------------------------------------------
    !
    !> \file    CompThermoData.f90
    !> \brief   Compute thermodynamic data
    !> \author  M.H.A. Piro
    !> \date    Apr. 24, 2012
    !> \sa      Thermochimica.f90
    !> \sa      CompGibbsMagnetic.f90
    !
    !
    ! Revisions:
    ! ==========
    !
    !   Date            Programmer          Description of change
    !   ----            ----------          ---------------------
    !   10/11/2011      M.H.A. Piro         Original code
    !   05/01/2012      M.H.A. Piro         Corrected storage of mixing component indices when solution species
    !                                        are no longer considered.  This affected the variables iSpeciesPass
    !                                        and iRegularParam.
    !   01/31/2012      M.H.A. Piro         The standard molar Gibbs energy coefficient of gaseous species is
    !                                        adjusted by an arbitrary quantity (see comments below).
    !   01/17/2013      M.H.A. Piro         Additional dummy species are included to the system for each electron
    !                                        representing a system component.
    !   02/11/2013      M.H.A. Piro         Mixing parameters for SUBL phases have been incorporated.  These are
    !                                        different than other phases because the constituent indices on each
    !                                        sublattice are stored, rather than the species indices.
    !   07/16/2013      M.H.A. Piro         Fix bug in constructing the iRegularParam vector when dealing with SUBL
    !                                        phases when the mixing constituents are not on the first sublattice.
    !   03/19/2018      M.H.A. Piro         Added capability to handle SUBG phases. This includes a special case
    !                                        to define the stoichiometry of pairs, rather than pure species. Otherwise
    !                                        they would be undefined.
    !
    !
    ! Purpose:
    ! ========
    !
    !> \details The purpose of this subroutine is to compute thermodynamic data for all substances in the
    !! system.  The specific variables that are computed include:
    !!   - the standard molar Gibbs energy for each pure substance using the specified temperature and pressure,
    !!   - the excess molar Gibbs energy of mixing for each sub-system,
    !!   - the total number of atoms per formula mass for each compound,
    !!   - the atomic fraction of each element in a particular compound, and
    !!   - construction of the Hessian matrix, which will be used later by the GEMNewton subroutine.
    !!
    !! The coefficients for the standard molar Gibbs energy equations for pure species originate from a ChemSage
    !! data-file that was parsed from the ParseCSDataFile program.  The format for the coefficients follow:
    !!
    !! \f$ g_i^{\circ} = A + BT + CTln(T) + DT^2 + ET^3 + F/T + (GT^U + HT^V + IT^W + JT^X + KT^Y + LT^Z) \f$
    !!
    !! Note that the terms in paranteses in the above equation are additional terms.  Some of the
    !! exponents used in the additional terms may be 99.  This corresponds to the natural logarithm.
    !!
    !! For some very strange reason, the B coefficient in the above equation is modified by ~0.10945 J/mol for
    !! only gaseous species when FactSage generates a ChemSage data-file.  Another very peculiar observation is
    !! that if a database is constructed in FactSage using the new ChemSage data-file, FactSage will automatically
    !! remove this quantity from the B coefficient.  The B coefficient for gaseous species is corrected in this
    !! subroutine.
    !
    !
    ! Pertinent variables:
    ! ====================
    !
    ! nSpecies              The number of species in the system.
    ! nSpeciesCS            The number of species in the original ChemSage data-file.  Note: this may be
    !                        higher than what may be considered in this system.
    ! nSpeciesPhase         An integer vector representing the number of species in each solution phase.
    ! nSpeciesPhaseCS       An integer vector representing the number of species in each solution phase
    !                        from the ChemSage data-file (note: this may be larger than nSpeciesPhase).
    ! dChemicalPotential    A double real vector representing the chemical potential of each species.
    !                        To be precise, this is defined as the difference between the standard molar
    !                        Gibbs energy and the chemical potential defined by the element potentials
    !                        (represented in dimensionless units and per formula mass).
    ! dExcessGibbsParam     A double real vector representing the excess Gibbs energy of mixing
    !                        parameters.
    ! dAtomFractionSpecies  A double real matrix representing tha atom fraction of each element in each
    !                        species.
    ! dJacobianLong         A double real matrix that is used to construct the Jacobian matrix in
    !                        GEMNewton.f90. This matrix is constructed once and never modified throughout
    !                        a calculation.
    ! iParticlesPerMole     The number of particles per mole of constituent.
    ! iSpeciesAtoms         An integer matrix representing the number of atoms of a particular element
    !                        for a particular species.
    ! iSpeciesAtomsCS       An integer matrix representing the number of atoms of a particular element
    !                        for a particular species from the database (note: this may be larger than
    !                        iSpecesAtoms).
    ! iSpeciesTotalAtoms    The total number of atoms per formula mass of a species.
    ! iRegularParam         An integer vector representing important mixing terms for a regular solution
    !                        model.
    ! iRegularParamCS       An integer vector representing important mixing terms for a regular solution
    !                        model from the ChemSage data-file (Note: this may be larger than
    !                        iRegularParam).
    ! dGibbsCoeffSpeciesTemp  A double real array containing the coefficients of Gibbs energy equations
    !                        of species in the database.
    !
    !-------------------------------------------------------------------------------------------------------------

subroutine CompThermoData

    USE ModuleParseCS
    USE ModuleThermo
    USE ModuleThermoIO

    implicit none

    integer                            :: i, j, k, l, m, n, s, iCounterGibbsEqn, nCounter, l1, l2, nn, c
    integer                            :: ii, jj, kk, ll, ka, la, iax, iay, ibx, iby, ia2x2, ia2y2, ib2x2, ib2y2
    integer                            :: iSublPhaseIndex, iFirst, iLast, nRemove, nA2X2, nTempSublattice
    integer                            :: iMixStart, iMixLength, nMixSets
    integer, dimension(nElementsCS**2) :: iRemove
    real(8)                            :: dLogT, dLogP, dTemp, dQx, dQy, dZa, dZb, dZx, dZy, dCoax
    real(8)                            :: dStdEnergyTemp
    real(8), dimension(6)              :: dGibbsCoeff
    real(8), dimension(nSpeciesCS)     :: dChemicalPotentialTemp
    character(12), dimension(:),     allocatable :: cElementNameTemp
    real(8),       dimension(:),     allocatable :: dMolesElementTemp
    real(8),       dimension(:,:),   allocatable :: dAtomFractionSpeciesTemp, dStoichSpeciesTemp
    real(8),       dimension(:,:,:), allocatable :: dStoichPairsTemp

    ! Initialize variables:
    j                = 0
    iCounterGibbsEqn = 0
    nDummySpecies    = 0
    nCounter         = 0
    dTemp            = 1D0 / (dIdealConstant * dTemperature)
    iSUBIMixType     = 0
    nTempSublattice  = 0

    ! Compute Gibbs energy coefficients:
    dGibbsCoeff(1)   = 1D0                             ! A
    dGibbsCoeff(2)   = dTemperature                    ! B
    dGibbsCoeff(3)   = dTemperature*DLOG(dTemperature) ! C
    dGibbsCoeff(4)   = dTemperature**2                 ! D
    dGibbsCoeff(5)   = dTemperature**3                 ! E
    dGibbsCoeff(6)   = 1D0 / dTemperature              ! F
    dLogT            = DLOG(dTemperature)              ! ln(T)
    dLogP            = DLOG(dPressure)                 ! ln(P)

    if (.NOT. lHeatCapacityCurrent) dTemperatureForLimits = dTemperature

    ! Loop through all species in the system:
    LOOP_nPhasesCS: do n = 1, nSolnPhasesSysCS
        iFirst = nSpeciesPhaseCS(n-1) + 1
        iLast  = nSpeciesPhaseCS(n)
        l = MAXVAL(iSpeciesPass(iFirst:iLast))
        if (l > 0) then
            if ((cSolnPhaseTypeCS(n) == 'SUBL').OR.(cSolnPhaseTypeCS(n) == 'SUBLM').OR. &
                 (cSolnPhaseTypeCS(n) == 'SUBG').OR.(cSolnPhaseTypeCS(n) == 'SUBQ').OR. &
                 (cSolnPhaseTypeCS(n) == 'SUBI')) then
                nTempSublattice = nTempSublattice + 1
            end if
        end if

        if ((cSolnPhaseTypeCS(n) == 'SUBG') .OR. (cSolnPhaseTypeCS(n) == 'SUBQ')) then
            iSublPhaseIndex = iPhaseSublatticeCS(n)
            iFirst = nSpeciesPhaseCS(n - 1) + 1
            dChemicalPotentialTemp = 0D0
            jj = 0
            LOOP_SROPairs: do i = iFirst, iFirst - 1 + nPairsSROCS(iSublPhaseIndex,1)
                l = 0
                ! Loop through the Gibbs energy equations to figure out which one to use:
                do k = 1, nGibbsEqSpecies(i)
                    iCounterGibbsEqn = iCounterGibbsEqn + 1
                    if ((dTemperatureForLimits <= dGibbsCoeffSpeciesTemp(1,iCounterGibbsEqn)).AND.(l == 0)) then
                        l = k
                    end if
                end do

                if (l == 0) l = nGibbsEqSpecies(i)

                l = l + iCounterGibbsEqn - nGibbsEqSpecies(i)

                do k = 2, 7
                    dChemicalPotentialTemp(i) = dChemicalPotentialTemp(i) + dGibbsCoeffSpeciesTemp(k,l) * dGibbsCoeff(k-1)
                end do

                ! Compute additional standard molar Gibbs energy terms:
                if (dGibbsCoeffSpeciesTemp(9,l) .EQ. 99) then
                    dChemicalPotentialTemp(i) = dChemicalPotentialTemp(i) + dGibbsCoeffSpeciesTemp(8,l) * dLogT
                else
                    dChemicalPotentialTemp(i) = dChemicalPotentialTemp(i) + dGibbsCoeffSpeciesTemp(8,l) &
                        *dTemperature**dGibbsCoeffSpeciesTemp(9,l)
                end if

                if (dGibbsCoeffSpeciesTemp(11,l) .EQ. 99) then
                    dChemicalPotentialTemp(i) = dChemicalPotentialTemp(i) + dGibbsCoeffSpeciesTemp(10,l) * dLogT
                else
                    dChemicalPotentialTemp(i) = dChemicalPotentialTemp(i) + dGibbsCoeffSpeciesTemp(10,l) &
                        * dTemperature**dGibbsCoeffSpeciesTemp(11,l)
                end if

                if (dGibbsCoeffSpeciesTemp(13,l) .EQ. 99) then
                    dChemicalPotentialTemp(i) = dChemicalPotentialTemp(i) + dGibbsCoeffSpeciesTemp(12,l) * dLogT
                else
                    dChemicalPotentialTemp(i) = dChemicalPotentialTemp(i) + dGibbsCoeffSpeciesTemp(12,l) &
                        * dTemperature**dGibbsCoeffSpeciesTemp(13,l)
                end if

                ! Convert chemical potentials to dimensionless units:
                dChemicalPotentialTemp(i) = dChemicalPotentialTemp(i) * dTemp * DFLOAT(iParticlesPerMoleCS(i))

                do k = 1, nElemOrComp
                    if ((dStoichPairsCS(iSublPhaseIndex,i - iFirst + 1,k) > 0).AND.(iElementSystem(k) == 0)) then
                        cycle LOOP_SROPairs
                    end if
                end do

                ! Check if pair should be saved - only worry about zeta here, LOOP_nSUBGQCS will take care
                ! of only using the necessary reference energy terms.
                nPairsSRO(nTempSublattice,1) = nPairsSRO(nTempSublattice,1) + 1
                jj = jj + 1
                dZetaSpecies(nTempSublattice,jj) = dZetaSpeciesCS(iSublPhaseIndex,i - iFirst + 1)
                dConstituentCoefficients(nTempSublattice,jj,1:5) = &
                                    dConstituentCoefficientsCS(iSublPhaseIndex,i - iFirst + 1,1:5)
                cPairName(nTempSublattice,jj) = cPairNameCS(iSublPhaseIndex,i - iFirst + 1)
                m = 0
                do k = 1, nElemOrComp
                    if (iElementSystem(k) /= 0) then
                        m = m + 1
                        dStoichPairs(nTempSublattice,jj,m) = dStoichPairsCS(iSublPhaseIndex,i - iFirst + 1,k)
                    end if
                end do

                ! I have set these to be g_A2/X2 * Z^A_A2/X2 because it helps me to follow,
                ! this is not the most computationally efficient option.
                dCoax = dConstituentCoefficientsCS(iSublPhaseIndex,i - iFirst + 1,1)
                dChemicalPotentialTemp(i) = dChemicalPotentialTemp(i) * 2D0 / dCoax
            end do LOOP_SROPairs

            LOOP_nSUBGQCS: do i = nSpeciesPhaseCS(n - 1) + 1, nSpeciesPhaseCS(n)
                if (iSpeciesPass(i) == 0) cycle LOOP_nSUBGQCS

                j = j + 1   ! New species index

                cSpeciesName(j)            = cSpeciesNameCS(i)
                iPhase(j)                  = iPhaseCS(i)
                iParticlesPerMole(j)       = iParticlesPerMoleCS(i)
                dCoeffGibbsMagnetic(j,1:4) = dGibbsMagneticCS(i,1:4)

                m = 0
                do k = 1, nElemOrComp
                    if (iElementSystem(k) /= 0) then
                        m = m + 1
                        dStoichSpecies(j,m) = dStoichSpeciesCS(i,k)
                    end if
                end do

                dZa = dCoordinationNumberCS(iSublPhaseIndex,i - iFirst + 1,1)
                dZb = dCoordinationNumberCS(iSublPhaseIndex,i - iFirst + 1,2)
                dZx = dCoordinationNumberCS(iSublPhaseIndex,i - iFirst + 1,3)
                dZy = dCoordinationNumberCS(iSublPhaseIndex,i - iFirst + 1,4)

                ii = iPairIDCS(iSublPhaseIndex,i - iFirst + 1,1)
                jj = iPairIDCS(iSublPhaseIndex,i - iFirst + 1,2)
                kk = iPairIDCS(iSublPhaseIndex,i - iFirst + 1,3)
                ll = iPairIDCS(iSublPhaseIndex,i - iFirst + 1,4)
                ! Anion indices adjusted to start from 1
                ka = kk - nConstituentSublatticeCS(iSublPhaseIndex,1)
                la = ll - nConstituentSublatticeCS(iSublPhaseIndex,1)

                dQx = dSublatticeChargeCS(iSublPhaseIndex,2,ka)
                dQy = dSublatticeChargeCS(iSublPhaseIndex,2,la)

                ! Indices of reference energy equations for A/X, B/X, A/Y, B/Y
                iax = 0
                ibx = 0
                iay = 0
                iby = 0
                nA2X2 = nConstituentSublatticeCS(iSublPhaseIndex,1) * nConstituentSublatticeCS(iSublPhaseIndex,2)
                do k = 1, nA2X2
                    if   ((iConstituentSublatticeCS(iSublPhaseIndex,1,k) == ii) &
                    .AND. (iConstituentSublatticeCS(iSublPhaseIndex,2,k) == ka)) then
                        iax = k
                    end if
                    if   ((iConstituentSublatticeCS(iSublPhaseIndex,1,k) == jj) &
                    .AND. (iConstituentSublatticeCS(iSublPhaseIndex,2,k) == ka)) then
                        ibx = k
                    end if
                    if   ((iConstituentSublatticeCS(iSublPhaseIndex,1,k) == ii) &
                    .AND. (iConstituentSublatticeCS(iSublPhaseIndex,2,k) == la)) then
                        iay = k
                    end if
                    if   ((iConstituentSublatticeCS(iSublPhaseIndex,1,k) == jj) &
                    .AND. (iConstituentSublatticeCS(iSublPhaseIndex,2,k) == la)) then
                        iby = k
                    end if
                end do

                ia2x2 = ii + ((ka - 1) * (nConstituentSublatticeCS(iSublPhaseIndex,1) &
                                        * (nConstituentSublatticeCS(iSublPhaseIndex,1) + 1) / 2))
                ib2x2 = jj + ((ka - 1) * (nConstituentSublatticeCS(iSublPhaseIndex,1) &
                                        * (nConstituentSublatticeCS(iSublPhaseIndex,1) + 1) / 2))
                ia2y2 = ii + ((la - 1) * (nConstituentSublatticeCS(iSublPhaseIndex,1) &
                                        * (nConstituentSublatticeCS(iSublPhaseIndex,1) + 1) / 2))
                ib2y2 = jj + ((la - 1) * (nConstituentSublatticeCS(iSublPhaseIndex,1) &
                                        * (nConstituentSublatticeCS(iSublPhaseIndex,1) + 1) / 2))

                ! This is equation 17.38 in Pelton's book
                dChemicalPotential(j) = ((dQx * dChemicalPotentialTemp(iax + iFirst - 1) / (2D0 * dZa * dZx))  &
                                       + (dQx * dChemicalPotentialTemp(ibx + iFirst - 1) / (2D0 * dZb * dZx))  &
                                       + (dQy * dChemicalPotentialTemp(iay + iFirst - 1) / (2D0 * dZa * dZy))  &
                                       + (dQy * dChemicalPotentialTemp(iby + iFirst - 1) / (2D0 * dZb * dZy))) &
                                       / ((dQx/dZx) + (dQy/dZy))
            end do LOOP_nSUBGQCS
        else
            LOOP_nSpeciesCS: do i = nSpeciesPhaseCS(n - 1) + 1, nSpeciesPhaseCS(n)
                l1 = 0
                ! Loop through the Gibbs energy equations to figure out which one to use:
                do k = 1, nGibbsEqSpecies(i)
                    iCounterGibbsEqn = iCounterGibbsEqn + 1
                    if ((dTemperatureForLimits <= dGibbsCoeffSpeciesTemp(1,iCounterGibbsEqn)).AND.(l1 == 0)) then
                        l1 = k
                    end if
                end do

                ! This species will not be considered part of the system if it didn't pass.
                if (iSpeciesPass(i) == 0) cycle LOOP_nSpeciesCS

                l2 = l1
                if (l2 == 0) l2 = nGibbsEqSpecies(i)

                l2 = l2 + iCounterGibbsEqn - nGibbsEqSpecies(i)
                j = j + 1   ! New species index

                cSpeciesName(j)            = cSpeciesNameCS(i)
                iPhase(j)                  = iPhaseCS(i)
                iParticlesPerMole(j)       = iParticlesPerMoleCS(i)
                dCoeffGibbsMagnetic(j,1:4) = dGibbsMagneticCS(i,1:4)
                dQKTOParams(j,1:2)         = dQKTOParamsCS(i,1:2)

                m = 0
                do k = 1, nElemOrComp
                    if (iElementSystem(k) /= 0) then
                        m = m + 1
                        dStoichSpecies(j,m) = dStoichSpeciesCS(i,k)
                        if (cSolnPhaseTypeCS(n) == 'SUBI') then
                            s = iPhaseSublatticeCS(n)
                            c = i - nSpeciesPhaseCS(n-1)
                            if (dSublatticeChargeCS(s,2,iConstituentSublatticeCS(s,2,c)) == 0D0) then
                                dStoichSpecies(j,m) = dStoichSpecies(j,m) &
                                                                  / dSublatticeChargeCS(s,1,iConstituentSublatticeCS(s,1,c))
                            end if
                        end if
                    end if
                end do

                dStdEnergyTemp = dChemicalPotential(j)

                do k = 2, 7
                    dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(k,l2) * dGibbsCoeff(k-1)
                end do

                ! Compute additional standard molar Gibbs energy terms:
                if (dGibbsCoeffSpeciesTemp(9,l2) .EQ. 99) then
                    dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(8,l2) * dLogT
                else
                    dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(8,l2) &
                        *dTemperature**dGibbsCoeffSpeciesTemp(9,l2)
                end if

                if (dGibbsCoeffSpeciesTemp(11,l2) .EQ. 99) then
                    dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(10,l2) * dLogT
                else
                    dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(10,l2) &
                        * dTemperature**dGibbsCoeffSpeciesTemp(11,l2)
                end if

                if (dGibbsCoeffSpeciesTemp(13,l2) .EQ. 99) then
                    dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(12,l2) * dLogT
                else
                    dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(12,l2) &
                        * dTemperature**dGibbsCoeffSpeciesTemp(13,l2)
                end if

                ! If there are multiple Standard Gibbs Energy equations, check which one to use in the
                ! case of repeated temperature ranges. The rule appears to be to use the greater (less
                ! negative) energy for pure condensed phases, but just the first expression listed for
                ! solution phases. Why? Because. Thus, for now we might as well skip this block. I am
                ! therefore commenting it out.

                ! if (nGibbsEqSpecies(i) > 1) then
                !     dChemPot1 = dChemicalPotential(j)
                !     if (l1 > 0) then
                !         if (dGibbsCoeffSpeciesTemp(1,l2) == dGibbsCoeffSpeciesTemp(1,l2+1)) then
                !             l2 = l2 + 1
                !         end if
                !     else if (l1 == 0) then
                !         if (dGibbsCoeffSpeciesTemp(1,l2) == dGibbsCoeffSpeciesTemp(1,l2-1)) then
                !             l2 = l2 - 1
                !         end if
                !     end if
                !
                !     dChemicalPotential(j) = dStdEnergyTemp
                !     do k = 2, 7
                !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(k,l2) * dGibbsCoeff(k-1)
                !     end do
                !     ! Compute additional standard molar Gibbs energy terms:
                !     if (dGibbsCoeffSpeciesTemp(9,l2) .EQ. 99) then
                !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(8,l2) * dLogT
                !     else
                !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(8,l2) &
                !             *dTemperature**dGibbsCoeffSpeciesTemp(9,l2)
                !     end if
                !
                !     if (dGibbsCoeffSpeciesTemp(11,l2) .EQ. 99) then
                !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(10,l2) * dLogT
                !     else
                !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(10,l2) &
                !             * dTemperature**dGibbsCoeffSpeciesTemp(11,l2)
                !     end if
                !
                !     if (dGibbsCoeffSpeciesTemp(13,l2) .EQ. 99) then
                !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(12,l2) * dLogT
                !     else
                !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(12,l2) &
                !             * dTemperature**dGibbsCoeffSpeciesTemp(13,l2)
                !     end if
                !
                !     dChemPot2 = dChemicalPotential(j)
                !     dChemicalPotential(j) = dChemPot1 !MIN(dChemPot1,dChemPot2)
                ! end if

                ! Compute the magnetic terms (if applicable):
                if ((dGibbsMagneticCS(i,1) /= 0D0)) then
                    if (iPhase(j) == 0) then
                        call CompGibbsMagnetic(i,j)
                    else
                        call CompGibbsMagneticSolnInit(j)
                    end if
                end if

                ! Convert chemical potentials to dimensionless units:
                dChemicalPotential(j) = dChemicalPotential(j) * dTemp * DFLOAT(iParticlesPerMoleCS(i))

                ! Add pressure dependence term to the chemical potential term:
                if (iPhaseCS(i) == 1) then
                    ! For real??? Using the phase name here? I'm going to lose my mind.
                    ! if (cSolnPhaseNameCS(1) == 'gas_ideal') then
                    if (cSolnPhaseTypeCS(1) == 'IDMX') then
                        ! Note: If an ideal gas is included in a ChemSage data-file, then it is always
                        ! the first solution phase in the data-file.
                        dChemicalPotential(j) = dChemicalPotential(j) + dLogP + (0.10945D0 / dIdealConstant)
                    end if
                end if

            end do LOOP_nSpeciesCS ! End loop of species (i)
        end if
    end do LOOP_nPhasesCS

    LOOP_nPureConSpeciesCS: do i = nSpeciesPhaseCS(nSolnPhasesSysCS) + 1, nSpeciesCS
        l1 = 0
        ! Loop through the Gibbs energy equations to figure out which one to use:
        do k = 1, nGibbsEqSpecies(i)
            iCounterGibbsEqn = iCounterGibbsEqn + 1
            if ((dTemperatureForLimits <= dGibbsCoeffSpeciesTemp(1,iCounterGibbsEqn)).AND.(l1 == 0)) then
                l1 = k
            end if
        end do

        ! This species will not be considered part of the system.
        if (iSpeciesPass(i) == 0) cycle LOOP_nPureConSpeciesCS

        l2 = l1
        if (l2 == 0) l2 = nGibbsEqSpecies(i)

        l2 = l2 + iCounterGibbsEqn - nGibbsEqSpecies(i)
        j = j + 1   ! New species index

        cSpeciesName(j)            = cSpeciesNameCS(i)
        iPhase(j)                  = iPhaseCS(i)
        iParticlesPerMole(j)       = iParticlesPerMoleCS(i)
        dCoeffGibbsMagnetic(j,1:4) = dGibbsMagneticCS(i,1:4)

        m = 0
        do k = 1, nElemOrComp
            if (iElementSystem(k) /= 0) then
                m = m + 1
                dStoichSpecies(j,m) = dStoichSpeciesCS(i,k)
            end if
        end do

        dStdEnergyTemp = dChemicalPotential(j)

        do k = 2, 7
            dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(k,l2) * dGibbsCoeff(k-1)
        end do

        ! Compute additional standard molar Gibbs energy terms:
        if (dGibbsCoeffSpeciesTemp(9,l2) .EQ. 99) then
            dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(8,l2) * dLogT
        else
            dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(8,l2) &
                *dTemperature**dGibbsCoeffSpeciesTemp(9,l2)
        end if

        if (dGibbsCoeffSpeciesTemp(11,l2) .EQ. 99) then
            dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(10,l2) * dLogT
        else
            dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(10,l2) &
                * dTemperature**dGibbsCoeffSpeciesTemp(11,l2)
        end if

        if (dGibbsCoeffSpeciesTemp(13,l2) .EQ. 99) then
            dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(12,l2) * dLogT
        else
            dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(12,l2) &
                * dTemperature**dGibbsCoeffSpeciesTemp(13,l2)
        end if

        ! If there are multiple Standard Gibbs Energy equations, check which one to use in the
        ! case of repeated temperature ranges. The rule appears to be to use the greater (less
        ! negative) energy for pure condensed phases only!
        ! As of 8.0 this seems not to be happening
        ! if (nGibbsEqSpecies(i) > 1) then
        !     dChemPot1 = dChemicalPotential(j)
        !     if ((l1 > 0) .AND. (l1 < nGibbsEqSpecies(i))) then
        !         if (dGibbsCoeffSpeciesTemp(1,l2) == dGibbsCoeffSpeciesTemp(1,l2+1)) then
        !             l2 = l2 + 1
        !         end if
        !     else if (l1 == 0) then
        !         if (dGibbsCoeffSpeciesTemp(1,l2) == dGibbsCoeffSpeciesTemp(1,l2-1)) then
        !             l2 = l2 - 1
        !         end if
        !     end if
        !     dChemicalPotential(j) = dStdEnergyTemp
        !     do k = 2, 7
        !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(k,l2) * dGibbsCoeff(k-1)
        !     end do
        !     ! Compute additional standard molar Gibbs energy terms:
        !     if (dGibbsCoeffSpeciesTemp(9,l2) .EQ. 99) then
        !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(8,l2) * dLogT
        !     else
        !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(8,l2) &
        !             *dTemperature**dGibbsCoeffSpeciesTemp(9,l2)
        !     end if
        !
        !     if (dGibbsCoeffSpeciesTemp(11,l2) .EQ. 99) then
        !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(10,l2) * dLogT
        !     else
        !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(10,l2) &
        !             * dTemperature**dGibbsCoeffSpeciesTemp(11,l2)
        !     end if
        !
        !     if (dGibbsCoeffSpeciesTemp(13,l2) .EQ. 99) then
        !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(12,l2) * dLogT
        !     else
        !         dChemicalPotential(j) = dChemicalPotential(j) + dGibbsCoeffSpeciesTemp(12,l2) &
        !             * dTemperature**dGibbsCoeffSpeciesTemp(13,l2)
        !     end if
        !
        !     dChemPot2 = dChemicalPotential(j)
        !     dChemicalPotential(j) = dChemPot1!(dChemPot1+dChemPot2)/2D0
        ! end if

        ! Compute the magnetic terms (if applicable):
        if ((dGibbsMagneticCS(i,1) /= 0D0).AND.(iPhase(j) == 0)) then
            call CompGibbsMagnetic(i,j)
        end if

        ! Convert chemical potentials to dimensionless units:
        dChemicalPotential(j) = dChemicalPotential(j) * dTemp * DFLOAT(iParticlesPerMoleCS(i))

        if (iPhaseCS(i) == -1) then
            ! Explicitly set dummy species chemical potentials
            dChemicalPotential(j) = 0D0
        end if
    end do LOOP_nPureConSpeciesCS ! End loop of species (i)

    ! Update all the excess properties:
    n = 0
    nn = 0
    LOOP_SolnPhases: do i = 1, nSolnPhasesSysCS

        if ((cSolnPhaseTypeCS(i) == 'SUBL').OR.(cSolnPhaseTypeCS(i) == 'SUBLM') &
        .OR.(cSolnPhaseTypeCS(i) == 'SUBG').OR.(cSolnPhaseTypeCS(i) == 'SUBQ' ) &
        .OR.(cSolnPhaseTypeCS(i) == 'SUBI')) nCounter = nCounter + 1

        LOOP_Param: do j = nParamPhaseCS(i-1) + 1, nParamPhaseCS(i)

            ! Proceed if the parameter passed:
            IF_ParamPass: if (iParamPassCS(j) /= 0) then
                n = n + 1

                ! Convert the excess Gibbs energy parameters to dimensionless units:
                iRegularParam(n,1:nParamMax*2+3) = iRegularParamCS(j,1:nParamMax*2+3)
                cRegularParam(n) = cRegularParamCS(j)

                select case (cSolnPhaseTypeCS(i))
                    case ('QKTO', 'RKMP', 'RKMPM')

                        do k = 1, 6
                            dExcessGibbsParam(n) = dExcessGibbsParam(n) + dRegularParamCS(j,k) * dGibbsCoeff(k)
                        end do
                        dExcessGibbsParam(n) = dExcessGibbsParam(n) * dTemp

                        ! Loop through species involved in mixing parameter:
                        do k = 1, iRegularParamCS(j,1)
                            m                    = iRegularParamCS(j,k+1) + nSpeciesPhaseCS(i-1)
                            iRegularParam(n,k+1) = iSpeciesPass(m)
                        end do
                        if (iRegularParamCS(j,1) == 3) then
                            m                    = iRegularParamCS(j,iRegularParamCS(j,1)+2) + nSpeciesPhaseCS(i-1)
                            iRegularParam(n,iRegularParamCS(j,1)+2) = iSpeciesPass(m)
                        end if

                    case ('SUBG','SUBQ')

                        ! Must remove unused constituents from iRegularParam
                        iSublPhaseIndex = iPhaseSublatticeCS(i)
                        nRemove = 0
                        iRemove = 0
                        do k = nSublatticePhaseCS(iSublPhaseIndex), 1, -1
                            do l = nConstituentSublatticeCS(iSublPhaseIndex,k), 1, -1
                                if (iConstituentPass(iSublPhaseIndex,k,l) <= 0) then
                                    nRemove = nRemove + 1
                                    iRemove(nRemove) = l + ((k - 1) * nConstituentSublatticeCS(iSublPhaseIndex,1))
                                end if
                            end do
                        end do

                        do k = 1, nRemove
                            do l = 2, 5
                                if (iRegularParam(n,l) > iRemove(k)) then
                                    iRegularParam(n,l) = iRegularParam(n,l) - 1
                                end if
                            end do
                            do l = 10, 11
                                if (iRegularParam(n,l) > iRemove(k)) then
                                    iRegularParam(n,l) = iRegularParam(n,l) - 1
                                end if
                            end do
                        end do

                        ! Note that this is different for SUBG phases than QKTO, RKMP, or SUBL phases:
                        do k = 1, 4
                            dExcessGibbsParam(n) = dExcessGibbsParam(n) + dRegularParamCS(j,k) * dGibbsCoeff(k)
                        end do
                        dExcessGibbsParam(n) = dExcessGibbsParam(n) * dTemp

                    case ('SUBL', 'SUBLM')

                        do k = 1, 6
                            dExcessGibbsParam(n) = dExcessGibbsParam(n) + dRegularParamCS(j,k) * dGibbsCoeff(k)
                        end do
                        dExcessGibbsParam(n) = dExcessGibbsParam(n) * dTemp

                        ! Loop through constituents involved in mixing parameter:
                        do k = 1, iRegularParamCS(j,1)

                            ! The constituent numbering scheme from ChemSage does not consider the sublattice #, but just
                            ! a continuing count of the constituents.
                            m = iRegularParamCS(j,k+1)

                            ! Figure out the sublattice (l) and constituent (m) indices:
                            LOOP_SUBL: do s = 1, nSublatticePhaseCS(nCounter)
                                l = s
                                if (m > nConstituentSublatticeCS(nCounter,s)) then
                                    m = m - nConstituentSublatticeCS(nCounter,s)
                                else
                                    exit LOOP_SUBL
                                end if

                            end do LOOP_SUBL

                            ! Apply indexing scheme (l is sublattice index, iCounstituentPass is constituent index on
                            ! sublattice l):
                            iRegularParam(n,k+1) = (10000 * l) + iConstituentPass(nCounter,s,m)

                        end do

                        ! Loop through constituents involved in mixing parameter to see if they need to be shuffled.
                        ! ChemSage files do not order the constituents based on which ones mix.  For instance, there
                        ! may be three constituents mixing where the first constituent is on the first sublattice and
                        ! the second and third are on the second sublattice.

                        nMixSets = 0
                        k = 2
                        LOOP_SUBL_Check: do while (k <= iRegularParamCS(j,1))

                            l = MOD(iRegularParam(n,k), 10000)
                            l = (iRegularParam(n,k) - l) / 10000

                            iMixLength = 1
                            iMixStart = 0

                            LOOP_SUBL_MIXING: do ii = k + 1, iRegularParamCS(j,1) + 1
                                m = MOD(iRegularParam(n,ii), 10000)
                                m = (iRegularParam(n,ii) - m) / 10000
                                if (l == m) then
                                    iMixLength = iMixLength + 1
                                    if (ii - k == 1) then
                                        nMixSets = nMixSets + 1
                                        iMixStart = k
                                    end if
                                else
                                    exit LOOP_SUBL_MIXING
                                end if
                            end do LOOP_SUBL_MIXING

                            if (iMixLength > 1) then
                                iSUBLParamData(n,nMixSets*2) = iMixStart
                                iSUBLParamData(n,nMixSets*2+1) = iMixLength
                            end if

                            k = k + iMixLength
                        end do LOOP_SUBL_Check

                        iSUBLParamData(n,1) = nMixSets

                    case ('SUBI')
                        ! Populating iSUBIMixType from parsed CS data
                        iSUBIMixType(n) = iSUBIMixTypeCS(j)
                        ! Compute mixing terms L
                        do k = 1, 6
                            dExcessGibbsParam(n) = dExcessGibbsParam(n) + dRegularParamCS(j,k) * dGibbsCoeff(k)
                        end do

                        dExcessGibbsParam(n) = dExcessGibbsParam(n) * dTemp

                        do k = 1, iRegularParamCS(j,1)
                            ! The constituent numbering scheme from ChemSage does not consider the sublattice #, but just
                            ! a continuing count of the constituents.
                            m = iRegularParamCS(j,k+1)

                            ! Figure out the sublattice (l) and constituent (m) indices:
                            LOOP_SUBI: do s = 1, nSublatticePhaseCS(nCounter)
                                l = s
                                if (m > nConstituentSublatticeCS(nCounter,s)) then
                                    m = m - nConstituentSublatticeCS(nCounter,s)
                                else
                                    exit LOOP_SUBI
                                end if

                            end do LOOP_SUBI

                            ! Apply indexing scheme (l is sublattice index, iCounstituentPass is constituent index on
                            ! sublattice l):
                            iRegularParam(n,k+1) = (10000 * l) + iConstituentPass(nCounter,s,m)
                        end do
                end select
            end if IF_ParamPass
        end do LOOP_Param
        LOOP_MagParam: do j = nMagParamPhaseCS(i-1) + 1, nMagParamPhaseCS(i)

            ! Proceed if the parameter passed:
            IF_MagParamPass: if (iMagParamPassCS(j) /= 0) then
                nn = nn + 1
                iMagneticParam(nn,1:nParamMax*2+3) = iMagneticParamCS(j,1:nParamMax*2+3)
                dMagneticParam(nn,1:2) = dMagneticParamCS(j,1:2)

                if (cSolnPhaseTypeCS(i) == 'RKMPM') then
                    ! Loop through species involved in mixing parameter:
                    do k = 1, iMagneticParamCS(j,1)
                        m                    = iMagneticParamCS(j,k+1) + nSpeciesPhaseCS(i-1)
                        iMagneticParam(nn,k+1) = iSpeciesPass(m)
                    end do
                    if (iMagneticParamCS(j,1) == 3) then
                        m                    = iMagneticParamCS(j,iMagneticParamCS(j,1)+2) + nSpeciesPhaseCS(i-1)
                        iMagneticParam(nn,iMagneticParamCS(j,1)+2) = iSpeciesPass(m)
                    end if
                else if (cSolnPhaseTypeCS(i) == 'SUBLM') then
                    ! Loop through constituents involved in mixing parameter:
                    do k = 1, iMagneticParamCS(j,1)

                        ! The constituent numbering scheme from ChemSage does not consider the sublattice #, but just
                        ! a continuing count of the constituents.
                        m = iMagneticParamCS(j,k+1)

                        ! Figure out the sublattice (l) and constituent (m) indices:
                        LOOP_SUBLmag: do s = 1, nSublatticePhaseCS(nCounter)
                            l = s
                            if (m > nConstituentSublatticeCS(nCounter,s)) then
                                m = m - nConstituentSublatticeCS(nCounter,s)
                            else
                                exit LOOP_SUBLmag
                            end if

                        end do LOOP_SUBLmag

                        ! Apply indexing scheme (l is sublattice index, iCounstituentPass is constituent index on
                        ! sublattice l):
                        iMagneticParam(nn,k+1) = (10000 * l) + iConstituentPass(nCounter,s,m)

                    end do

                    ! Loop through constituents involved in mixing parameter to see if they need to be shuffled.
                    ! ChemSage files do not order the constituents based on which ones mix.  For instance, there
                    ! may be three constituents mixing where the first constituent is on the first sublattice and
                    ! the second and third are on the second sublattice.
                    LOOP_SUBLmag_Check: do k = 2, iMagneticParamCS(j,1)

                        l = MOD(iMagneticParam(nn,k), 10000)
                        l = (iMagneticParam(nn,k) - l) / 10000

                        m = MOD(iMagneticParam(nn,k+1), 10000)
                        m = (iMagneticParam(nn,k+1) - m) / 10000

                        ! If these two constituents are on the same sublattice and they correspond to the
                        ! first two mixing constituents, then exit:
                        if (l == m) then

                            ! Check the order of constituents:
                            if (k == 2) then
                                ! These constituent indices are correctly placed:
                                exit LOOP_SUBLmag_Check
                            elseif (k == 3) then
                                ! Shuffle the vector:
                                l = iMagneticParam(nn,2)
                                iMagneticParam(nn,2) = iMagneticParam(nn,3)
                                iMagneticParam(nn,3) = iMagneticParam(nn,4)
                                iMagneticParam(nn,4) = l

                                exit LOOP_SUBLmag_Check
                            elseif (k == 4) then
                                ! Shuffle the vector:
                                l = iMagneticParam(nn,2)
                                m = iMagneticParam(nn,3)
                                iMagneticParam(nn,2) = iMagneticParam(nn,4)
                                iMagneticParam(nn,3) = iMagneticParam(nn,5)
                                iMagneticParam(nn,4) = l
                                iMagneticParam(nn,5) = m

                                exit LOOP_SUBLmag_Check
                            else
                                ! Report an error and exit:
                                INFOThermo = 36
                                exit LOOP_SUBLmag_Check
                            end if
                        end if
                    end do LOOP_SUBLmag_Check
                end if
            end if IF_MagParamPass
        end do LOOP_MagParam
    end do LOOP_SolnPhases

    ! Update the phase index vector (iPhase):
    do i = 1, nSolnPhasesSys
        j = nSpeciesPhase(i-1) + 1
        k = nSpeciesPhase(i)
        iPhase(j:k) = i
    end do

    ! If necessary, add dummy species for each ionic phase:
    if (nChargedConstraints > 0) then
        j = 0

        do i = nSpecies, nSpecies - nChargedConstraints + 1, -1
            j = j + 1
            cSpeciesName(i)       = 'e-'
            iPhase(i)             = -1
            dChemicalPotential(i) = 5D6 / (dIdealConstant * dTemperature)
            iParticlesPerMole(i)  = 1
            dStoichSpecies(i,nElements - j + 1) = 1D0
        end do
    end if

    i = nSpeciesPhase(nSolnPhasesSys) + 1
    nDummySpecies = ABS(SUM(iPhase(i:nSpecies)))

    ! Compute the total number of atoms per formula mass of each species:
    dSpeciesTotalAtoms = SUM(ABS(dStoichSpecies),DIM=2)

    ! Compute the atomic fraction of each element in a species:
    do i = 1, nSpecies

        ! SUBG phases will be allocated with more species than appear in the data-file
        ! because it'll work with pair fractions rather than species. Therefore, dummy
        ! entries will be created with no stoichiometric values. To avoid dividing by
        ! zero, the following is included:
        !dSpeciesTotalAtoms(i) = DMAX1(dSpeciesTotalAtoms(i),0.1D0)
        if (dSpeciesTotalAtoms(i) == 0) then
            dSpeciesTotalAtoms(i) = 0.00001D0
        end if

        ! Convert chemical potentials from [J/mol] to [J/g-at]:
        dTemp                 = 1D0 / dSpeciesTotalAtoms(i)
        dChemicalPotential(i) = dChemicalPotential(i) * dTemp

        do j = 1, nElements
            dAtomFractionSpecies(i,j) = dStoichSpecies(i,j) * dTemp
        end do
    end do

    ! Store the standard molar Gibbs energies:
    do i = 1, nSpecies
        dStdGibbsEnergy(i) = dChemicalPotential(i) * dSpeciesTotalAtoms(i) / DFLOAT(iParticlesPerMole(i))
        dChemicalPotential(i) = dChemicalPotential(i) + dMagGibbsEnergy(i)
    end do

    nRemove = 0
    iRemove = 0
    LOOP_ElementUse: do i = nElements, 1, -1
        if (dElementMass(i) == 0D0) then
            do j = 1, nSpecies - nDummySpecies
                if (DABS(dStoichSpecies(j,i)) > 0D0) cycle LOOP_ElementUse
            end do
            ! If there is none of an element (i.e. it's an electron) and no species use it
            ! then set it for removal
            nRemove = nRemove + 1
            iRemove(nRemove) = i
        end if
    end do LOOP_ElementUse

    ! Do removal of unused electron elements
    if (nRemove > 0) then
        ! Create temp variables for copying
        allocate(cElementNameTemp(nElements),dMolesElementTemp(nElements))
        allocate(dAtomFractionSpeciesTemp(nSpecies,nElements),dStoichSpeciesTemp(nSpecies,nElements))
        if (nCountSublattice > 0) allocate(dStoichPairsTemp(nCountSublattice,MAXVAL(nPairsSROCS(:,1)),nElements))
        cElementNameTemp = cElementName
        dMolesElementTemp = dMolesElement
        dAtomFractionSpeciesTemp = dAtomFractionSpecies
        dStoichSpeciesTemp = dStoichSpecies
        if (nCountSublattice > 0) dStoichPairsTemp = dStoichPairs

        ! Set new number of elements
        nElements = nElements - nRemove
        nChargedConstraints = nChargedConstraints - nRemove

        ! De- and re-allocate permanent variables
        deallocate(cElementName,dMolesElement,dAtomFractionSpecies,dStoichSpecies)
        if (nCountSublattice > 0) deallocate(dStoichPairs)
        allocate(cElementName(nElements),dMolesElement(nElements))
        allocate(dAtomFractionSpecies(nSpecies,nElements),dStoichSpecies(nSpecies,nElements))
        if (nCountSublattice > 0) allocate(dStoichPairs(nCountSublattice,MAXVAL(nPairsSROCS(:,1)),nElements))

        ! Copy required elements
        k = 0
        LOOP_CopyElements: do i = 1, (nElements + nRemove)
            do j = 1, nRemove
                if (i == iRemove(j)) cycle LOOP_CopyElements
            end do
            k = k + 1
            cElementName(k) = cElementNameTemp(i)
            dMolesElement(k) = dMolesElementTemp(i)
            dAtomFractionSpecies(1:nSpecies,k) = dAtomFractionSpeciesTemp(1:nSpecies,i)
            dStoichSpecies(1:nSpecies,k) = dStoichSpeciesTemp(1:nSpecies,i)
            if (nCountSublattice > 0) dStoichPairs(1:nCountSublattice,1:MAXVAL(nPairsSROCS(:,1)),k) = &
                        dStoichPairsTemp(1:nCountSublattice,1:MAXVAL(nPairsSROCS(:,1)),i)
        end do LOOP_CopyElements

        ! Cleanup
        deallocate(cElementNameTemp,dMolesElementTemp,dAtomFractionSpeciesTemp,dStoichSpeciesTemp)
        if (nCountSublattice > 0) deallocate(dStoichPairsTemp)
    end if

    ! Store an integer vector representing the component index when the phase is ionic.
    if (allocated(iPhaseElectronID)) deallocate(iPhaseElectronID)
    allocate(iPhaseElectronID(nSolnPhasesSys))
    iPhaseElectronID = 0

    ! Loop through all solution phases:
    LOOP_EID_Phase: do i = 1, nSolnPhasesSys
        ! Check if this phase may be ionic:
        if ((cSolnPhaseType(i) == 'SUBL').OR.(cSolnPhaseType(i) == 'SUBLM')) then
            ! Loop through species in phase:
            LOOP_EID_Species: do j = nSpeciesPhase(i-1) + 1, nSpeciesPhase(i)
                ! Loop through system components represented by electrons:
                do k = nElements - nChargedConstraints + 1, nElements
                    ! If a species in this phase is represented by this electron, store in the index and cycle:

                    if (dStoichSpecies(j,k) /= 0D0) then
                        iPhaseElectronID(i) = k
                        cycle LOOP_EID_Phase
                    end if
                end do
            end do LOOP_EID_Species
        end if
    end do LOOP_EID_Phase

    return

end subroutine CompThermoData