Appendix F: FORTRAN example code (cademo1.f)

The program cademo1.f whose source code is included below is part of every FORTRAN distribution of ChemApp. When receiving ChemApp it is useful to compile and run this program, and compare its results to the master results, which for most versions of ChemApp are included on the distribution disks, and which are also printed here. The program's main purpose is to call virtually every ChemApp routine to verify that the ChemApp installation is working properly. For implementations of this program in languages other than FORTRAN (see for instance the C example in Appendix G) it also serves to verify the interface to ChemApp in this language. In contrast to the other code examples in this manual it also demonstrates how every ChemApp subroutine should be checked for the value of the error variable it returns, and act accordingly if it indicates an error. If you are using ChemApp with a language other than FORTRAN or C, your distribution disk contains a program similar to this one, even tough it might not be included in the appendix of this manual.


C Note that this program contains target calculations, C which cannot be performed with the 'light' version of ChemApp
C Please also read the notes regarding the use of ChemApp with FORTRAN, C which are contained in a separate file. C C Please direct any questions about this program or the results it C produces on your machine to: C C GTT-Technologies C Kaiserstrasse 100 C 52134 Herzogenrath C Germany C C Phone: +49-2407-59533 C Fax: +49-2407-59661 C E-mail: support@gtt-technologies.de C WWW: http://www.gtt-technologies.de/

C --------------------------------------------------------------------- C Subroutine : CABORT C --------------------------------------------------------------------- C Subject : If a ChemApp error occurs, this function reports the error C number and the routine it occurred in before it exits the C program. C C --------------------------------------------------------------------- SUBROUTINE CABORT(LINENO, SRNAME, NOERR) IMPLICIT NONE INTEGER LINENO CHARACTER SRNAME*(*) INTEGER NOERR
41 FORMAT(1X,A,I4,2A)
WRITE(UNIT=*,FMT=41) 'ChemApp error no. ',NOERR, * ' occurred when calling ', srname WRITE(UNIT=*,FMT=41) 'Aborting on line ',LINENO,' of ', 'CADEMO1' STOP END

C --------------------------------------------------------------------- C Function: STRLEN C --------------------------------------------------------------------- C Subject : The support function STRLEN returns the total number of C characters in a string, not counting trailing blanks C C --------------------------------------------------------------------- INTEGER FUNCTION STRLEN(INSTR) CHARACTER INSTR*(*)
STRLEN = LEN(INSTR) DO WHILE(INSTR(STRLEN:STRLEN) .EQ. ' ') STRLEN = STRLEN - 1 IF (STRLEN .EQ. 1) RETURN ENDDO RETURN END

PROGRAM CADEMO IMPLICIT NONE
C Declaration of variables to be used in the program
C Declarations of integer variables for use with ChemApp:
C Declaring the return variable for the function STRLEN INTEGER STRLEN
C An all-purpose integer variable INTEGER I
C Integer variables which hold the index number of phases, phase C constituents, and sublattices INTEGER INDEXP, INDEXC, INDEXL
C The standard error code return variable INTEGER NOERR
C A variable that is set to the version number of ChemApp INTEGER CAVERS
C This variable is used in conjunction with TQLITE to determine whether C the program is linked to the regular or the 'light' version of C ChemApp. In the latter case, target calculations are left out. INTEGER ISLITE
C A number of integer variables to return the information from TQSIZE INTEGER LIA,LIB,LIC,LID,LIE,LIF,LIG,LIH,LII,LIJ,LIK
C A similar set of integer variables to return the information from C TQUSED INTEGER LDIA,LDIB,LDIC,LDID,LDIE,LDIF,LDIG,LDIH,LDII,LDIJ,LDIK
C An integer that will contain the index number of a condition set INTEGER NUMCON
C An integer that will contain the number of phases in the data-file C loaded INTEGER NPHASE
C An integer that will hold the number of phase constituents in the gas C phase INTEGER NPCGAS
C An integer that will contain the number of system components INTEGER NSCOM
C Integers that will contain the number of sublattices, and the number of C sublattice contituents of a particular sublattice INTEGER NSL, NSLCON
C Two integers used for one-dimensional phase mapping calculations, one C to indicate whether we need to make more calls to TQMAP/TQMAPL INTEGER ICONT
C And another one to keep track of the number of results we obtained INTEGER RESNO
C An integer that holds the FORTRAN unit number under which the C thermochemical data-files will be opened INTEGER UNITNO
C An integer that holds the FORTRAN unit number to which the C ChemApp error messages are written by default INTEGER EUNIT
C Variables to store information from various transparent file header C fields INTEGER TFHVER, TFHVNW(3), TFHVNR(3), TFHDTC(6), TFHDTE(6)
C These integers will be used to store the HASP dongle id, plus the C ChemApp expiration month and year INTEGER HASPID, EDMON, EDYEAR
C Declarations of real variables for use with ChemApp:
C An all-purpose real variable DOUBLE PRECISION D1
C An array of two real variables... DOUBLE PRECISION DARR2(2)
C ...and one of 30 DOUBLE PRECISION DARR30(30)
C A variable to hold the molecular mass of a system component (see C TQSTSC) or a phase constituent (see TQSTPC) DOUBLE PRECISION WMASS
C This particular array of two real variables will contain the C temperature and pressure of a stream to be defined with TQSTTP DOUBLE PRECISION TP(2)

C Declarations of string variables for use with ChemApp:

C Every 'name' used in conjuction with ChemApp, whether it is the name C of a phase, a system component, a phase constituent, or a stream, can C be up to 24 characters in length. Exceptions are the subroutine TQERR, C as well as the subroutines TQGTRH, TQGTID, TQGTNM, TQGTPI, and TQWSTR, C which all pass or return strings longer than 24 characters.
C An all-purpose string for ChemApp CHARACTER DSTR*24
C Character strings that will hold names of phases, phase constituents, C system components, mixture model names, as well as sublattices and C their constituents
CHARACTER*24 PNAME, CNAME, MNAME
C The following variable is an array of 3 strings. It will be used with C TQCSC, which is the only ChemApp subroutine that takes a character C array as input. CHARACTER NEWSC(3)*24
C The following character array will hold a ChemApp error message. In C this case, it will be used to retrieve the copyright string. CHARACTER MESS(3)*80
C The next three character strings are used to hold information C about the copy of ChemApp used. CHARACTER ID*255, NAME*80, PID*24
C Variables to store information from various transparent file header C fields CHARACTER TFHID*255, TFHUSR*80,TFHNWP*40,TFHNRP*40,TFHREM*80
C This character variable will be used to store the HASP dongle type CHARACTER HASPT*5
C Format declarations 51 FORMAT(1X,A) 52 FORMAT(1X,A,I4) 53 FORMAT(1X,A,//) 54 FORMAT(1X,A,I4,//) 55 FORMAT(1X,//,1X,A,//) 56 FORMAT(1X,A,I3,A,I3,A) 57 FORMAT(1X,A,G12.5) 58 FORMAT(1X,A,G12.5,//) 59 FORMAT(1X,A,F8.2,A) 60 FORMAT(1X,A,I3,A) 61 FORMAT(1X,A,I4.4,5(A1,I2.2)) 62 FORMAT(1X,A,I2,A,I4,//)
C Print the program title WRITE(UNIT=*,FMT=51) 'This is cademo1, a FORTRAN program ' // * 'to test ChemApp and to demonstrate ' WRITE(UNIT=*,FMT=53) 'the use of every ChemApp subroutine.'

Output:
 This is cademo1, a FORTRAN program to test ChemApp and to demonstrate 
 the use of every ChemApp subroutine.

C Initialise ChemApp CALL TQINI(NOERR)
C After all following invocations of ChemApp subroutines, the error C variable will be checked for a non-zero value, in which case an C error has occurred and the program is aborted IF (NOERR .NE. 0) CALL CABORT(224,'TQINI',NOERR)

C Print the copyright string. The way to do this is to first C call the subroutine TQCPRT, which copies the copyright string into C the ChemApp-internal error buffer string: CALL TQCPRT(NOERR) IF (NOERR .NE. 0) CALL CABORT(231,'TQCPRT',NOERR)
C Then do what is usually done to retrieve an error message: C call TQERR to retrieve it. CALL TQERR(MESS,NOERR) IF (NOERR .NE. 0) CALL CABORT(236,'TQERR',NOERR)
C To print out the message, print all the three strings of C MESS : WRITE(UNIT=*,FMT='(3(A,/))') MESS(1), MESS(2), MESS(3)

Output:
 This program contains ChemApp                                                  
 Copyright GTT-Technologies, Kaiserstrasse 100, D-52134 Herzogenrath, Germany   
 http://www.gtt-technologies.de


C Get version number of ChemApp CALL TQVERS(CAVERS,NOERR) IF (NOERR .NE. 0) CALL CABORT(246,'TQVERS',NOERR) WRITE(UNIT=*,FMT=54) 'ChemApp version is: ', CAVERS
C Get array sizes CALL TQSIZE(LIA,LIB,LIC,LID,LIE,LIF,LIG,LIH,LII,LIJ,LIK,NOERR) IF (NOERR .NE. 0) CALL CABORT(251,'TQSIZE',NOERR) WRITE(UNIT=*,FMT=51) 'Internal array sizes of this version '// * 'of ChemApp,' WRITE(UNIT=*,FMT=51) 'maximum number of:' WRITE(UNIT=*,FMT=52) 'constituents: ',LIA WRITE(UNIT=*,FMT=52) 'system components: ',LIB WRITE(UNIT=*,FMT=52) 'mixture phases: ',LIC WRITE(UNIT=*,FMT=52) 'excess Gibbs energy '// * 'coefficients for a mixture phase: ',LID WRITE(UNIT=*,FMT=52) 'excess magnetic '// * 'coefficients for a mixture phase: ',LIE WRITE(UNIT=*,FMT=52) 'sublattices for a '// * 'mixture phase: ',LIF WRITE(UNIT=*,FMT=52) 'constituents of a '// * 'sublattice: ',LIG WRITE(UNIT=*,FMT=51) 'oxide constituents of a phase described '// * 'by the Gaye-Kapoor-Frohberg ' WRITE(UNIT=*,FMT=52) ' or modified quasichemical '// * 'formalisms: ',LIH WRITE(UNIT=*,FMT=52) 'Gibbs energy/heat '// * 'capacity equations for a constituent: ',LII WRITE(UNIT=*,FMT=52) 'Gibbs energy/heat '// * 'capacity equations: ',LIJ WRITE(UNIT=*,FMT=54) 'constituents with '// * 'P,T-dependent molar volumes: ',LIK

Output:
 ChemApp version is:  502


Internal array sizes of this version of ChemApp, maximum number of: constituents: 300 system components: 15 mixture phases: 20 excess Gibbs energy coefficients for a mixture phase: 199 excess magnetic coefficients for a mixture phase: 99 sublattices for a mixture phase: 5 constituents of a sublattice: 12 oxide constituents of a phase described by the Gaye-Kapoor-Frohberg or modified quasichemical formalisms: 15 Gibbs energy/heat capacity equations for a constituent: 20 Gibbs energy/heat capacity equations: 900 constituents with P,T-dependent molar volumes: 300


C Get value for input/output option
C Determine which FORTRAN unit is used by default by TQRFIL for C reading the thermochemical data-file, and use the returned number to C subsequently open and close the data-file CALL TQGIO('FILE ', UNITNO, NOERR) IF (NOERR .NE. 0) CALL CABORT(285,'TQGIO',NOERR) WRITE(UNIT=*,FMT=54) 'The thermochemical data will be read '// * 'from the file associated with unit ', UNITNO

Output:
 The thermochemical data will be read from the file associated with unit   10


C Open ASCII data-file for reading CALL TQOPNA('cosi.dat', UNITNO, NOERR) IF (NOERR .NE. 0) CALL CABORT(293,'TQOPNA',NOERR)
C Read data-file CALL TQRFIL(NOERR) IF (NOERR .NE. 0) CALL CABORT(297,'TQRFIL',NOERR)
C Close data-file CALL TQCLOS(UNITNO, NOERR) IF (NOERR .NE. 0) CALL CABORT(301,'TQCLOS',NOERR)

C TQUSED is used to find out to what extent the thermochemical data C storage space in the internal arrays of ChemApp (as reported by C TQSIZE) are occupied by the currently loaded thermochemical system. CALL TQUSED(LDIA,LDIB,LDIC,LDID,LDIE,LDIF,LDIG,LDIH,LDII, * LDIJ,LDIK,NOERR) IF (NOERR .NE. 0) CALL CABORT(309,'TQUSED',NOERR)
C As an example, check how many Gibbs energy/heat capacity equations C are currently being used. WRITE(UNIT=*,FMT=52) 'Actual number of Gibbs energy/heat ' // * 'capacity equations used (TQUSED):', LDIJ WRITE(UNIT=*,FMT=52) 'Maximum number in this version ' // * 'of ChemApp (TQSIZE):', LIJ

Output:
 Actual number of Gibbs energy/heat capacity equations used (TQUSED):  35
 Maximum number in this version of ChemApp (TQSIZE): 900


C Get system units CALL TQGSU('Pressure ',DSTR,NOERR) IF (NOERR .NE. 0) CALL CABORT(322,'TQGSU',NOERR) WRITE(UNIT=*,FMT=51) 'Pressure unit: ' // DSTR
CALL TQGSU('Volume ',DSTR,NOERR) IF (NOERR .NE. 0) CALL CABORT(326,'TQGSU',NOERR) WRITE(UNIT=*,FMT=51) 'Volume unit: ' // DSTR
CALL TQGSU('Temperature ',DSTR,NOERR) IF (NOERR .NE. 0) CALL CABORT(330,'TQGSU',NOERR) WRITE(UNIT=*,FMT=51) 'Temperature unit: ' // DSTR
CALL TQGSU('Energy ',DSTR,NOERR) IF (NOERR .NE. 0) CALL CABORT(334,'TQGSU',NOERR) WRITE(UNIT=*,FMT=51) 'Energy unit: ' // DSTR
C Change system unit C Here we change the unit for the quantity 'amount' to gram, C mainly so that when we call TQSTSC further down, we get the C molecular mass expressed in the unit g/mol CALL TQCSU('Amount ','gram ', NOERR) IF (NOERR .NE. 0) CALL CABORT(342,'TQCSU',NOERR) WRITE(UNIT=*,FMT=53) 'Amount unit set to gram'

Output:
 Pressure unit: bar                     
 Volume unit: dm3                     
 Temperature unit: K                       
 Energy unit: J                       
 Amount unit set to gram



C Get number of system components CALL TQNOSC(NSCOM,NOERR) IF (NOERR .NE. 0) CALL CABORT(350,'TQNOSC',NOERR) WRITE(UNIT=*,FMT=52) 'Number of system components: ',NSCOM

C Get stoichiometry of system component C Note that the array used (DARR30) is much bigger than needed, C which doesn't matter to ChemApp, as long as it is big enough. CALL TQSTSC(1,DARR30,WMASS,NOERR) IF (NOERR .NE. 0) CALL CABORT(358,'TQSTSC',NOERR)
C Here the array elements are printed, which describe the stoichiometry C of the system component. This particular piece of information is C not too useful though, since the array contains all zeros, except C for a '1.0' at the position of the system component we specified. WRITE(UNIT=*,FMT='(1X,A,3F4.1)') 'Stoichiometry of system '// * 'component 1: ', * DARR30(1), DARR30(2), DARR30(3)
C Had we not changed the amount unit from the default value 'mol' to C 'gram', we would receive a value of 1.0 for wmass. The reason is that C the molecular mass is always output using the unit [current amount C unit]/mol. If the current amount unit happens to be mol, the value for C wmass will be 1.0. WRITE(UNIT=*,FMT=57) 'Molecular mass in g/mol: ', WMASS

C Get name of system component #1 CALL TQGNSC(1,DSTR,NOERR) IF (NOERR .NE. 0) CALL CABORT(378,'TQGNSC',NOERR) WRITE(UNIT=*,FMT=51) 'Name of system component 1: ' // DSTR
C Get index number of system component C This is the reverse of the above, using the name of the C system component to get its index number CALL TQINSC(DSTR, I, NOERR) IF (NOERR .NE. 0) CALL CABORT(385,'TQINSC',NOERR) WRITE(UNIT=*,FMT=54) 'Index number of system component '// * DSTR // ' is: ', I

Output:
 Number of system components:    3
 Stoichiometry of system component 1:  1.0 0.0 0.0
 Molecular mass in g/mol:   12.011    
 Name of system component 1: C                       
 Index number of system component C                        is:    1


C Change system components C TQCSC is the only ChemApp subroutine which takes a character C array as input. The array contains as many strings as there C are system components in the data-file. The system components are C changed to SiO, SiC, and CO.
NEWSC(1) = 'SiO' NEWSC(2) = 'SiC' NEWSC(3) = 'CO'
CALL TQCSC(NEWSC,NOERR) IF (NOERR .NE. 0) CALL CABORT(402,'TQCSC',NOERR) WRITE(UNIT=*,FMT=51) 'System components changed to SiO-SiC-CO.'

C Print names of new system components
DO I=1, NSCOM
CALL TQGNSC(I,DSTR,NOERR) IF (NOERR .NE. 0) CALL CABORT(411,'TQGNSC',NOERR) WRITE(UNIT=*,FMT='(1X,A,I4,A)') 'Name of new system '// * 'component ', I, ': '// DSTR ENDDO
C Get stoichiometry of system component #1 CALL TQSTSC(1,DARR30,D1,NOERR) IF (NOERR .NE. 0) CALL CABORT(418,'TQSTSC',NOERR)
C Note that the stoichiometry of the system component #1 does not C change, as compared to the first call to TQSTSC before the C system components have been altered:
WRITE(UNIT=*,FMT='(1X,A,3F4.1)') 'Stoichiometry of system '// * 'component 1: ', * DARR30(1), DARR30(2), DARR30(3)
C What changes of course is the molecular mass: WRITE(UNIT=*,FMT=58) 'Molecular mass in g/mol: ', D1

Output:
 System components changed to SiO-SiC-CO.
 Name of new system component    1: SiO                     
 Name of new system component    2: SiC                     
 Name of new system component    3: CO                      
 Stoichiometry of system component 1:  1.0 0.0 0.0
 Molecular mass in g/mol:   44.085



C Get number of phases CALL TQNOP(NPHASE,NOERR) IF (NOERR .NE. 0) CALL CABORT(436,'TQNOP',NOERR) WRITE(UNIT=*,FMT=52) 'Number of phases: ',NPHASE

C Get name of phase #1 CALL TQGNP(1,DSTR,NOERR) IF (NOERR .NE. 0) CALL CABORT(442,'TQGNP',NOERR) WRITE(UNIT=*,FMT=51) 'Name of phase 1: ' // DSTR

C Get index number of phase C This is the reverse of the above, using the name of the C phase to get its index number CALL TQINP(DSTR,I,NOERR) IF (NOERR .NE. 0) CALL CABORT(450,'TQINP',NOERR) WRITE(UNIT=*,FMT=52) 'Index number of phase ' // DSTR // * ' is: ', I

C Get model name of phase CALL TQMODL(1, DSTR, NOERR) IF (NOERR .NE. 0) CALL CABORT(457,'TQMODL',NOERR) WRITE(UNIT=*,FMT=51) 'Model name of phase 1: ' // DSTR
C Whereas above we asked for the model name of the first phase in C the data-file, which is normally a mixture phase (gas phase), C we now ask for the model name of the last phase in the C data-file (NPHASE, as determined by a call to TQNOP). This is C typically a stoichiometric condensed phase, in which case the C string 'PURE' is returned as the model name CALL TQMODL(NPHASE, DSTR, NOERR) IF (NOERR .NE. 0) CALL CABORT(467,'TQMODL',NOERR) WRITE(UNIT=*,FMT='(1X,A,I2,A,//)') 'Model name of phase ', * nphase,': ' // DSTR

C Get number of phase constituents of phase #1 CALL TQNOPC(1, NPCGAS, NOERR) IF (NOERR .NE. 0) CALL CABORT(474,'TQNOPC',NOERR) WRITE(UNIT=*,FMT=52) 'Number of phase constituents in phase 1: ', * npcgas
C Get name of phase constituent CALL TQGNPC(1, 1, DSTR, NOERR) IF (NOERR .NE. 0) CALL CABORT(480,'TQGNPC',NOERR) WRITE(UNIT=*,FMT=51) 'Name of phase constituent 1 of phase 1: '// * dstr

C Get index number of phase constituent C This is the reverse of the above, using the name of the phase C constituent to get its index number CALL TQINPC(DSTR, 1, I, NOERR) IF (NOERR .NE. 0) CALL CABORT(489,'TQINPC',NOERR) WRITE(UNIT=*,FMT=52) 'Index number of phase constituent ' // * DSTR // ' is: ', I
C Check whether the phase constituent can be used as incoming species C (note that there are only few cases, and only for certain mixture C phase models, where this might not be permitted) CALL TQPCIS(1, 1, I, NOERR) IF (NOERR .NE. 0) CALL CABORT(497,'TQPCIS',NOERR) IF (I .EQ. 0) THEN WRITE(*,FMT=51) 'Usage of ' // DSTR // * ' as incoming species is not permitted.' ELSE WRITE(*,FMT=51) 'Usage of ' // DSTR // * ' as incoming species is o.k.' ENDIF
C Get stoichiometry of phase constituent C TQSTPC is a subroutine similar to TQSTSC, except that with C TQSTPC, the returned stoichiometry array contains valuable C information. This stoichiometry array represents the C stoichiometry matrix for a phase constituent, as described in C the ChemApp manual CALL TQSTPC(1, 1, DARR30, WMASS, NOERR) IF (NOERR .NE. 0) CALL CABORT(513,'TQSTPC',NOERR)
WRITE(UNIT=*,FMT='(1X,A,3F4.1)') 'Stoichiometry of phase '// * 'constituent 1 in phase 1 is: ', * DARR30(1), DARR30(2), DARR30(3) WRITE(UNIT=*,FMT=58) 'Molecular mass: ', WMASS

Output:
 Number of phases:    8
 Name of phase 1: GAS                     
 Index number of phase GAS                      is:    1
 Model name of phase  1: IDMX                    
 Model name of phase  8: PURE                    


Number of phase constituents in phase 1: 15 Name of phase constituent 1 of phase 1: C Index number of phase constituent C is: 1 Usage of C as incoming species is o.k. Stoichiometry of phase constituent 1 in phase 1 is: -0.5 0.5 0.5 Molecular mass: 12.011

C The system components are changed back to the original set of C, O, C and Si, since we don't need the modified set of system components C (SiO, SiC, and CO) any longer. NEWSC(1) = 'C' NEWSC(2) = 'O' NEWSC(3) = 'Si'
CALL TQCSC(NEWSC,NOERR) IF (NOERR .NE. 0) CALL CABORT(529,'TQCSC',NOERR) WRITE(UNIT=*,FMT=51) 'System components changed back to C-O-Si.'

Output:
 System components changed back to C-O-Si.


C Changing the status of phases and phase constituents C Change status of phase #1 CALL TQCSP(1, 'eliminated ', NOERR) IF (NOERR .NE. 0) CALL CABORT(537,'TQCSP',NOERR) WRITE(UNIT=*,FMT=51) 'Status of phase 1 set to ''eliminated'''
C Get status of phase #1 CALL TQGSP(1, DSTR, NOERR) IF (NOERR .NE. 0) CALL CABORT(542,'TQGSP',NOERR) WRITE(UNIT=*,FMT=51) 'Status of phase 1 is: ' // DSTR
C Change status of phase #1 CALL TQCSP(1, 'entered ', NOERR) IF (NOERR .NE. 0) CALL CABORT(547,'TQCSP',NOERR) WRITE(UNIT=*,FMT=51) 'Status of phase 1 set to ''entered'''
C Get status of phase #1 CALL TQGSP(1, DSTR, NOERR) IF (NOERR .NE. 0) CALL CABORT(552,'TQGSP',NOERR) WRITE(UNIT=*,FMT=51) 'Status of phase 1 is: ' // DSTR

Output:
 Status of phase 1 set to 'eliminated'
 Status of phase 1 is: ELIMINATED              
 Status of phase 1 set to 'entered'
 Status of phase 1 is: ENTERED


C Change status of phase constituent #1 of phase #1 CALL TQCSPC(1, 1, 'dormant ', NOERR) IF (NOERR .NE. 0) CALL CABORT(559,'TQCSPC',NOERR) WRITE(UNIT=*,FMT=51) 'Status of constituent 1 of '// * 'phase 1 set to ''dormant'''
C Get status of phase constituent #1 of phase #1 CALL TQGSPC(1, 1, DSTR, NOERR) IF (NOERR .NE. 0) CALL CABORT(565,'TQGSPC',NOERR) WRITE(UNIT=*,FMT=51) 'Status of constituent 1 of ' // * 'phase 1 is: '//DSTR
C Change status of phase constituent #1 of phase #1 CALL TQCSPC(1, 1, 'entered ', NOERR) IF (NOERR .NE. 0) CALL CABORT(571,'TQCSPC',NOERR) WRITE(UNIT=*,FMT=51) 'Status of constituent 1 of phase 1 set '// * 'to ''entered'''
C Get status of phase constituent #1 of phase #1 CALL TQGSPC(1, 1, DSTR, NOERR) IF (NOERR .NE. 0) CALL CABORT(577,'TQGSPC',NOERR) WRITE(UNIT=*,FMT=53) 'Status of constituent 1 of ' // * 'phase 1 is: '//DSTR

Output:
 Status of constituent 1 of phase 1 set to 'dormant'
 Status of constituent 1 of phase 1 is: DORMANT                 
 Status of constituent 1 of phase 1 set to 'entered'
 Status of constituent 1 of phase 1 is: ENTERED

WRITE(UNIT=*,FMT=53) 'Calculations using global conditions:'
C Set equilibrium conditions
C When setting equilibrium conditions with TQSETC, the currently C set units for the various quantities (pressure, amount, etc.) C apply. Thus, in the following call, the incoming amount ('ia') C of phase constituent #4 of phase #1 is set to 1.0 gram. CALL TQSETC('ia ', 1, 4, 1.D0,NUMCON,NOERR) IF (NOERR .NE. 0) CALL CABORT(591,'TQSETC',NOERR)
C Change system unit C Now we change the amount unit back to 'mol'. This does not C affect any previously set values. The amount of phase C constituent #4 of phase #1 we set above is still 1.0 gram, and C is not changed to 1.0 mol! CALL TQCSU('Amount ','mol ', NOERR) IF (NOERR .NE. 0) CALL CABORT(599,'TQCSU',NOERR) WRITE(UNIT=*,FMT=53) 'Amount unit set to mol'

Output:
 Calculations using global conditions:


Amount unit set to mol

C Since we changed the amount unit to 'mol', the call to TQSETC C below inputs 3.0 mol of constituent #12 of phase #1 CALL TQSETC('ia ', 1, 12, 3.D0,NUMCON,NOERR) IF (NOERR .NE. 0) CALL CABORT(606,'TQSETC',NOERR)
C Now we input 8 mol of constituent #8 of phase #1 CALL TQSETC('ia ', 1, 8, 2.D0,NUMCON,NOERR) IF (NOERR .NE. 0) CALL CABORT(610,'TQSETC',NOERR)
C Remove equilibrium condition C This call to TQREMC demonstrates how the parameter 'NUMCON' is C used, which is returned with every call to TQSETC: The C previously defined condition is removed again CALL TQREMC(NUMCON, NOERR) IF (NOERR .NE. 0) CALL CABORT(617,'TQREMC',NOERR)
C Set temperature C Here TQSETC is used to enter the temperature. Since we didn't change C the default temperature unit, it is still Kelvin ('K') CALL TQSETC('t ', 0, 0, 1800.D0,NUMCON,NOERR) IF (NOERR .NE. 0) CALL CABORT(623,'TQSETC',NOERR)

C Change limit of target variables C The next call demonstrates how TQCLIM is used. Normally though C calls to TQCLIM are only necessary under special circumstances CALL TQCLIM('plow ', 1.0D-49, NOERR) IF (NOERR .NE. 0) CALL CABORT(630,'TQCLIM',NOERR)

C Display present settings C TQSHOW outputs the conditions which are active at this point WRITE(UNIT=*,FMT=51) 'Currently active conditions:' WRITE(UNIT=*,FMT=51) '**** Begin output table produced by TQSHOW' CALL TQSHOW(NOERR) WRITE(UNIT=*,FMT=53) '**** End output table produced by TQSHOW' IF (NOERR .NE. 0) CALL CABORT(639,'TQSHOW',NOERR)

Output:
 Currently active conditions:
 **** Begin output table produced by TQSHOW

TARGET LIMITS: Pressure/bar : 1.00000E-49 1.00000E+07 Volume/dm3 : 1.00000E-07 1.00000E+50 Temperature/K: 298.15 6000.00
T = 1800.00 K P = 1.00000E+00 bar
STREAM CONSTITUENT CO/GAS/ AMOUNT/mol = 3.57010E-02 Si3/GAS/ AMOUNT/mol = 3.00000E+00 **** End output table produced by TQSHOW



C Calculate equilibrium C A simple call to TQCE or TQCEL is sufficient to calculate the C equilibrium. All input parameters to TQCE and TQCEL are only C significant if target calculations are performed. DARR2(1)=0.0 CALL TQCE(' ',0,0,DARR2,NOERR) IF (NOERR .NE. 0) CALL CABORT(650,'TQCE',NOERR)

C Calculate and list equilibrium C TQCEL calculates the equilibrium just like TQCE, with the only C difference that a ChemSage-type result table is written to the C file/unit associated with 'LIST' (see TQGIO/TQCIO) WRITE(UNIT=*,FMT=51) '**** Begin output table produced by TQCEL' CALL TQCEL(' ',0,0,DARR2,NOERR) WRITE(UNIT=*,FMT=53) '**** End output table produced by TQCEL' IF (NOERR .NE. 0) CALL CABORT(660,'TQCEL',NOERR)

Output:
 **** Begin output table produced by TQCEL

T = 1800.00 K P = 1.00000E+00 bar V = 0.00000E+00 dm3
STREAM CONSTITUENTS AMOUNT/mol CO/GAS/ 3.5701E-02 Si3/GAS/ 3.0000E+00
EQUIL AMOUNT MOLE FRACTION FUGACITY PHASE: GAS mol bar SiO 0.0000E+00 9.1172E-01 4.1405E-02 CO 0.0000E+00 8.8214E-02 4.0062E-03 Si 0.0000E+00 6.7137E-05 3.0490E-06 CO2 0.0000E+00 1.1335E-06 5.1478E-08 Si2C 0.0000E+00 7.0003E-07 3.1792E-08 Si2 0.0000E+00 4.5031E-07 2.0451E-08 SiO2 0.0000E+00 1.1236E-07 5.1026E-09 Si3 0.0000E+00 1.0170E-07 4.6188E-09 SiC 0.0000E+00 7.3913E-12 3.3567E-13 O 0.0000E+00 7.2737E-12 3.3033E-13 C 0.0000E+00 1.2567E-13 5.7072E-15 O2 0.0000E+00 1.6031E-16 7.2803E-18 C2 0.0000E+00 9.7381E-17 4.4225E-18 C3 0.0000E+00 4.0495E-17 1.8391E-18 O3 0.0000E+00 8.7034E-33 3.9526E-34 TOTAL: 0.0000E+00 1.0000E+00 4.5415E-02 mol ACTIVITY Si 8.9464E+00 1.0000E+00 SiC 3.5701E-02 1.0000E+00 SiO2(cristobalite) 1.7851E-02 1.0000E+00 SiO2(tridymite) 0.0000E+00 9.9970E-01 SiO2(liquid) 0.0000E+00 9.3907E-01 SiO2(quartz) T 0.0000E+00 9.1491E-01 C 0.0000E+00 2.2526E-02 ******************************************************************** Cp_EQUIL H_EQUIL S_EQUIL G_EQUIL V_EQUIL J.K-1 J J.K-1 J dm3 ******************************************************************** 2.46520E+02 7.84578E+05 8.41507E+02 -7.30135E+05 0.00000E+00
Mole fraction of system components: GAS C 4.4110E-02 O 4.9998E-01 Si 4.5591E-01
Data on 1 constituent marked with 'T' are extrapolated outside their valid temperature range **** End output table produced by TQCEL


C Once an equilibrium has been calculated using TQCE/TQCEL, subsequent C calculations can be performed using TQCEN/TQCENL, which speeds up the C equilibrium calculation by taking results from the previous C calculation as initial estimates CALL TQSETC('t ', 0, 0, 1850.D0,NUMCON,NOERR) IF (NOERR .NE. 0) CALL CABORT(669,'TQSETC',NOERR)
CALL TQCEN(' ', 0, 0, DARR2, NOERR) IF (NOERR .NE. 0) CALL CABORT(672,'TQCEN',NOERR)
CALL TQSETC('t ', 0, 0, 1900.D0,NUMCON,NOERR) IF (NOERR .NE. 0) CALL CABORT(675,'TQSETC',NOERR)
C Similar to TQCEL, TQCENL also provides a ChemSage output table WRITE(UNIT=*,FMT=51) '**** Begin output table produced by TQCENL' CALL TQCENL(' ', 0, 0, DARR2, NOERR) WRITE(UNIT=*,FMT=53) '**** End output table produced by TQCENL' IF (NOERR .NE. 0) CALL CABORT(681,'TQCENL',NOERR)


C Change value of input/output option C Redirect the output from TQCEL to a file using TQCIO with the C option 'LIST'. C The following call to TQCIO will redirect the output of ChemApp C subroutines like TQCEL and TQSHOW to the FORTRAN unit number C 21. Note that only certain values for the unit number are C permitted (see the ChemApp manual entry for TQCIO). CALL TQCIO('LIST ', 21, NOERR) IF (NOERR .NE. 0) CALL CABORT(693,'TQCIO',NOERR)
C Note that at this point, no file has been associated with unit C number 21 yet! This is done with the subsequent call to TQOPEN:
CALL TQOPEN('result', 21, NOERR) IF (NOERR .NE. 0) CALL CABORT(699,'TQOPEN',NOERR)
C TQWSTR can be used to write user-defined text to the unit numbers C associated with 'LIST' and 'ERROR'. This subroutine is especially C useful for non-FORTRAN programs. CALL TQWSTR('LIST', 'Output from TQCEL (ChemSage result table):', * NOERR) IF (NOERR .NE. 0) CALL CABORT(706,'TQWSTR',NOERR)

WRITE(UNIT=*,FMT=53) 'Output of next call to TQCEL will be '// * ' written to file ''result'''
C Now the ChemSage-type result table output by the following C call to TQCEL is written to the file 'result': CALL TQCEL(' ',0,0,DARR2,NOERR) IF (NOERR .NE. 0) CALL CABORT(715,'TQCEL',NOERR)

C Close the file associated with unit number 21 again... CALL TQCLOS(21, NOERR) IF (NOERR .NE. 0) CALL CABORT(720,'TQCLOS',NOERR)
C ...and redirect the output of routines like TQCEL and TQSHOW C back to unit number 6, which is the default and associated with C the standard output unit: CALL TQCIO('LIST ', 6, NOERR) IF (NOERR .NE. 0) CALL CABORT(726,'TQCIO',NOERR) WRITE(UNIT=*,FMT=53) 'Output of subsequent calls to TQCEL/'// * 'TQSHOW/TQMAPL will go to UNIT 6 again'

Output:
 **** Begin output table produced by TQCENL

T = 1900.00 K P = 1.00000E+00 bar V = 0.00000E+00 dm3
STREAM CONSTITUENTS AMOUNT/mol CO/GAS/ 3.5701E-02 Si3/GAS/ 3.0000E+00
EQUIL AMOUNT MOLE FRACTION FUGACITY PHASE: GAS mol bar SiO 0.0000E+00 8.5125E-01 1.2301E-01 CO 0.0000E+00 1.4866E-01 2.1481E-02 Si 0.0000E+00 8.4025E-05 1.2142E-05 CO2 0.0000E+00 3.7718E-06 5.4503E-07 Si2C 0.0000E+00 1.4416E-06 2.0832E-07 Si2 0.0000E+00 7.4859E-07 1.0817E-07 SiO2 0.0000E+00 2.7664E-07 3.9975E-08 Si3 0.0000E+00 1.6036E-07 2.3173E-08 SiC 0.0000E+00 3.6804E-11 5.3182E-12 O 0.0000E+00 2.9557E-11 4.2711E-12 C 0.0000E+00 7.5814E-13 1.0955E-13 O2 0.0000E+00 1.3988E-15 2.0213E-16 C2 0.0000E+00 1.3301E-15 1.9220E-16 C3 0.0000E+00 7.2884E-16 1.0532E-16 O3 0.0000E+00 6.6967E-31 9.6770E-32 TOTAL: 0.0000E+00 1.0000E+00 1.4450E-01 mol ACTIVITY Si 8.9464E+00 1.0000E+00 SiC 3.5701E-02 1.0000E+00 SiO2(cristobalite) 1.7851E-02 1.0000E+00 SiO2(tridymite) 0.0000E+00 9.9918E-01 SiO2(liquid) 0.0000E+00 9.7125E-01 SiO2(quartz) T 0.0000E+00 9.0635E-01 C 0.0000E+00 3.4700E-02 ******************************************************************** Cp_EQUIL H_EQUIL S_EQUIL G_EQUIL V_EQUIL J.K-1 J J.K-1 J dm3 ******************************************************************** 2.46541E+02 8.09231E+05 8.54837E+02 -8.14958E+05 0.00000E+00
Mole fraction of system components: GAS C 7.4334E-02 O 4.9998E-01 Si 4.2569E-01
Data on 1 constituent marked with 'T' are extrapolated outside their valid temperature range **** End output table produced by TQCENL

Output of next call to TQCEL will be written to file 'result'

Output of subsequent calls to TQCEL/TQSHOW/TQMAPL will go to UNIT 6 again

C Get and display results
C First, we get the fraction of the first phase constituent of the first C phase. We get also its name with TQGNPC. Since the current amount C unit is 'mol', we get the 'mole fraction' (as opposed to 'mass C fraction') CALL TQGNSC(1, DSTR, NOERR) IF (NOERR .NE. 0) CALL CABORT(738,'TQGNSC',NOERR) CALL TQGETR('xp ', 1, 1, D1, NOERR) IF (NOERR .NE. 0) CALL CABORT(740,'TQGETR',NOERR) WRITE(UNIT=*,FMT=51) 'Mole fraction of system component '//DSTR WRITE(UNIT=*,FMT=57) 'in the GAS phase: ', D1

Output:
 Mole fraction of system component C                       
 in the GAS phase:  0.74334E-01

C Now get the equilibrium amount of the same constituent CALL TQGETR('a ',1,1,D1,NOERR) IF (NOERR .NE. 0) CALL CABORT(747,'TQGETR',NOERR) WRITE(UNIT=*,FMT=51) 'Equilibrium amount of '//DSTR WRITE(UNIT=*,FMT=57) 'in the GAS phase: ', D1

Output:
 Equilibrium amount of C                       
 in the GAS phase:   0.0000

C Note that the result is zero, because the gas phase is not C considered stable. The reason is that we are using the default C ambient pressure of 1 bar for the calculations, and from the C result table produced with the previous call to TQCEL or the C subsequent call to TQGETR we see that the activity of the gas C phase is less than unity, which means it is not considered C stable, and thus the equilibrium amounts of all its phase C constituents are zero.
C Using the option 'ac', TQGETR retrieves activities, in this case the C activity (fugacity) of the gas phase CALL TQGETR('ac ',1, 0, D1, NOERR) IF (NOERR .NE. 0) CALL CABORT(764,'TQGETR',NOERR) WRITE(UNIT=*,FMT=58) 'Activity of the GAS phase: ', D1
C Now TQGETR is used in a way that causes it to return an array. In C this case, we are going to retrieve an array that contains the C fugacities of all constituents of the gas phase.

C Get the fugacities: CALL TQGETR('ac ', 1, -1, DARR30, NOERR) IF (NOERR .NE. 0) CALL CABORT(774,'TQGETR',NOERR)
C Print them out, together with the names of the constituents C (note that these values correspond to the last column of the C result table previously output with TQCEL):
WRITE(UNIT=*,FMT=51) 'Fugacities of all gas phase constituents:'
DO I=1, NPCGAS CALL TQGNPC(1, I, DSTR, NOERR) IF (NOERR .NE. 0) CALL CABORT(784,'TQGNPC',NOERR)
WRITE(UNIT=*,FMT=57) DSTR // ': ', DARR30(I) ENDDO
WRITE(UNIT=*,FMT=51) ' '

Output:
 Activity of the GAS phase:  0.14450    


Fugacities of all gas phase constituents: C : 0.10955E-12 C2 : 0.19220E-15 C3 : 0.10532E-15 CO : 0.21481E-01 CO2 : 0.54503E-06 O : 0.42711E-11 O2 : 0.20213E-15 O3 : 0.96770E-31 Si : 0.12142E-04 Si2 : 0.10817E-06 Si2C : 0.20832E-06 Si3 : 0.23173E-07 SiC : 0.53182E-11 SiO : 0.12301 SiO2 : 0.39975E-07


C Get thermodynamic data of a phase constituent CALL TQGDPC('G ', 1, 1, D1, NOERR) IF (NOERR .NE. 0) CALL CABORT(795,'TQGDPC',NOERR) WRITE(UNIT=*,FMT=58) 'The (dimensionless) value of G for ' // * 'constituent 1 of phase 1 is: ', D1

Output:
 The (dimensionless) value of G for constituent 1 of phase 1 is:   23.826

C Note that when retrieving CP, H, S, or G using TQGDPC, the values are C returned 'dimensionless', which means they might have to be multiplied C by R*T. Since the default amount unit is mol, results are returned for C 1 mol. Note also that care has to be taken if a temperature unit C different from Kelvin has been used. Refer to the manual entry of C TQGDPC for a more detailed example.
C Check if we are working with the 'light' version. C If we do, omit the following target calculation(s). CALL TQLITE(ISLITE, NOERR) IF (ISLITE .EQ. 1) THEN WRITE(*,FMT='(3(1X,A,/))') * '*** Target calculations have been omitted here,', * '*** since they are not possible with the ', * '*** ''light'' version of ChemApp.'
ELSE
C Perform a target calculation. This example demonstrates how to C have ChemApp find the temperature at which a certain phase becomes C stable.
C We want to have ChemApp calculate the temperature at which the C phase 'SiO2(liquid)' becomes stable. First, determine the index C number of this phase (note that it is allowed to abbreviate its C name, as long as it is unambiguous): CALL TQINP('SiO2(liq', I, NOERR) IF (NOERR .NE. 0) CALL CABORT(827,'TQINP',NOERR)
C Next, define the formation target : CALL TQSETC('a ', I, 0, 0.D0, NUMCON, NOERR) IF (NOERR .NE. 0) CALL CABORT(831,'TQSETC',NOERR)
C Then call TQCEL to calculate the equilibrium, and pass it the C information about the target variable (we want ChemApp to vary C the _temperature_ until the phase is stable):
WRITE(UNIT=*,FMT=51) 'Performing a phase target calculation' WRITE(UNIT=*,FMT=51) 'Phase target: formation of SiO2(liquid)' WRITE(UNIT=*,FMT=53) 'Target variable: Temperature'
DARR2(1) = 2000.D0 WRITE(UNIT=*,FMT=51) '**** Begin output table produced by TQCEL' CALL TQCEL('t ',0,0,DARR2,NOERR) WRITE(UNIT=*,FMT=53) '**** End output table produced by TQCEL' IF (NOERR .NE. 0) CALL CABORT(845,'TQCEL',NOERR) DARR2(1) = 0.D0

Output:
 Performing a phase target calculation
 Phase target: formation of SiO2(liquid)
 Target variable: Temperature


**** Begin output table produced by TQCEL
*T = 1995.99 K P = 1.00000E+00 bar V = 0.00000E+00 dm3
STREAM CONSTITUENTS AMOUNT/mol CO/GAS/ 3.5701E-02 Si3/GAS/ 3.0000E+00
EQUIL AMOUNT MOLE FRACTION FUGACITY PHASE: GAS mol bar SiO 0.0000E+00 7.7442E-01 3.1442E-01 CO 0.0000E+00 2.2547E-01 9.1542E-02 Si 0.0000E+00 9.8748E-05 4.0093E-05 CO2 0.0000E+00 1.0313E-05 4.1873E-06 Si2C 0.0000E+00 2.5977E-06 1.0547E-06 Si2 0.0000E+00 1.1230E-06 4.5595E-07 SiO2 0.0000E+00 5.8253E-07 2.3651E-07 Si3 0.0000E+00 2.2928E-07 9.3089E-08 SiC 0.0000E+00 1.4260E-10 5.7898E-11 O 0.0000E+00 9.6352E-11 3.9120E-11 C 0.0000E+00 3.4767E-12 1.4116E-12 C2 0.0000E+00 1.2342E-14 5.0109E-15 O2 0.0000E+00 8.8186E-15 3.5804E-15 C3 0.0000E+00 8.5734E-15 3.4809E-15 O3 0.0000E+00 2.7760E-29 1.1271E-29 TOTAL: 0.0000E+00 1.0000E+00 4.0601E-01 mol ACTIVITY Si 8.9464E+00 1.0000E+00 SiC 3.5701E-02 1.0000E+00 SiO2(cristobalite) 1.7851E-02 1.0000E+00 SiO2(liquid) 0.0000E+00 1.0000E+00 SiO2(tridymite) T 0.0000E+00 9.9861E-01 SiO2(quartz) T 0.0000E+00 8.9882E-01 C 0.0000E+00 5.0396E-02 ******************************************************************** Cp_EQUIL H_EQUIL S_EQUIL G_EQUIL V_EQUIL J.K-1 J J.K-1 J dm3 ******************************************************************** 2.46559E+02 8.32897E+05 8.66988E+02 -8.97600E+05 0.00000E+00
Mole fraction of system components: GAS C 1.1274E-01 O 4.9998E-01 Si 3.8728E-01
Data on 2 constituents marked with 'T' are extrapolated outside their valid temperature range **** End output table produced by TQCEL

C Note that in the resulting ChemSage table, the temperature is marked C with an asterisk, meaning that it has been calculated by ChemApp, and C the activity of SiO2(liquid) is unity, whereas its equilibrium amount C is zero, meaning it just forms at this temperature
C Retrieve the calculated temperature: CALL TQGETR('t ', 0, 0, D1, NOERR) IF (NOERR .NE. 0) CALL CABORT(856,'TQGETR',NOERR) WRITE(UNIT=*,FMT=58) 'Calculated formation temperature of '// * 'SiO2(liquid): ', D1

Output:
 Calculated formation temperature of SiO2(liquid):   1996.0


C One-dimensional phase mapping calculations. The example below C demonstrates how to locate all phase transitions in a given C temperature range.
C One-dimensional phase mapping calculations are only possible with C version 3.x and later of ChemApp. WRITE(UNIT=*,FMT=51) 'Performing one-dimensional phase' // * ' mapping calculations'
C Initialise the variable to zero that holds the number of times TQMAP C found results RESNO = 0
C Remove all previous conditions CALL TQREMC(-2, NOERR) IF (NOERR .NE. 0) CALL CABORT(877,'TQREMC',NOERR)
C Determine the index number for the phase SiO2(quartz) CALL TQINP('SiO2(quartz) ', I, NOERR) IF (NOERR .NE. 0) CALL CABORT(881,'TQINP',NOERR)
C Enter one mol of SiO2 CALL TQSETC('IA ', I, 0, 1.0D0, NUMCON, NOERR) IF (NOERR .NE. 0) CALL CABORT(885,'TQSETC',NOERR)
C The temperature search interval is supposed to range from 300 to C 3000 K: DARR2(1) = 300.D0 DARR2(2) = 3000.D0
C First call to TQMAP, note the 'F' ('First') in the option parameter CALL TQMAP('TF ',0,0,DARR2,ICONT,NOERR) IF (NOERR .NE. 0) CALL CABORT(894,'TQMAP',NOERR)
C The variable RESNO keeps track of the number of times we called TQMAP: RESNO = RESNO + 1
C Retrieve and print the temperature, which we know is DARR2(1) CALL TQGETR('T ',0,0,D1,NOERR) IF (NOERR .NE. 0) CALL CABORT(901,'TQGETR',NOERR)
WRITE(*,59) '*** Lower interval boundary:', D1, ' K'
30 CONTINUE
C TQMAP is called again. Note the 'N' ('Next') in the option C parameter. If we are at the first phase boundary (RESNO is still 2), C we call TQMAPL for a change to produce a ChemSage output table...
IF (RESNO .EQ. 2) THEN WRITE(*,FMT=51) * '*** ChemSage result table for the first phase ' // * 'boundary found:' WRITE(UNIT=*,FMT=51) * '**** Begin output table produced by TQMAPL' CALL TQMAPL('TN ',0,0,DARR2,ICONT,NOERR) IF (NOERR .NE. 0) CALL CABORT(918,'TQMAPL',NOERR) WRITE(UNIT=*,FMT=51) * '**** End output table produced by TQMAPL' IF (NOERR .NE. 0) CALL CABORT(921,'TQMAPL',NOERR) WRITE(*,*) C ...otherwise we just call TQMAP: ELSE CALL TQMAP('TN ',0,0,DARR2,ICONT,NOERR) IF (NOERR .NE. 0) CALL CABORT(926,'TQMAP',NOERR) ENDIF RESNO = RESNO +1
C Get the temperature... CALL TQGETR('T ',0,0,D1,NOERR) IF (NOERR .NE. 0) CALL CABORT(932,'TQGETR',NOERR)
C ... and print the entry for the result table. If we have called TQMAP C twice already, we know that we found a phase boundary. If not, we have C retrieved the temperature value of the upper interval boundary C DARR2(2): IF (RESNO .GT. 2) THEN WRITE(*,59) '*** Phase boundary found at ', D1, ' K' ELSE WRITE(*,59) '*** Upper interval boundary: ', D1, ' K' ENDIF
C For as long as ICONT is positive, we need to make further calls to C TQMAP IF (ICONT .GT. 0) GOTO 30

Output:
 Performing one-dimensional phase mapping calculations
 *** Lower interval boundary:  300.00 K
 *** Upper interval boundary:  3000.00 K
 *** ChemSage result table for the first phase boundary found:
 **** Begin output table produced by TQMAPL

*T = 1140.10 K P = 1.00000E+00 bar V = 0.00000E+00 dm3
STREAM CONSTITUENTS AMOUNT/mol SiO2(quartz) 1.0000E+00
EQUIL AMOUNT MOLE FRACTION FUGACITY PHASE: GAS mol bar SiO 0.0000E+00 6.2502E-01 4.4829E-16 O2 0.0000E+00 2.5088E-01 1.7994E-16 O 0.0000E+00 1.2327E-01 8.8415E-17 SiO2 0.0000E+00 8.3208E-04 5.9681E-19 Si 0.0000E+00 4.4651E-15 3.2026E-30 O3 0.0000E+00 2.4636E-19 1.7670E-34 Si2 0.0000E+00 5.5659E-36 3.9921E-51 Si3 0.0000E+00 NOT CALCD. 6.4667E-69 TOTAL: 0.0000E+00 1.0000E+00 7.1725E-16 mol ACTIVITY SiO2(quartz) 1.0000E+00 1.0000E+00 SiO2(tridymite) 0.0000E+00 1.0000E+00 SiO2(cristobalite) 0.0000E+00 9.9398E-01 SiO2(liquid) 0.0000E+00 6.4439E-01 Si 0.0000E+00 2.5802E-17 ******************************************************************** Cp_EQUIL H_EQUIL S_EQUIL G_EQUIL V_EQUIL J.K-1 J J.K-1 J dm3 ******************************************************************** 7.04991E+01 -8.55398E+05 1.25413E+02 -9.98382E+05 0.00000E+00
Mole fraction of system components: GAS O 6.6667E-01 Si 3.3333E-01 **** End output table produced by TQMAPL
*** Phase boundary found at 1140.10 K *** Phase boundary found at 1738.28 K *** Phase boundary found at 1995.99 K

C With the above example the temperatures of all phase boundaries in a C system which contains 1 mol of SiO2 have been calculated. Thus the C phase boundaries determined reflect the stability ranges of the C various modifications of SiO2. Also note again that the first two C temperatures determined are _no_ phase boundaries, but the lower and C upper limit of the search interval (DARR2(1) and DARR2(2)).
ENDIF
C Since from now on we will be using streams instead of global C conditions, we first remove all conditions and targets set up to this C point. Note that using a '-1' instead of the '-2' in the call to C TQREMC would also reset ChemApp to default units and values. CALL TQREMC(-2, NOERR) IF (NOERR .NE. 0) CALL CABORT(963,'TQREMC',NOERR)
WRITE(UNIT=*,FMT=53) 'Calculations using streams:'
C Set name and temperature for a stream via the array TP, which C will be passed to TQSTTP when each of the streams is created. TP(1) = 1000.D0 TP(2) = 1.D0
C Create 3 streams: CALL TQSTTP('stream1 ', TP, NOERR) IF (NOERR .NE. 0) CALL CABORT(974,'TQSTTP',NOERR) CALL TQSTTP('stream2 ', TP, NOERR) IF (NOERR .NE. 0) CALL CABORT(976,'TQSTTP',NOERR) CALL TQSTTP('stream3 ', TP, NOERR) IF (NOERR .NE. 0) CALL CABORT(978,'TQSTTP',NOERR)
C Set constituent amounts for each stream CALL TQSTCA('stream1 ', 1, 4, 1.D0,NOERR) IF (NOERR .NE. 0) CALL CABORT(982,'TQSTCA',NOERR) CALL TQSTCA('stream2 ', 1, 12, 3.D0,NOERR) IF (NOERR .NE. 0) CALL CABORT(984,'TQSTCA',NOERR) CALL TQSTCA('stream3 ', 1, 8, 2.D0,NOERR) IF (NOERR .NE. 0) CALL CABORT(986,'TQSTCA',NOERR)
C Remove the last stream CALL TQSTRM('stream3 ', NOERR) IF (NOERR .NE. 0) CALL CABORT(990,'TQSTRM',NOERR)
C Set equilibrium temperature CALL TQSTEC('t ', 0, 1800.D0, NOERR) IF (NOERR .NE. 0) CALL CABORT(994,'TQSTEC',NOERR)
C Calculate and list equilibrium WRITE(UNIT=*,FMT=51) '**** Begin output table produced by TQCEL' CALL TQCEL(' ',0,0,DARR2,NOERR) WRITE(UNIT=*,FMT=53) '**** End output table produced by TQCEL' IF (NOERR .NE. 0) CALL CABORT(1000,'TQCEL',NOERR)

Output:
 Calculations using streams:


**** Begin output table produced by TQCEL
T = 1800.00 K P = 1.00000E+00 bar V = 0.00000E+00 dm3
STREAM CONSTITUENTS AMOUNT/mol TEMPERATURE/K PRESSURE/bar STREAM CO/GAS/ 1.0000E+00 1000.00 1.0000E+00 2 Si3/GAS/ 3.0000E+00 1000.00 1.0000E+00 3
EQUIL AMOUNT MOLE FRACTION FUGACITY PHASE: GAS mol bar SiO 0.0000E+00 9.1172E-01 4.1405E-02 CO 0.0000E+00 8.8214E-02 4.0062E-03 Si 0.0000E+00 6.7137E-05 3.0490E-06 CO2 0.0000E+00 1.1335E-06 5.1478E-08 Si2C 0.0000E+00 7.0003E-07 3.1792E-08 Si2 0.0000E+00 4.5031E-07 2.0451E-08 SiO2 0.0000E+00 1.1236E-07 5.1026E-09 Si3 0.0000E+00 1.0170E-07 4.6188E-09 SiC 0.0000E+00 7.3913E-12 3.3567E-13 O 0.0000E+00 7.2737E-12 3.3033E-13 C 0.0000E+00 1.2567E-13 5.7072E-15 O2 0.0000E+00 1.6031E-16 7.2803E-18 C2 0.0000E+00 9.7381E-17 4.4225E-18 C3 0.0000E+00 4.0495E-17 1.8391E-18 O3 0.0000E+00 8.7034E-33 3.9526E-34 TOTAL: 0.0000E+00 1.0000E+00 4.5415E-02 mol ACTIVITY Si 7.5000E+00 1.0000E+00 SiC 1.0000E+00 1.0000E+00 SiO2(cristobalite) 5.0000E-01 1.0000E+00 SiO2(tridymite) 0.0000E+00 9.9970E-01 SiO2(liquid) 0.0000E+00 9.3907E-01 SiO2(quartz) T 0.0000E+00 9.1491E-01 C 0.0000E+00 2.2526E-02 ******************************************************************** DELTA_Cp DELTA_H DELTA_S DELTA_G DELTA_V J.K-1 J J.K-1 J dm3 ******************************************************************** 7.76666E+01 -1.68039E+06 -3.77532E+02 -2.00297E+06 -3.32581E+02
Mole fraction of system components: GAS C 4.4110E-02 O 4.9998E-01 Si 4.5591E-01
Data on 1 constituent marked with 'T' are extrapolated outside their valid temperature range **** End output table produced by TQCEL


C Get thermodynamic properties of a stream (note that as of version C 3.2.0 of ChemApp, TQSTXP can also be called _before_ the equilibrium C calculation) CALL TQSTXP('stream1 ', 'H ', D1, NOERR) IF (NOERR .NE. 0) CALL CABORT(1008,'TQSTXP',NOERR)
WRITE(UNIT=*,FMT=58) 'Enthalpy of ''stream1'': ', D1

Output:
 Enthalpy of 'stream1':  -88842.



WRITE(UNIT=*, FMT=53) 'Demonstrating the use of ChemApp ' // * 'subroutines involving sublattices.' C To demonstrate the subroutines dealing with sublattices, a different C data-file (subl-ex.dat) needs to be loaded. This data-file contains an C extract of the system Co-Cr-Fe: the phase SIGMA:30, which is modelled C according to the sublattice formalism, and the BCC phase, described by C a Redlich-Kister polynomial. Both phases are each included twice, to C account for miscibility gaps. CALL TQOPNA('subl-ex.dat', UNITNO, NOERR) IF (NOERR .NE. 0) CALL CABORT(1024,'TQOPNA',NOERR)
C Read data-file CALL TQRFIL(NOERR) IF (NOERR .NE. 0) CALL CABORT(1028,'TQRFIL',NOERR)
C Close data-file CALL TQCLOS(UNITNO, NOERR) IF (NOERR .NE. 0) CALL CABORT(1032,'TQCLOS',NOERR)
C The first of the two identical copies of the SIGMA:30 phase, which C ChemApp calls SIGMA:30#1, will be investigated with respect to the C number of sublattices, the number of sublattice constituents on each C sublattice, and the names of the sublattice constituents.
C Get the index number for the phase SIGMA:30#1 PNAME = 'SIGMA:30#1' CALL TQINP(PNAME, INDEXP, NOERR) IF (NOERR .NE. 0) CALL CABORT(1042,'TQINP',NOERR)
C Get the number of sublattices CALL TQNOSL(INDEXP, NSL, NOERR) IF (NOERR .NE. 0) CALL CABORT(1046,'TQNOSL',NOERR) WRITE(UNIT=*,FMT=60) PNAME // ' has ', NSL, ' sublattices'

Output:
 Demonstrating the use of ChemApp subroutines involving sublattices.


SIGMA:30#1 has 3 sublattices

C Loop over all sublattices DO INDEXL=1, NSL
C Get the number of sublattice constituents CALL TQNOLC(INDEXP, INDEXL, NSLCON, NOERR) IF (NOERR .NE. 0) CALL CABORT(1055,'TQNOLC',NOERR) WRITE(UNIT=*,FMT=56) 'Sublattice ', INDEXL, ' has ', NSLCON, * ' constituents with the following names:'
C Get the name of each sublattice constituent DO INDEXC=1, NSLCON CALL TQGNLC(INDEXP, INDEXL ,INDEXC, CNAME, NOERR) IF (NOERR .NE. 0) CALL CABORT(1062,'TQGNLC',NOERR)
C The reverse (getting the index number for the name of the sublattice C constituent just retrieved), is rather superfluous here, and only used C to demonstrate the call to TQINLC: CALL TQINLC(CNAME, INDEXP, INDEXL, I, NOERR) IF (NOERR .NE. 0) CALL CABORT(1068,'TQINLC',NOERR) WRITE(UNIT=*,FMT=60) ' ',I,':' // CNAME ENDDO ENDDO

Output:
 Sublattice   1 has   2 constituents with the following names:
     1:Co                      
     2:Fe                      
 Sublattice   2 has   1 constituents with the following names:
     1:Cr                      
 Sublattice   3 has   3 constituents with the following names:
     1:Co                      
     2:Cr                      
     3:Fe


C Set the temperature to 1000 K CALL TQSETC('T ', 0, 0, 1000.D0, NUMCON, NOERR) IF (NOERR .NE. 0) CALL CABORT(1077,'TQSETC',NOERR)
C Set the incoming amounts to 0.25 mol Co, 0.25 mol Cr, and 0.5 mol Fe CALL TQINSC('Co ', INDEXC, NOERR) IF (NOERR .NE. 0) CALL CABORT(1081,'TQINSC',NOERR) CALL TQSETC('IA ', 0, INDEXC, 0.25D0, NUMCON, NOERR) IF (NOERR .NE. 0) CALL CABORT(1083,'TQSETC',NOERR)
CALL TQINSC('Cr ', INDEXC, NOERR) IF (NOERR .NE. 0) CALL CABORT(1086,'TQINSC',NOERR) CALL TQSETC('IA ', 0, INDEXC, 0.25D0, NUMCON, NOERR) IF (NOERR .NE. 0) CALL CABORT(1088,'TQSETC',NOERR)
CALL TQINSC('Fe ', INDEXC, NOERR) IF (NOERR .NE. 0) CALL CABORT(1091,'TQINSC',NOERR) CALL TQSETC('IA ', 0, INDEXC, 0.50D0, NUMCON, NOERR) IF (NOERR .NE. 0) CALL CABORT(1093,'TQSETC',NOERR)

WRITE(UNIT=*, FMT=55) 'Calculate the equilibrium, ' // * 'get information on the stable sublattice phases.'
C Calculate the equilibrium DARR2(1) = 0.D0 CALL TQCE(' ', 0, 0, DARR2, NOERR) IF (NOERR .NE. 0) CALL CABORT(1102,'TQCE',NOERR)
C Get information on the stable sublattice phases.
C Loop over all phases, and if a phase is stable and a sublattice phase, C print information about the sublattices CALL TQNOP(NPHASE, NOERR) IF (NOERR .NE. 0) CALL CABORT(1109,'TQNOP',NOERR)
DO INDEXP=1, NPHASE
C Check if the phase is stable, otherwise we don't need to consider it CALL TQGETR('A ', INDEXP, 0, D1, NOERR) IF (NOERR .NE. 0) CALL CABORT(1115,'TQGETR',NOERR) IF (D1 .GT. 0.D0) THEN
C Get the name and model of the phase CALL TQGNP(INDEXP, PNAME, NOERR) IF (NOERR .NE. 0) CALL CABORT(1120,'TQGNP',NOERR) CALL TQMODL(INDEXP, MNAME, NOERR) IF (NOERR .NE. 0) CALL CABORT(1122,'TQMODL',NOERR)
C If the model name of the phase starts with 'SUB' we have a sublattice C phase IF (INDEX(MNAME,'SUB') .EQ. 1) THEN
C Print a header similar to the one in the ChemSage output table WRITE(UNIT=*,FMT=51) 'Mole fraction of the ' // * 'sublattice constituents in ' // PNAME
C Get the number of sublattices CALL TQNOSL(INDEXP, NSL, NOERR) IF (NOERR .NE. 0) CALL CABORT(1134,'TQNOSL',NOERR)
C Loop over all sublattices DO INDEXL=1, NSL
write(UNIT=*,FMT=60) 'Sublattice ',INDEXL,' :' C Get the number of sublattice constituents CALL TQNOLC(INDEXP, INDEXL, NSLCON, NOERR) IF (NOERR .NE. 0) CALL CABORT(1142,'TQNOLC', * NOERR)
DO INDEXC=1, NSLCON C Get the name of each sublattice constituent... CALL TQGNLC(INDEXP, INDEXL, INDEXC, CNAME, NOERR) IF (NOERR .NE. 0) CALL CABORT(1148,'TQGNLC', * NOERR)
C ... and its mole fraction in the sublattice CALL TQGTLC(INDEXP, INDEXL, INDEXC, D1, NOERR) IF (NOERR .NE. 0) CALL CABORT(1153,'TQGTLC', * NOERR)
WRITE(UNIT=*,FMT=57) CNAME, D1
ENDDO ENDDO ENDIF ENDIF ENDDO WRITE(UNIT=*,FMT=51) ' '

Output:

Calculate the equilibrium, get information on the stable sublattice phases.

Mole fraction of the sublattice constituents in SIGMA:30#1 Sublattice 1 : Co 0.54238 Fe 0.45762 Sublattice 2 : Cr 1.0000 Sublattice 3 : Co 0.10926 Cr 0.69192 Fe 0.19882


C The following section demonstrates the use of subroutines that were C introduced in ChemApp V4.0.0 to support 'transparent data-files'.
C Retrieve the licensee's user ID CALL TQGTID(ID, NOERR) IF (NOERR .NE. 0) CALL CABORT(1172,'TQGTID',NOERR)
C Retrieve the licensee's name CALL TQGTNM(NAME, NOERR) IF (NOERR .NE. 0) CALL CABORT(1176,'TQGTNM',NOERR)
C Retrieve the program ID CALL TQGTPI(PID, NOERR) IF (NOERR .NE. 0) CALL CABORT(1180,'TQGTPI',NOERR)
C Print all three WRITE(UNIT=*,FMT=51) 'Licensee''s user ID: ' // * ID(1:STRLEN(ID)) WRITE(UNIT=*,FMT=51) 'Licensee''s name : ' // * NAME(1:STRLEN(NAME)) WRITE(UNIT=*,FMT=53) 'Program ID : ' // * PID(1:STRLEN(PID))

Output:
 Licensee's user ID: 5001
 Licensee's name   : GTT - Technologies
 Program ID        : CAFU

C The following pieces of information are only meaningful if a version C of ChemApp is used that requires a dongle (hardware key). C Get the HASP dongle type and id CALL TQGTHI(HASPT, HASPID, NOERR) IF (NOERR .NE. 0) CALL CABORT(1195,'TQGTHI',NOERR)
C Get the ChemApp license expiration date (month and year) CALL TQGTED(EDMON, EDYEAR, NOERR) IF (NOERR .NE. 0) CALL CABORT(1199,'TQGTED',NOERR)
C Print info if HASP dongle is used: IF (HASPID .NE. 0) THEN WRITE(UNIT=*,FMT=51) 'HASP dongle type : ' // * HASPT(1:STRLEN(HASPT)) WRITE(UNIT=*,FMT=*) 'HASP dongle id : ', HASPID
WRITE(UNIT=*,FMT=62) 'ChemApp license expiration date ' // * '(month/year): ', EDMON, '/', EDYEAR ELSE WRITE(UNIT=*,FMT=53) 'This ChemApp version does not ' // * 'require a HASP hardware key (dongle)' ENDIF

Output:
 This ChemApp version does not require a HASP hardware key (dongle)

C Get the default unit number associated with error output CALL TQGIO('ERROR ', EUNIT, NOERR) IF (NOERR .NE. 0) CALL CABORT(1217,'TQGIO',NOERR)
C Turn off automatic error reporting for the next call to TQCIO. CALL TQCIO('ERROR ', 0, NOERR) IF (NOERR .NE. 0) CALL CABORT(1221,'TQCIO',NOERR)
C The following section is only executed if the file "cosiex.cst", C a sample thermochemical data-file in transparent format, is C present. It is not included in regular distributions of ChemApp CALL TQOPNT('cosiex.cst', UNITNO, I)
C Turn automatic error reporting back on. CALL TQCIO('ERROR ', EUNIT, NOERR) IF (NOERR .NE. 0) CALL CABORT(1230,'TQCIO',NOERR)
IF (I .NE. 0) THEN WRITE(UNIT=*,FMT=*) 'Skipping transparent file ' // * 'subroutines, since file "cosiex.cst" cannot be read' ELSE
C Read data-file CALL TQRCST(NOERR) IF (NOERR .NE. 0) CALL CABORT(1239,'TQRCST',NOERR)
C Close data-file CALL TQCLOS(UNITNO, NOERR) IF (NOERR .NE. 0) CALL CABORT(1243,'TQCLOS',NOERR)
C Once the transparent data-file has been read, information on its C header can be retrieved. CALL TQGTRH(TFHVER, TFHNWP, TFHVNW, TFHNRP, TFHVNR, TFHDTC, * TFHDTE, TFHID, TFHUSR, TFHREM, NOERR) IF (NOERR .NE. 0) CALL CABORT(1249,'TQGTRH',NOERR)
WRITE(UNIT=*, FMT=*) 'Version number of the transparent ' // * 'file header format:', TFHVER WRITE(UNIT=*, FMT=*) 'Name of the program which wrote the ' // * 'data-file: ', TFHNWP(1:STRLEN(TFHNWP)) WRITE(UNIT=*, FMT=*) 'Version number of the writing program: ', * TFHVNW(1), '.', TFHVNW(2), '.', TFHVNW(3) WRITE(UNIT=*, FMT=*) 'Programs which are permitted to read ' // * 'the data-file: ', TFHNRP(1:STRLEN(TFHNRP)) WRITE(UNIT=*, FMT=*) 'Min. version number of the reading ' // * 'program: ', TFHVNR(1), '.', TFHVNR(2), '.', TFHVNR(3) WRITE(UNIT=*, FMT=61) 'File was created on ', * TFHDTC(1),'/', TFHDTC(2),'/',TFHDTC(3),' ', * TFHDTC(4),':',TFHDTC(5),':',TFHDTC(6) WRITE(UNIT=*, FMT=61) 'File will expire on ', * TFHDTE(1),'/', TFHDTE(2),'/',TFHDTE(3),' ', * TFHDTE(4),':',TFHDTE(5),':',TFHDTE(6) WRITE(UNIT=*, FMT=*) 'Licensee''s user ID(s): ', * TFHID(1:STRLEN(TFHID)) WRITE(UNIT=*, FMT=*) 'Licensee''s name: ', * TFHUSR(1:STRLEN(TFHUSR)) WRITE(UNIT=*, FMT=*) 'Remarks: ', * TFHREM(1:STRLEN(TFHREM))

Output:
 Version number of the transparent file header format: 1
 Name of the program which wrote the data-file: FactSage
 Version number of the writing program:  5. 0. 0
 Programs which are permitted to read the data-file: CAFU,CALI
 Min. version number of the reading program:  -1. -1. -1
 File was created on 2003/06/05 12:11:39
 File will expire on 2036/12/31 12:00:00
 Licensee's user ID(s): ????
 Licensee's name: GTT - Technologies
 Remarks: For use with ChemApp example programs

ENDIF

WRITE(UNIT=*,FMT=55) 'End of output from cademo1.'

Output:

End of output from cademo1.

END


ChemApp Programmer's Manual, Edition 3.6© GTT-Technologies, 2003