|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？注册
x
今天从网上搜到了一段chemkin的源程序，贴出来给大家看看，我觉得主要是读取链接文件，以及计算反应速率的程序，请高手鉴别一下。我想知道这有没有包括化学平衡反应的计算
C----------------------------------------------------------------------C
C C
SUBROUTINE CKLEN (LINC, LOUT, LI, LR, LC, IFLAG)
C
C START PROLOGUE
C
C SUBROUTINE CKLEN (LINC, LOUT, LENI, LENR, LENC, IFLAG)
C Returns the lengths required for the work arrays.
C
C INPUT
C
C LINC - Logical file number for the binary file.
C Data type - integer scalar
C LOUT - Output file for printed diagnostics.
C Data type - integer scalar
C
C OUTPUT
C LENI - Minimum length required for the integer work array.
C Data type - integer scalar
C LENR - Minimum length required for the real work array.
C Data type - integer scalar
C LENC - Minimum length required for the character work array.
C IFLAG - IFLAG=0 indicates successful reading of
C binary linking file; IFLAG>0 indicates
C error type.
C Data type - integer
C
C END PROLOGUE
C
IMPLICIT NONE
C
INTEGER IFLAG,II,KK,LC,LENC,LENCCK,LENI,LENICK,LENR,LENRCK,LI,
1 LINC,LOUT,LR,MAXORD,MAXSP,MAXTB,MAXTP,MM,N,NCHRG,NFL,
2 NIFAR,NIPAR,NITAR,NLIST,NLT,NOR,NRL,NRV,NSTO,NTB,NTHCF,
3 NW
DOUBLE PRECISION
1 CKMN
PARAMETER (NLIST = 5)
LOGICAL KERR, VOK, POK
CHARACTER LIST(NLIST)*16, PREC*16, VERS*16
COMMON /CKCONS/ PREC, VERS, KERR, LENI, LENR, LENC
C DATA LIST/';1.9';,';2.0';,';2.1';,';2.2';,';2.3';,';2.4';,';2.5';,';2.6';,
C 1 ';2.7';,';2.8';,';2.9';,';3.0';,';3.1';,';3.2';,';3.3';/
C DATA LIST /';3.4';,';3.5';,';3.6';/
DATA LIST /';3.6b';,';3.6c';,';3.6d';,';3.8';, ';3.9';/
C
VERS = '; ';
PREC = '; ';
LENI = 0
LENR = 0
LENC = 0
C
KERR = .FALSE.
IFLAG = 0
REWIND LINC
READ (LINC, ERR=999) VERS, PREC, KERR
C
VOK = .FALSE.
DO 5 N = 1, NLIST
IF (VERS .EQ. LIST(N)) VOK = .TRUE.
5 CONTINUE
C
POK = .FALSE.
C*****precision > double
IF (INDEX(PREC, ';DOUB';) .GT. 0) POK = .TRUE.
C*****END precision > double
C*****precision > single
C IF (INDEX(PREC, ';SING';) .GT. 0) POK = .TRUE.
C*****END precision > single
C
IF (KERR .OR. (.NOT.POK) .OR. (.NOT.VOK)) THEN
IF (KERR) THEN
WRITE (LOUT,';(/A,/A)';)
1 '; There is an error in the Chemkin binary file...';,
2 '; Check CHEMKIN INTERPRETER output for error conditions.';
ENDIF
IF (.NOT. VOK) THEN
WRITE (LOUT,';(/A,A)';)
1 '; Chemkin binary file is incompatible with Chemkin';,
2 '; Library Version 4.5';
ENDIF
IF (.NOT. POK) THEN
WRITE (LOUT,';(/A,A)';)
1 '; Precision of Chemkin binary file does not agree with';,
2 '; precision of Chemkin library';
ENDIF
IFLAG = 20
RETURN
ENDIF
C
READ (LINC, ERR=1111) LENICK, LENRCK, LENCCK, MM, KK, II,
1 MAXSP, MAXTB, MAXTP, NTHCF, NIPAR, NITAR,
2 NIFAR, NRV, NFL, NTB, NLT, NRL, NW, NCHRG,
3 NSTO, NOR, MAXORD, CKMN
REWIND LINC
C
LENI = LENICK
LENR = LENRCK
LENC = LENCCK
LI = LENI
LR = LENR
LC = LENC
RETURN
C
999 CONTINUE
WRITE (LOUT, 50)
50 FORMAT ('; Error reading Chemkin binary file.';)
IFLAG = 1
RETURN
1111 CONTINUE
WRITE (LOUT, 50)
IFLAG = 2
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKINIT (LENIWK, LENRWK, LENCWK, LINC, LOUT, ICKWRK,
1 RCKWRK, CCKWRK, IFLAG)
C
C START PROLOGUE
C
C SUBROUTINE CKINIT (LENIWK, LENRWK, LENCWK, LINC, LOUT, ICKWRK,
C RCKWRK, CCKWRK, IFLAG)*
C Reads the binary file and creates the internal work arrays
C ICKWRK, CCKWRK, and RCKWORK. CKINIT must be called before any
C other CHEMKIN subroutine is called. The work arrays must then
C be made available as input to the other CHEMKIN subroutines.
C
C INPUT
C LENIWK - Length of the integer work array, ICKWRK.
C Data type - integer scalar
C LENCWK - Length of the character work array, CCKWRK.
C The minimum length of CCKWRK(*) is MM + KK.
C Data type - integer scalar
C LENRWK - Length of the real work array, WORK.
C Data type - integer scalar
C LINC - Logical file number for the binary file.
C Data type - integer scalar
C LOUT - Output file for printed diagnostics.
C Data type - integer scalar
C
C OUTPUT
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C CCKWRK - Array of character work space.
C Data type - CHARACTER*16 array
C Dimension CCKWRK(*) at least LENCWK.
C IFLAG - IFLAG=0 indicates successful reading of
C binary linking file; IFLAG>0 indicates
C error type.
C Data type - integer
C END PROLOGUE
C
IMPLICIT NONE
C
INCLUDE ';ckstrt.h';
C
INTEGER I,ICKWRK,IFLAG,II,ITOC,ITOT,K,KK,L,LC,LENC,LENCWK,
1 LENI,LENIWK,LENR,LENRWK,LI,LINC,LOUT,LR,M,MAXORD,MAXSP,
2 MAXTB,MAXTP,MM,MXOR,N,NCHRG,NEI,NEX,NF1,NFL,NIF1,NIFAR,
3 NIJAN,NIPAR,NITAR,NJA,NLT,NOR,NRL,NRV,NSTO,NTB,NTHCF,
4 NTOT,NW
DOUBLE PRECISION
1 BIG,CKMIN,CKMN,EXPARG,PA,RCKWRK,RU,RUC,SMALL
DIMENSION ICKWRK(*), RCKWRK(*)
CHARACTER*(*) CCKWRK(*)
CHARACTER*16 VERS, PREC
LOGICAL IOK, ROK, COK, KERR
COMMON /CKCONS/ PREC, VERS, KERR, LENI, LENR, LENC
COMMON /MACH/ SMALL,BIG,EXPARG
COMMON /CMIN/ CKMIN
C
C Data about the machine dependent constants is carried in
C
C COMMON/MACH/SMALL,BIG,EXPARG
C
C Set the gas constant to the 1986 CODATA recommended
C value - Gas constant in cals/mole is determined by
C division of the later value by 4.184.
C
C The number for conversion to one atmosphere is exact.
C
PARAMETER (RU=8.314510D+07, RUC=RU/4.184D+07, PA=1.01325D+06)
C*****exponent range > +/-300
SMALL = 10.0D00**(-300)
BIG = 10.0D00**(+300)
C*****END exponent range > +/-300
EXPARG = DLOG(BIG)
C
WRITE (LOUT,15)
15 FORMAT (/1X,'; CKLIB: Chemical Kinetics Library';,
1 /1X,'; CHEMKIN-II Version 4.5, January 1995';,
2 /1X,'; DOUBLE PRECISION';)
C
CALL CKLEN (LINC, LOUT, LI, LR, LC, IFLAG)
C
IF (IFLAG .GT. 0) RETURN
IOK = (LENIWK .GE. LI)
ROK = (LENRWK .GE. LR)
COK = (LENCWK .GE. LC)
IF (.NOT.IOK .OR. .NOT.ROK .OR. .NOT.COK) THEN
IF (.NOT. IOK) WRITE (LOUT, 300) LI
IF (.NOT. ROK) WRITE (LOUT, 350) LR
IF (.NOT. COK) WRITE (LOUT, 375) LC
IFLAG = 20
RETURN
ENDIF
C
REWIND LINC
READ (LINC, ERR=110) VERS, PREC, KERR
READ (LINC, ERR=110) LENI, LENR, LENC, MM, KK, II,
1 MAXSP, MAXTB, MAXTP, NTHCF, NIPAR, NITAR,
2 NIFAR, NRV, NFL, NTB, NLT, NRL, NW, NCHRG,
3 NEI, NJA, NIJAN, NF1, NIF1, NEX,
4 NSTO, NOR, MAXORD, CKMN
C
IF (LEN(CCKWRK(1)) .LT. 16) THEN
WRITE (LOUT,475)
IFLAG = 21
RETURN
ENDIF
C
NMM = MM
NKK = KK
NII = II
MXSP = MAXSP
MXTB = MAXTB
MXTP = MAXTP
MXOR = MAXORD
NCP = NTHCF
NCP1 = NTHCF+1
NCP2 = NTHCF+2
NCP2T = NCP2*(MAXTP-1)
NPAR = NIPAR
NLAR = NITAR
NFAR = NIFAR
NTHB = NTB
NLAN = NLT
NFAL = NFL
NREV = NRV
NRLT = NRL
NWL = NW
NEIM = NEI
NJAR = NJA
NJAN = NIJAN
NFT1 = NIF1
NF1R = NF1
NEXC = NEX
NRNU= NSTO
NORD = NOR
MXORD= MAXORD
CKMIN= CKMN
C
C APPORTION work arrays
C
C SET ICKWRK(*)=1 TO FLAG THAT CKINIT HAS BEEN CALLED
C
ICKWRK(1) = 1
C
C STARTING LOCATIONS OF INTEGER SPACE
C
C! elemental composition of species
IcNC = 2
C! species phase array
IcPH = IcNC + KK*MM
C! species charge array
IcCH = IcPH + KK
C! # of temperatures for fit
IcNT = IcCH + KK
C! stoichiometric coefficients
IcNU = IcNT + KK
C! species numbers for the coefficients
IcNK = IcNU + MAXSP*II
C! # of non-zero coefficients (<0=reversible, >0=irreversible)
IcNS = IcNK + MAXSP*II
C! # of reactants
IcNR = IcNS + II
C! Landau-Teller reaction numbers
IcLT = IcNR + II
C! Reverse Landau-Teller reactions
IcRL = IcLT + NLAN
C! Fall-off reaction numbers
IcFL = IcRL + NRLT
C! Fall-off option numbers
IcFO = IcFL + NFAL
C! Fall-off enhanced species
IcKF = IcFO + NFAL
C! Third-body reaction numbers
IcTB = IcKF + NFAL
C! number of 3rd bodies for above
IcKN = IcTB + NTHB
C! array of species #';s for above
IcKT = IcKN + NTHB
C! Reverse parameter reaction numbers
IcRV = IcKT + MAXTB*NTHB
C! Radiation wavelength reactions
IcWL = IcRV + NREV
C! Electon-impact reaction numbers
IcEI = IcWL + NWL
C! Electron-impact temperature dependence flags
IcTD = IcEI + NEIM
C! Janev-Langer_Evans& ost type reaction numbers
IcJN = IcTD + NEIM
C! Reaction numbers using fit#1
IcF1 = IcJN + NJAN
C! Reaction numbers for excitation-only reactions
IcEX = IcF1 + NFT1
C! Real stoichometry reactions
IcRNU= IcEX + NEXC
C! Change of order reactions
IcORD= IcRNU + NRNU
C! Species for which there is a change of order
IcKOR= IcORD + NORD
C
ITOT = IcKOR + NORD*MXORD - 1
C
C STARTING LOCATIONS OF CHARACTER SPACE
C
C! start of element names
IcMM = 1
C! start of species names
IcKK = IcMM + MM
ITOC = IcKK + KK - 1
C
C STARTING LOCATIONS OF REAL SPACE
C
C! atomic weights
NcAW = 1
C! molecular weights
NcWT = NcAW + MM
C! temperature fit array for species
NcTT = NcWT + KK
C! thermodynamic coefficients
NcAA = NcTT + MAXTP*KK
C! Arrhenius coefficients (3)
NcCO = NcAA + (MAXTP-1)*NCP2*KK
C! Reverse coefficients
NcRV = NcCO + (NPAR+1)*II
C! Landau-Teller #';s for NLT reactions
NcLT = NcRV + (NPAR+1)*NREV
C! Reverse Landau-Teller #';s
NcRL = NcLT + NLAR*NLAN
C! Fall-off parameters for NFL reactions
NcFL = NcRL + NLAR*NRLT
C! 3rd body coef';nts for NTHB reactions
NcKT = NcFL + NFAR*NFAL
C! wavelength
NcWL = NcKT + MAXTB*NTHB
C! Janev-type coefficients
NcJN = NcWL + NWL
C! Fit#1 parameters
NcF1 = NcJN + NJAR*NJAN
C! Excitation-only reaction energy loss
NcEX = NcF1 + NF1R*NFT1
C! real stoichometric coefficients
NcRNU= NcEX + NEXC
C! change of order for species/reactions
NcKOR= NcRNU + NRNU*MXSP
C! universal gas constant
NcRU = NcKOR + NORD*MXORD
C! universal gas constant in units
NcRC = NcRU + 1
C! pressure of one atmosphere
NcPA = NcRC + 1
C! intermediate temperature-dependent forward rates
NcKF = NcPA + 1
C! intermediate temperature-dependent reverse rates
NcKR = NcKF + II
C! internal work space of length kk
NcK1 = NcKR + II
C! ';ditto';
NcK2 = NcK1 + KK
C! ';ditto';
NcK3 = NcK2 + KK
C! ';ditto';
NcK4 = NcK3 + KK
NcI1 = NcK4 + KK
NcI2 = NcI1 + II
NcI3 = NcI2 + II
NcI4 = NcI3 + II
NTOT = NcI4 + II - 1
C
C SET UNIVERSAL CONSTANTS IN CGS UNITS
C
RCKWRK(NcRU) = RU
RCKWRK(NcRC) = RUC
RCKWRK(NcPA) = PA
C
C!element names, !atomic weights
READ (LINC,err=111) (CCKWRK(IcMM+M-1), RCKWRK(NcAW+M-1), M=1,MM)
C
C!species names, !composition, !phase, !charge, !molec weight,
C!# of fit temps, !array of temps, !fit coeff';nts
READ (LINC,err=222) (CCKWRK(IcKK+K-1),
1 (ICKWRK(IcNC+(K-1)*MM + M-1),M=1,MM),
2 ICKWRK(IcPH+K-1),
3 ICKWRK(IcCH+K-1),
4 RCKWRK(NcWT+K-1),
5 ICKWRK(IcNT+K-1),
6 (RCKWRK(NcTT+(K-1)*MAXTP + L-1),L=1,MAXTP),
7 ((RCKWRK(NcAA+(L-1)*NCP2+(K-1)*NCP2T+N-1),
8 N=1,NCP2), L=1,(MAXTP-1)), K = 1,KK)
C
IF (II .EQ. 0) RETURN
C
C!# spec,reactants, !Arr. coefficients, !stoic coef, !species numbers
READ (LINC,end=100,err=333)
1 (ICKWRK(IcNS+I-1), ICKWRK(IcNR+I-1),
2 (RCKWRK(NcCO+(I-1)*(NPAR+1)+N-1), N=1,NPAR),
3 (ICKWRK(IcNU+(I-1)*MAXSP+N-1),
4 ICKWRK(IcNK+(I-1)*MAXSP+N-1), N=1,MAXSP),
5 I = 1,II)
C
C PERTURBATION FACTOR
C
DO 10 I = 1, II
RCKWRK(NcCO + (I-1)*(NPAR+1) + NPAR) = 1.D00
10 CONTINUE
C
IF (NREV .GT. 0) READ (LINC,err=444)
1 (ICKWRK(IcRV+N-1), (RCKWRK(NcRV+(N-1)*(NPAR+1)+L-1),
1 L=1,NPAR), N = 1,NREV)
C
IF (NFAL .GT. 0) READ (LINC,err=555)
1 (ICKWRK(IcFL+N-1), ICKWRK(IcFO+N-1), ICKWRK(IcKF+N-1),
2 (RCKWRK(NcFL+(N-1)*NFAR+L-1),L=1,NFAR),N=1,NFAL)
C
IF (NTHB .GT. 0) READ (LINC,err=666)
1 (ICKWRK(IcTB+N-1), ICKWRK(IcKN+N-1),
2 (ICKWRK(IcKT+(N-1)*MAXTB+L-1),
3 RCKWRK(NcKT+(N-1)*MAXTB+L-1),L=1,MAXTB),N=1,NTHB)
C
IF (NLAN .GT. 0) READ (LINC,err=777)
1 (ICKWRK(IcLT+N-1), (RCKWRK(NcLT+(N-1)*NLAR+L-1),L=1,NLAR),
2 N=1,NLAN)
C
IF (NRLT .GT. 0) READ (LINC,err=888)
1 (ICKWRK(IcRL+N-1), (RCKWRK(NcRL+(N-1)*NLAR+L-1),L=1,NLAR),
2 N=1,NRLT)
C
IF (NWL .GT. 0) READ (LINC,err=999)
1 (ICKWRK(IcWL+N-1), RCKWRK(NcWL+N-1), N=1,NWL)
C
IF (NEIM .GT. 0) READ (LINC,err=1001)
1 (ICKWRK(IcEI+N-1), ICKWRK(IcTD+N-1), N=1,NEIM)
IF (NJAN .GT. 0) READ (LINC,err=1002)
1 (ICKWRK(IcJN+N-1), (RCKWRK(NcJN+(N-1)*NJAR+L-1),L=1,NJAR),
2 N=1,NJAN)
IF (NFT1 .GT. 0) READ (LINC,err=1003)
1 (ICKWRK(IcF1+N-1), (RCKWRK(NcF1+(N-1)*NF1R+L-1),L=1,NF1R),
2 N=1,NFT1)
IF (NEXC .GT. 0) READ (LINC,err=1004)
1 (ICKWRK(IcEX+N-1), RCKWRK(NcEX+N-1), N=1,NEXC)
C
IF (NRNU .GT. 0) READ (LINC,err=1111)
1 (ICKWRK(IcRNU+N-1), (RCKWRK(NcRNU+(N-1)*MAXSP+L-1),L=1,MAXSP),
2 N=1,NRNU)
C
IF (NORD .GT. 0) READ (LINC,err=2222)
1 (ICKWRK(IcORD+N-1), (ICKWRK(IcKOR+(N-1)*MXORD+L-1),
2 RCKWRK(NcKOR+(N-1)*MXORD+L-1),
3 L = 1, MXORD), N=1, NORD)
C
100 CONTINUE
RETURN
C
110 WRITE (LOUT,*) '; Error reading binary file...';
IFLAG = 1
RETURN
111 WRITE (LOUT,*) '; Error reading element data...';
IFLAG = 2
RETURN
222 WRITE (LOUT,*) '; Error reading species data...';
IFLAG = 3
RETURN
333 WRITE (LOUT,*) '; Error reading reaction data...';
IFLAG = 4
RETURN
444 WRITE (LOUT,*) '; Error reading reverse Arrhenius parameters...';
IFLAG = 5
RETURN
555 WRITE (LOUT,*) '; Error reading Fall-off data...';
IFLAG = 6
RETURN
666 WRITE (LOUT,*) '; Error reading third-body data...';
IFLAG = 7
RETURN
777 WRITE (LOUT,*) '; Error reading Landau-Teller data...';
IFLAG = 8
RETURN
888 WRITE (LOUT,*) '; Error reading reverse Landau-Teller data...';
IFLAG = 9
RETURN
999 WRITE (LOUT,*) '; Error reading Wavelength data...';
IFLAG = 10
RETURN
1001 WRITE (LOUT,*) '; Error reading Electron-impact data...';
IFLAG = 11
RETURN
1002 WRITE (LOUT,*) '; Error reading Janev-Langer-Evans-Post data...';
IFLAG = 12
RETURN
1003 WRITE (LOUT,*) '; Error reading Fit # 1 data...';
IFLAG = 13
RETURN
1004 WRITE (LOUT,*) '; Error reading excitation-reaction data...';
IFLAG = 14
RETURN
1111 WRITE (LOUT,*) '; Error reading real stoichometric data...';
IFLAG = 15
RETURN
2222 WRITE (LOUT,*) '; Error reading order data...';
IFLAG = 16
RETURN
C
300 FORMAT (10X,';ICKWRK MUST BE DIMENSIONED AT LEAST ';,I5)
350 FORMAT (10X,';RCKWRK MUST BE DIMENSIONED AT LEAST ';,I5)
375 FORMAT (10X,';CCKWRK MUST BE DIMENSIONED AT LEAST ';,I5)
475 FORMAT (10X,';CHARACTER LENGTH OF CCKWRK MUST BE AT LEAST 16 ';)
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKINDX (ICKWRK, RCKWRK, MM, KK, II, NFIT)
C
C START PROLOGUE
C
C SUBROUTINE CKINDX (ICKWRK, RCKWRK, MM, KK, II, NFIT)*
C Returns a group of indices defining the size of the particular
C reaction mechanism
C
C INPUT
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C MM - Total number of elements in mechanism.
C Data type - integer scalar
C KK - Total number of species in mechanism.
C Data type - integer scalar
C II - Total number of reactions in mechanism.
C Data type - integer scalar
C NFIT - number of coefficients in fits to thermodynamic data
C for one temperature range; NFIT = number of
C coefficients in polynomial fits to CP/R + 2.
C Data type - integer scalar
C
C END PROLOGUE
C
IMPLICIT NONE
C
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,II,KK,MM,NFIT
DOUBLE PRECISION
1 RCKWRK
DIMENSION ICKWRK(*), RCKWRK(*)
C
MM = NMM
KK = NKK
II = NII
NFIT = NCP2
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKRHOY (P, T, Y, ICKWRK, RCKWRK, RHO)
C
C START PROLOGUE
C
C SUBROUTINE CKRHOY (P, T, Y, ICKWRK, RCKWRK, RHO)
C Returns the mass density of the gas mixture given the pressure,
C temperature and mass fractions; see Eq. (2).
C
C INPUT
C P - Pressure.
C cgs units - dynes/cm**2
C Data type - real scalar
C T - Temperature.
C cgs units - K
C Data type - real scalar
C Y - Mass fractions of the species.
C cgs units - none
C Data type - real array
C Dimension Y(*) at least KK, the total number of
C species.
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C RHO - Mass density.
C cgs units - gm/cm**3
C Data type - real scalar
C
C END PROLOGUE
C
IMPLICIT NONE
C
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K
DOUBLE PRECISION
1 P,RCKWRK,RHO,SUMYOW,T,Y
DIMENSION ICKWRK(*), RCKWRK(*), Y(*)
C
SUMYOW = 0.D00
DO 150 K = 1, NKK
SUMYOW = SUMYOW + Y(K)/RCKWRK(NcWT + K - 1)
150 CONTINUE
RHO = P / (SUMYOW*T*RCKWRK(NcRU))
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKYTX (Y, ICKWRK, RCKWRK, X)
C
C START PROLOGUE
C
C SUBROUTINE CKYTX (Y, ICKWRK, RCKWRK, X)
C Returns the mole fractions given the mass fractions; see Eq. (6).
C
C INPUT
C Y - Mass fractions of the species.
C cgs units - none
C Data type - real array
C Dimension Y(*) at least KK, the total number of
C species.
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C X - Mole fractions of the species.
C cgs units - none
C Data type - real array
C Dimension X(*) at least KK, the total number of
C species.
C
C END PROLOGUE
C
IMPLICIT NONE
C
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K
DOUBLE PRECISION
1 RCKWRK,SUMYOW,X,Y
DIMENSION ICKWRK(*), RCKWRK(*), Y(*), X(*)
C
SUMYOW = 0.D00
DO 150 K = 1, NKK
SUMYOW = SUMYOW + Y(K)/RCKWRK(NcWT + K - 1)
150 CONTINUE
DO 200 K = 1, NKK
X(K) = Y(K)/(SUMYOW*RCKWRK(NcWT + K - 1))
200 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKWYP (P, T, Y, ICKWRK, RCKWRK, WDOT)
C
C START PROLOGUE
C
C SUBROUTINE CKWYP (P, T, Y, ICKWRK, RCKWRK, WDOT)
C Returns the molar production rates of the species given the
C pressure, temperature and mass fractions; see Eq. (49).
C
C INPUT
C P - Pressure.
C cgs units - dynes/cm**2
C Data type - real scalar
C T - Temperature.
C cgs units - K
C Data type - real scalar
C Y - Mass fractions of the species.
C cgs units - none
C Data type - real array
C Dimension Y(*) at least KK, the total number of
C species.
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C WDOT - Chemical molar production rates of the species.
C cgs units - moles/(cm**3*sec)
C Data type - real array
C Dimension WDOT(*) at least KK, the total number of
C species.
C
C END PROLOGUE
C
IMPLICIT NONE
C
INCLUDE ';ckstrt.h';
C
INTEGER I,ICKWRK,K,L,N
DOUBLE PRECISION
1 P,RCKWRK,ROP,T,WDOT,Y
DIMENSION ICKWRK(*), RCKWRK(*), Y(*), WDOT(*)
C
CALL CKRATT (RCKWRK, ICKWRK, NII, MXSP, RCKWRK(NcRU),
1 RCKWRK(NcPA), T, ICKWRK(IcNS), ICKWRK(IcNU),
2 ICKWRK(IcNK), NPAR+1, RCKWRK(NcCO), NREV,
3 ICKWRK(IcRV), RCKWRK(NcRV), NLAN, NLAR, ICKWRK(IcLT),
4 RCKWRK(NcLT), NRLT, ICKWRK(IcRL), RCKWRK(NcRL),
5 RCKWRK(NcK1), NRNU, ICKWRK(IcRNU), RCKWRK(NcRNU),
6 NEIM, ICKWRK(IcEI), ICKWRK(IcTD), NJAN, NJAR,
7 ICKWRK(IcJN), RCKWRK(NcJN), NFT1, NF1R,
8 ICKWRK(IcF1), RCKWRK(NcF1),
9 RCKWRK(NcKF), RCKWRK(NcKR), RCKWRK(NcI1))
C
CALL CKYTCP (P, T, Y, ICKWRK, RCKWRK, RCKWRK(NcK1))
C
CALL CKRATX (NII, NKK, MXSP, MXTB, T, RCKWRK(NcK1),
1 ICKWRK(IcNU), ICKWRK(IcNK), NPAR+1, RCKWRK(NcCO),
2 NFAL, ICKWRK(IcFL), ICKWRK(IcFO), ICKWRK(IcKF), NFAR,
3 RCKWRK(NcFL), NTHB, ICKWRK(IcTB), ICKWRK(IcKN),
4 RCKWRK(NcKT), ICKWRK(IcKT), RCKWRK(NcKF),
5 RCKWRK(NcKR), RCKWRK(NcI1), RCKWRK(NcI2),
6 RCKWRK(NcI3), NRNU, ICKWRK(IcRNU), RCKWRK(NcRNU),
7 NORD, ICKWRK(IcORD), MXORD, ICKWRK(IcKOR),
8 RCKWRK(NcKOR))
C
DO 50 K = 1, NKK
WDOT(K) = 0.D00
50 CONTINUE
DO 100 N = 1, MXSP
DO 100 I = 1, NII
K = ICKWRK(IcNK + (I-1)*MXSP + N - 1)
IF (K .NE. 0) THEN
ROP = RCKWRK(NcI1+I-1) - RCKWRK(NcI2+I-1)
WDOT(K) = WDOT(K) + ROP *
1 DBLE(ICKWRK(IcNU + (I-1)*MXSP + N - 1))
ENDIF
100 CONTINUE
C
IF (NRNU .LE. 0) RETURN
DO 200 L = 1, NRNU
I = ICKWRK(IcRNU + L - 1)
ROP = RCKWRK(NcI1+I-1) - RCKWRK(NcI2+I-1)
DO 200 N = 1, MXSP
K = ICKWRK(IcNK + (I-1)*MXSP + N - 1)
IF (K .NE. 0) WDOT(K) = WDOT(K) + ROP *
1 RCKWRK(NcRNU + (L-1)*MXSP + N - 1)
200 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKCPBS (T, Y, ICKWRK, RCKWRK, CPBMS)
C
C START PROLOGUE
C
C SUBROUTINE CKCPBS (T, Y, ICKWRK, RCKWRK, CPBMS)
C Returns the mean specific heat at constant pressure;
C see Eq. (34).
C
C INPUT
C T - Temperature.
C cgs units - K
C Data type - real scalar
C Y - Mass fractions of the species.
C cgs units - none
C Data type - real array
C Dimension Y(*) at least KK, the total number of
C species.
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C CPBMS - Mean specific heat at constant pressure in mass units.
C cgs units - ergs/(gm*K)
C Data type - real scalar
C
C END PROLOGUE
C
IMPLICIT NONE
C
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K
DOUBLE PRECISION
1 CPBMS,RCKWRK,T,Y
DIMENSION ICKWRK(*), RCKWRK(*), Y(*)
C
CALL CKCPMS (T, ICKWRK, RCKWRK, RCKWRK(NcK1))
C
CPBMS = 0.D00
DO 100 K = 1, NKK
CPBMS = CPBMS + Y(K)*RCKWRK(NcK1 + K - 1)
100 CONTINUE
C
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKHBMS (T, Y, ICKWRK, RCKWRK, HBMS)
C
C START PROLOGUE
C
C SUBROUTINE CKHBMS (T, Y, ICKWRK, RCKWRK, HBMS)
C Returns the mean enthalpy of the mixture in mass units;
C see Eq. (38).
C
C INPUT
C T - Temperature.
C cgs units - K
C Data type - real scalar
C Y - Mass fractions of the species.
C cgs units - none
C Data type - real array
C Dimension Y(*) at least KK, the total number of
C species.
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C HBMS - Mean enthalpy in mass units.
C cgs units - ergs/gm
C Data type - real scalar
C
C END PROLOGUE
C
IMPLICIT NONE
C
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K
DOUBLE PRECISION
1 HBMS,RCKWRK,T,Y
DIMENSION ICKWRK(*), RCKWRK(*), Y(*)
C
CALL CKHMS (T, ICKWRK, RCKWRK, RCKWRK(NcK1))
C
HBMS = 0.D00
DO 100 K = 1, NKK
HBMS = HBMS + Y(K)*RCKWRK(NcK1 + K - 1)
100 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKMMWY (Y, ICKWRK, RCKWRK, WTM)
C
C START PROLOGUE
C
C SUBROUTINE CKMMWY (Y, ICKWRK, RCKWRK, WTM)
C Returns the mean molecular weight of the gas mixture given the
C mass fractions; see Eq. (3).
C
C INPUT
C Y - Mass fractions of the species.
C cgs units - none
C Data type - real array
C Dimension Y(*) at least KK, the total number of
C species.
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C WTM - Mean molecular weight of the species mixture.
C cgs units - gm/mole
C Data type - real scalar
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K
DOUBLE PRECISION
1 RCKWRK,SUMYOW,WTM,Y
DIMENSION Y(*), ICKWRK(*), RCKWRK(*)
C
SUMYOW=0.D00
DO 150 K = 1, NKK
SUMYOW = SUMYOW + Y(K)/RCKWRK(NcWT + K - 1)
150 CONTINUE
WTM = 1.D00/SUMYOW
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKRP (ICKWRK, RCKWRK, RU, RUC, PA)
C
C START PROLOGUE
C
C SUBROUTINE CKRP (ICKWRK, RCKWRK, RU, RUC, PA)
C Returns universal gas constants and the pressure of one standard
C atmosphere
C
C INPUT
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C RU - Universal gas constant.
C cgs units - 8.314510E7 ergs/(mole*K)
C Data type - real scalar
C RUC - Universal gas constant used only in conjuction with
C activation energy.
C preferred units - RU / 4.184 cal/(mole*K)
C Data type - real scalar
C PA - Pressure of one standard atmosphere.
C cgs units - 1.01325E6 dynes/cm**2
C Data type - real scalar
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK
DOUBLE PRECISION
1 PA,RCKWRK,RU,RUC
DIMENSION ICKWRK(*), RCKWRK(*)
C
RU = RCKWRK(NcRU)
RUC = RCKWRK(NcRC)
PA = RCKWRK(NcPA)
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKSYME (CCKWRK, LOUT, ENAME, KERR)
C
C START PROLOGUE
C
C SUBROUTINE CKSYME (CCKWRK, LOUT, ENAME, KERR)*
C Returns the character strings of element names.
C
C INPUT
C CCKWRK - Array of character work space.
C Data type - CHARACTER*16 array
C Dimension CCKWRK(*) at least LENCWK.
C LOUT - Output unit for printed diagnostics.
C Data type - integer scalar
C
C OUTPUT
C ENAME - Element names.
C Data type - CHARACTER*(*) array
C Dimension ENAME at least MM, the total number of
C elements in the problem.
C KERR - Error flag; character length error will result in
C KERR = .TRUE.
C Data type - logical
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ILASCH,ILEN,LOUT,LT,M
CHARACTER*(*) CCKWRK(*), ENAME(*)
LOGICAL KERR
C
KERR = .FALSE.
ILEN = LEN(ENAME(1))
DO 150 M = 1, NMM
LT = ILASCH(CCKWRK(IcMM+M-1))
ENAME(M) = '; ';
IF (LT .LE. ILEN) THEN
ENAME(M) = CCKWRK(IcMM+M-1)
ELSE
WRITE (LOUT,';(A)';)
1 '; Error in CKSYME...character string length too small ';
KERR = .TRUE.
ENDIF
150 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKRATT (RCKWRK, ICKWRK, II, MAXSP, RU, PATM, T, NSPEC,
1 NU, NUNK, NPAR, PAR, NREV, IREV, RPAR, NLAN,
2 NLAR, ILAN, PLT, NRLT, IRLT, RPLT, SMH, NRNU,
3 IRNU, RNU, NEIM, IEIM, ITDEP, NJAN, NJAR, IJAN,
4 PJAN, NFT1, NF1R, IFT1, PF1, RKFT, RKRT, EQK)
C
C START PROLOGUE
C
C SUBROUTINE CKRATT (RCKWRK, ICKWRK, II, MAXSP, RU, PATM, T, NSPEC,
C 1 NU, NUNK, NPAR, PAR, NREV, IREV, RPAR, NLAN,
C 2 NLAR, ILAN, PLT, NRLT, IRLT, RPLT, SMH, NRNU,
C 3 IRNU, RNU, NEIM, IEIM, ITDEP, NJAN, NJAR, IJAN,
C 4 PJAN, NFT1, NF1R, IFT1, PF1, RKFT, RKRT, EQK)
C
C END PROLOGUE
C
IMPLICIT NONE
C
INTEGER I,I2,ICKWRK,IEIM,IFT1,II,IJAN,ILAN,IREV,IRLT,IRNU,ITDEP,
1 J,L,MAXSP,N,N2,NEIM,NF1R,NFT1,NJAN,NJAR,NLAN,NLAR,NPAR,
2 NREV,NRLT,NRNU,NSPEC,NU,NUNK,NUSUMK
DOUBLE PRECISION
1 ALOGT,BIG,EQK,EXPARG,PAR,PATM,PF1,PFAC,PFAC2,PFR,PJAN,
2 PLT,POWR1,POWR2,RCKWRK,RKFT,RKRT,RNU,RNUSUM,ROOTJ,
3 RPAR,RPLT,RU,SMALL,SMH,SUMJ,SUMSMH,T,TEMP,TEV,TFAC
DIMENSION RCKWRK(*), ICKWRK(*), NSPEC(*), NU(MAXSP,*),
1 NUNK(MAXSP,*), PAR(NPAR,*), IREV(*), RPAR(NPAR,*),
2 ILAN(*), IRLT(*), PLT(NLAR,*), RPLT(NLAR,*), SMH(*),
3 RKFT(*), RKRT(*), EQK(*), IRNU(*), RNU(MAXSP,*),
4 T(*), IEIM(*), ITDEP(*), IJAN(*), PJAN(NJAR,*),
5 IFT1(*), PF1(NF1R,*)
C
COMMON /MACH/ SMALL,BIG,EXPARG
C
ALOGT = DLOG(T(1))
C
DO 20 I = 1, II
RKFT(I) = PAR(1,I) *DEXP(PAR(2,I)*ALOGT - PAR(3,I)/T(1))
20 CONTINUE
C
C Landau-Teller reactions
C
POWR1 = 1.D00/3.D00
POWR2 = 2.D00/3.D00
DO 25 N = 1, NLAN
I = ILAN(N)
TFAC = PLT(1,N)/T(1)**POWR1 + PLT(2,N)/T(1)**POWR2
RKFT(I) = RKFT(I) * DEXP(TFAC)
25 CONTINUE
C
CALL CKSMH (T(1), ICKWRK, RCKWRK, SMH)
DO 50 I = 1, II
SUMSMH = 0.D00
DO 40 N = 1, MAXSP
IF (NUNK(N,I).NE.0)
1 SUMSMH = SUMSMH + DBLE(NU(N,I))*SMH(NUNK(N,I))
40 CONTINUE
IF (SUMSMH .NE. 0.0) EQK(I) = DEXP(DMIN1(SUMSMH,EXPARG))
50 CONTINUE
C
DO 55 N = 1, NRNU
SUMSMH = 0.D00
I = IRNU(N)
DO 45 L = 1, MAXSP
IF (NUNK(L,I).NE.0) SUMSMH=SUMSMH+RNU(L,N)*SMH(NUNK(L,I))
45 CONTINUE
IF (SUMSMH .NE. 0.0) EQK(I) = DEXP(DMIN1(SUMSMH,EXPARG))
55 CONTINUE
C
PFAC = PATM / (RU*T(1))
DO 60 I = 1, II
NUSUMK = NU(1,I)+NU(2,I)+NU(3,I)+NU(4,I)+NU(5,I)+NU(6,I)
EQK(I) = EQK(I) * PFAC**NUSUMK
60 CONTINUE
DO 65 N = 1, NRNU
RNUSUM = RNU(1,N)+RNU(2,N)+RNU(3,N)+RNU(4,N)+RNU(5,N)+RNU(6,N)
I = IRNU(N)
PFR = PFAC ** RNUSUM
EQK(I) = EQK(I) * PFR
65 CONTINUE
C
DO 68 I = 1, II
C
C RKR=0.0 for irreversible reactions, else RKR=RKF/MAX(EQK,SMALL)
C
RKRT(I) = 0.D00
IF (NSPEC(I).GT.0) RKRT(I) = RKFT(I) / DMAX1(EQK(I),SMALL)
68 CONTINUE
C
C if reverse parameters have been given:
C
DO 70 N = 1, NREV
I = IREV(N)
RKRT(I) = RPAR(1,N) * DEXP(RPAR(2,N)*ALOGT - RPAR(3,N)/T(1))
EQK(I) = RKFT(I)/RKRT(I)
70 CONTINUE
C
C if reverse Landau-Teller parameters have been given:
C
POWR1 = 1.D00/3.D00
POWR2 = 2.D00/3.D00
DO 75 N = 1, NRLT
I = IRLT(N)
TFAC = RPLT(1,N)/T(1)**POWR1 + RPLT(2,N)/T(1)**POWR2
RKRT(I) = RKRT(I) * DEXP(TFAC)
EQK(I) = RKFT(I)/RKRT(I)
75 CONTINUE
C
C electron-impact reactions
C
DO 85 N = 1, NEIM
I = IEIM(N)
RKRT(I) = 0.D00
TEMP = T(ITDEP(N))
RKFT(I) = PAR(1,I) * DEXP(PAR(2,I) * DLOG(TEMP)
1 - PAR(3,I)/TEMP)
PFAC2 = PATM/(RU*TEMP)
NUSUMK = NU(1,I)+NU(2,I)+NU(3,I)+NU(4,I)+NU(5,I)+NU(6,I)
EQK(I) = EQK(I) * (PFAC2/PFAC)**NUSUMK
C
DO 80 L = 1, NRNU
IF (IRNU(L) .EQ. I) THEN
RNUSUM=RNU(1,L)+RNU(2,L)+RNU(3,L)+RNU(4,L)+RNU(5,L)+
1 RNU(6,L)
PFR = (PFAC2/PFAC) ** RNUSUM
EQK(I) = EQK(I) * PFR
ENDIF
80 CONTINUE
C
IF (NSPEC(I) .GT. 0) RKRT(I) = RKFT(I)/DMAX1(EQK(I),SMALL)
DO 83 N2 = 1, NREV
I2 = IREV(N2)
IF (I2.EQ.I) THEN
RKRT(I) = RPAR(1,N2) * DEXP(RPAR(2,N2) *DLOG(TEMP)
1 - RPAR(3,N2)/TEMP)
EQK(I) = RKFT(I)/DMAX1(RKRT(I),SMALL)
ENDIF
83 CONTINUE
85 CONTINUE
C
C jannev, langer, evans & post - type reactions
C
DO 95 N = 1, NJAN
I = IJAN(N)
RKRT(I) = 0.D00
C
C CONVERT E- TEMPERATURE TO eV';s
C
TEMP = T(2)
TEV = TEMP/ 11600.0D00
RKFT(I) = PAR(1,I) * DEXP(PAR(2,I) * DLOG(TEMP)
1 - PAR(3,I)/TEMP)
SUMJ = 0.D00
DO 90 J = 1, NJAR
SUMJ = SUMJ + PJAN(J,N) * (DLOG(TEV))**(J-1)
90 CONTINUE
RKFT(I) = RKFT(I) * DEXP(SUMJ)
95 CONTINUE
C
C reactions using fit#1: k = A * T^B * exp(v1/T+v2/T^2+v3/T^3...)
C
DO 105 N = 1, NFT1
I = IFT1(N)
RKRT(I) = 0.D00
C
C CHECK IF REACTION IS ALSO AN ELECTRON-IMPACT REAX
C
TEMP = T(1)
DO 100 N2 = 1, NEIM
IF (IEIM(N2) .EQ. I) TEMP = T(ITDEP(N2))
100 CONTINUE
C
RKFT(I) = PAR(1,I) * TEMP**PAR(2,I)
SUMJ = 0.D00
DO 103 J = 1, NF1R
ROOTJ = 1.0D00/DFLOAT(J)
IF (TEMP.GE.BIG**ROOTJ .OR. SUMJ .GE. DLOG(BIG)) THEN
SUMJ = DLOG(BIG)
ELSE
SUMJ = SUMJ + PF1(J,N)/TEMP**DFLOAT(J)
ENDIF
103 CONTINUE
IF (SUMJ .GE. DLOG(BIG)) SUMJ = DLOG(BIG)
RKFT(I) = DMIN1(BIG, RKFT(I) * DEXP(SUMJ))
105 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKYTCP (P, T, Y, ICKWRK, RCKWRK, C)
C
C START PROLOGUE
C
C SUBROUTINE CKYTCP (P, T, Y, ICKWRK, RCKWRK, C)
C Returns the molar concentrations given the pressure,
C temperature and mass fractions; see Eq. (7).
C
C INPUT
C P - Pressure.
C cgs units - dynes/cm**2
C Data type - real scalar
C T - Temperature.
C cgs units - K
C Data type - real scalar
C Y - Mass fractions of the species.
C cgs units - none
C Data type - real array
C Dimension Y(*) at least KK, the total number of
C species.
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C C - Molar concentrations of the species.
C cgs units - mole/cm**3
C Data type - real array
C Dimension C(*) at least KK, the total number of
C species.
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K
DOUBLE PRECISION
1 C,P,RCKWRK,SUMYOW,T,Y
DIMENSION ICKWRK(*), RCKWRK(*), Y(*), C(*)
C
SUMYOW = 0.D00
DO 150 K = 1, NKK
SUMYOW = SUMYOW + Y(K) / RCKWRK(NcWT + K - 1)
150 CONTINUE
SUMYOW = SUMYOW * T * RCKWRK(NcRU)
DO 200 K = 1, NKK
C(K) = P * Y(K) / (SUMYOW * RCKWRK(NcWT + K - 1))
200 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKRATX (II, KK, MAXSP, MAXTB, T, C, NU, NUNK, NPAR,
1 PAR, NFAL, IFAL, IFOP, KFAL, NFAR, FPAR, NTHB,
2 ITHB, NTBS, AIK, NKTB, RKFT, RKRT, RKF, RKR,
3 CTB, NRNU, IRNU, RNU, NORD, IORD, MXORD, KORD,
4 RORD)
C
C START PROLOGUE
C
C SUBROUTINE CKRATX (II, KK, MAXSP, MAXTB, T, C, NU, NUNK, NPAR,
C 1 PAR, NFAL, IFAL, IFOP, KFAL, NFAR, FPAR, NTHB,
C 2 ITHB, NTBS, AIK, NKTB, RKFT, RKRT, RKF, RKR,
C 3 CTB, NRNU, IRNU, RNU, NORD, IORD, MXORD, KORD,
C 4 RORD)
C
C END PROLOGUE
C
IMPLICIT NONE
C
INTEGER I,IFAL,IFOP,II,IORD,IRNU,ITHB,K,KFAL,KK,KORD,L,MAXSP,
1 MAXTB,MXORD,N,NFAL,NFAR,NK,NKTB,NORD,NPAR,NRNU,NTBS,
2 NTHB,NU,NUNK
DOUBLE PRECISION
1 AIK,ALOGT,BIG,C,C1,C2,C3,C4,C5,C6,CNK,CPRLOG,CTB,CTOT,
2 EXPARG,FC,FCENT,FCLOG,FLOG,FPAR,PAR,PCOR,PR,PRLOG,RKF,
3 RKFT,RKLOW,RKR,RKRT,RNU,RORD,SMALL,T,XN,XP
DIMENSION C(*), NU(MAXSP,*), NUNK(MAXSP,*), PAR(NPAR,*),
1 IFAL(*), IFOP(*), KFAL(*), FPAR(NFAR,*), ITHB(*),
2 NTBS(*), AIK(MAXTB,*), NKTB(MAXTB,*), RKFT(*),
3 RKRT(*), RKF(*), RKR(*), CTB(*), IRNU(*), RNU(MAXSP,*),
4 IORD(*), KORD(MXORD,*), RORD(MXORD,*)
C
COMMON /MACH/ SMALL,BIG,EXPARG
C
DO 20 I = 1, II
CTB(I) = 1.D00
RKF(I) = 0.D00
RKR(I) = 0.D00
20 CONTINUE
C
C third-body reactions
C
IF (NTHB .GT. 0) THEN
CTOT = 0.D00
DO 10 K = 1, KK
CTOT = CTOT + C(K)
10 CONTINUE
DO 80 N = 1, NTHB
CTB(ITHB(N)) = CTOT
DO 80 L = 1, NTBS(N)
CTB(ITHB(N)) = CTB(ITHB(N))+(AIK(L,N)-1.D00)*C(NKTB(L,N))
80 CONTINUE
ENDIF
C
C If fall-off (pressure correction):
C
IF (NFAL .GT. 0) THEN
ALOGT = DLOG(T)
C
DO 90 N = 1, NFAL
C
RKLOW = FPAR(1,N) * DEXP(FPAR(2,N)*ALOGT - FPAR(3,N)/T)
C
C CONCENTRATION OF THIRD BODY
C
IF (KFAL(N) .EQ. 0) THEN
PR = RKLOW * CTB(IFAL(N)) / RKFT(IFAL(N))
CTB(IFAL(N)) = 1.D00
ELSE
PR = RKLOW * C(KFAL(N)) / RKFT(IFAL(N))
ENDIF
C
PCOR = PR / (1.D00 + PR)
C
IF (IFOP(N) .GT. 1) THEN
PRLOG = DLOG10(DMAX1(PR,SMALL))
C
IF (IFOP(N) .EQ. 2) THEN
C
C 8-PARAMETER SRI FORM
C
XP = 1.D00/(1.D00 + PRLOG**2)
FC = ((FPAR(4,N)*DEXP(-FPAR(5,N)/T)
1 + DEXP(-T/FPAR(6,N))) **XP)
2 * FPAR(7,N) * T**FPAR(8,N)
C
ELSE
C
C 6-PARAMETER TROE FORM
C
FCENT = (1.D00-FPAR(4,N)) * DEXP(-T/FPAR(5,N))
1 + FPAR(4,N) * DEXP(-T/FPAR(6,N))
C
C 7-PARAMETER TROE FORM
C
IF (IFOP(N) .EQ. 4) FCENT = FCENT + DEXP(-FPAR(7,N)/T)
C
FCLOG = DLOG10(DMAX1(FCENT,SMALL))
XN = 0.75D00 - 1.27D00*FCLOG
CPRLOG= PRLOG - (0.4D00 + 0.67D00*FCLOG)
FLOG = FCLOG/(1.D00 + (CPRLOG/(XN-0.14D00*CPRLOG))**2)
FC = 10.D00**FLOG
ENDIF
PCOR = FC * PCOR
ENDIF
C
RKFT(IFAL(N)) = RKFT(IFAL(N)) * PCOR
RKRT(IFAL(N)) = RKRT(IFAL(N)) * PCOR
90 CONTINUE
ENDIF
C
C Multiply by the product of reactants and product of products
C PAR(4,I) is a perturbation factor
C
DO 150 I = 1, II
RKFT(I) = RKFT(I) * CTB(I) * PAR(4,I)
RKRT(I) = RKRT(I) * CTB(I) * PAR(4,I)
C
IF (NU(1,I) .NE. 0) THEN
RKF(I) = RKFT(I)*C(NUNK(1,I))**IABS(NU(1,I))
RKR(I) = RKRT(I)*C(NUNK(4,I))**NU(4,I)
IF (NUNK(2,I) .NE. 0) THEN
RKF(I)= RKF(I) * C(NUNK(2,I))**IABS(NU(2,I))
IF (NUNK(3,I) .NE. 0)
1 RKF(I) = RKF(I) * C(NUNK(3,I))**IABS(NU(3,I))
ENDIF
IF (NUNK(5,I) .NE. 0) THEN
RKR(I) = RKR(I) * C(NUNK(5,I))**NU(5,I)
IF (NUNK(6,I) .NE. 0) RKR(I)=RKR(I)*C(NUNK(6,I))**NU(6,I)
ENDIF
ENDIF
150 CONTINUE
C
DO 160 N = 1, NRNU
I = IRNU(N)
C1 = C(NUNK(1,I)) ** DABS(RNU(1,N))
C4 = C(NUNK(4,I)) ** RNU(4,N)
RKF(I) = RKFT(I) * C1
RKR(I) = RKRT(I) * C4
IF (NUNK(2,I) .NE. 0) THEN
C2 = C(NUNK(2,I)) ** DABS(RNU(2,N))
RKF(I) = RKF(I) * C2
IF (NUNK(3,I) .NE. 0) THEN
C3 = C(NUNK(3,I)) ** DABS(RNU(3,N))
RKF(I) = RKF(I) * C3
ENDIF
ENDIF
IF (NUNK(5,I) .NE. 0) THEN
C5 = C(NUNK(5,I)) ** RNU(5,N)
RKR(I) = RKR(I) * C5
IF (NUNK(6,I) .NE. 0) THEN
C6 = C(NUNK(6,I)) ** RNU(6,N)
RKR(I) = RKR(I) * C6
ENDIF
ENDIF
160 CONTINUE
C
DO 200 N = 1, NORD
I = IORD(N)
RKF(I) = RKFT(I)
RKR(I) = RKRT(I)
C
DO 190 L = 1, MXORD
NK = KORD(L,N)
IF (NK .LT. 0) THEN
NK = IABS(NK)
CNK = C(NK) ** RORD(L,N)
RKF(I) = RKF(I) * CNK
ELSEIF (NK .GT. 0) THEN
CNK = C(NK) ** RORD(L,N)
RKR(I) = RKR(I) * CNK
ENDIF
190 CONTINUE
200 CONTINUE
C
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKCPMS (T, ICKWRK, RCKWRK, CPMS)
C
C START PROLOGUE
C
C SUBROUTINE CKCPMS (T, ICKWRK, RCKWRK, CPMS)
C Returns the specific heats at constant pressure in mass units;
C see Eq. (26).
C
C INPUT
C T - Temperature.
C cgs units - K
C Data type - real scalar
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C CPMS - Specific heats at constant pressure in mass units
C for the species.
C cgs units - ergs/(gm*K)
C Data type - real array
C Dimension CPMS(*) at least KK, the total number of
C species.
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K,L,N,NA1
DOUBLE PRECISION
1 CPMS,RCKWRK,SUM,T,TEMP,TN
DIMENSION ICKWRK(*), RCKWRK(*), CPMS(*), TN(10)
C
TN(1) = 1.D00
DO 150 N = 2, NCP
TN(N) = T**(N-1)
150 CONTINUE
C
DO 250 K = 1, NKK
L = 1
DO 220 N = 2, ICKWRK(IcNT + K - 1)-1
TEMP = RCKWRK(NcTT + (K-1)*MXTP + N - 1)
IF (T .GT. TEMP) L = L+1
220 CONTINUE
C
NA1 = NcAA + (L-1)*NCP2 + (K-1)*NCP2T
SUM = 0.D00
DO 240 N = 1, NCP
SUM = SUM + TN(N)*RCKWRK(NA1 + N - 1)
240 CONTINUE
CPMS(K) = RCKWRK(NcRU) * SUM / RCKWRK(NcWT + K - 1)
C
250 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKHMS (T, ICKWRK, RCKWRK, HMS)
C
C START PROLOGUE
C
C SUBROUTINE CKHMS (T, ICKWRK, RCKWRK, HMS)
C Returns the enthalpies in mass units; see Eq. (27).
C
C INPUT
C T - Temperature.
C cgs units - K
C Data type - real scalar
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C OUTPUT
C HMS - Enthalpies in mass units for the species.
C cgs units - ergs/gm
C Data type - real array
C Dimension HMS(*) at least KK, the total number of
C species.
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K,L,N,NA1
DOUBLE PRECISION
1 HMS,RCKWRK,RUT,SUM,T,TEMP,TN
DIMENSION ICKWRK(*), RCKWRK(*), HMS(*), TN(10)
C
RUT = T*RCKWRK(NcRU)
TN(1)=1.D00
DO 150 N = 2, NCP
TN(N) = T**(N-1)/N
150 CONTINUE
C
DO 250 K = 1, NKK
L = 1
DO 220 N = 2, ICKWRK(IcNT + K - 1)-1
TEMP = RCKWRK(NcTT + (K-1)*MXTP + N - 1)
IF (T .GT. TEMP) L = L+1
220 CONTINUE
C
NA1 = NcAA + (L-1)*NCP2 + (K-1)*NCP2T
SUM = 0.D00
DO 225 N = 1, NCP
SUM = SUM + TN(N)*RCKWRK(NA1 + N - 1)
225 CONTINUE
HMS(K) = RUT * (SUM + RCKWRK(NA1 + NCP1 - 1)/T)
1 / RCKWRK(NcWT + K - 1)
250 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKSMH (T, ICKWRK, RCKWRK, SMH)
C
C START PROLOGUE
C
C SUBROUTINE CKSMH (T, ICKWRK, RCKWRK, SMH)*
C Returns the array of entropies minus enthalpies for the species.
C It is normally not called directly by the user.
C
C INPUT
C T - Temperature.
C cgs units - K
C Data type - real scalar
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C SMH - Entropy minus enthalpy for the species,
C SMH(K) = S(K)/R - H(K)/RT.
C cgs units - none
C Data type - real array
C Dimension SMH(*) at least KK, the total number of
C species.
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K,L,N,NA1
DOUBLE PRECISION
1 RCKWRK,SMH,SUM,T,TEMP,TN
DIMENSION ICKWRK(*), RCKWRK(*), SMH(*), TN(10)
C
TN(1) = DLOG(T) - 1.D00
DO 150 N = 2, NCP
TN(N) = T**(N-1)/((N-1)*N)
150 CONTINUE
C
DO 250 K = 1, NKK
L = 1
DO 220 N = 2, ICKWRK(IcNT + K - 1)-1
TEMP = RCKWRK(NcTT + (K-1)*MXTP + N - 1)
IF (T .GT. TEMP) L = L+1
220 CONTINUE
C
NA1 = NcAA + (L-1)*NCP2 + (K-1)*NCP2T
SUM = 0.D00
DO 225 N = 1, NCP
SUM = SUM + TN(N)*RCKWRK(NA1 + N - 1)
225 CONTINUE
SMH(K) = SUM + RCKWRK(NA1 + NCP2 - 1)
1 - RCKWRK(NA1 + NCP1 - 1)/T
C
250 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
FUNCTION IPPLEN (LINE)
C
C BEGIN PROLOGUE
C
C FUNCTION IPPLEN (LINE)
C Returns the effective length of a character string, i.e.,
C the index of the last character before an exclamation mark (!)
C indicating a comment.
C
C INPUT
C LINE - A character string.
C Data type - CHARACTER*(*)
C
C OUTPUT
C IPPLEN - The effective length of the character string.
C Data type - integer scalar
C
C END PROLOGUE
C
IMPLICIT NONE
C
INTEGER IFIRCH,ILASCH,IN,IPPLEN
CHARACTER LINE*(*)
C
IN = IFIRCH(LINE)
IF (IN.EQ.0 .OR. LINE(IN:IN) .EQ. ';!';) THEN
IPPLEN = 0
ELSE
IN = INDEX(LINE,';!';)
IF (IN .EQ. 0) THEN
IPPLEN = ILASCH(LINE)
ELSE
IPPLEN = ILASCH(LINE(:IN-1))
ENDIF
ENDIF
RETURN
END
C
C----------------------------------------------------------------------C
C C
FUNCTION IFIRCH (STRING)
C BEGIN PROLOGUE IFIRCH
C DATE WRITTEN 850626
C REVISION DATE 850626
C CATEGORY NO. M4.
C KEYWORDS CHARACTER STRINGS,SIGNIFICANT CHARACTERS
C AUTHOR CLARK,G.L.,GROUP C-3 LOS ALAMOS NAT';L LAB
C PURPOSE Determines first significant (non-blank) character
C in character variable
C DESCRIPTION
C
C-----------------------------------------------------------------------
C IFIRCH locates the first non-blank character in a string of
C arbitrary length. If no characters are found, IFIRCH is set = 0.
C When used with the companion routine ILASCH, the length of a string
C can be determined, and/or a concatenated substring containing the
C significant characters produced.
C-----------------------------------------------------------------------
C
C REFERENCES (NONE)
C ROUTINES CALLED (NONE)
C END PROLOGUE IFIRCH
IMPLICIT NONE
C
INTEGER I,IFIRCH,NLOOP
CHARACTER* (*)STRING
C
C FIRST EXECUTABLE STATEMENT IFIRCH
NLOOP = LEN(STRING)
C
IF (NLOOP.EQ.0 .OR. STRING.EQ.'; ';) THEN
IFIRCH = 0
RETURN
ENDIF
C
DO 100 I = 1, NLOOP
IF (STRING(I:I) .NE. '; ';) GO TO 120
100 CONTINUE
C
IFIRCH = 0
RETURN
120 CONTINUE
IFIRCH = I
END
C
FUNCTION ILASCH (STRING)
C BEGIN PROLOGUE ILASCH
C DATE WRITTEN 850626
C REVISION DATE 850626
C CATEGORY NO. M4.
C KEYWORDS CHARACTER STRINGS,SIGNIFICANT CHARACTERS
C AUTHOR CLARK,G.L.,GROUP C-3 LOS ALAMOS NAT';L LAB
C PURPOSE Determines last significant (non-blank) character
C in character variable
C DESCRIPTION
C
C-----------------------------------------------------------------------
C IFIRCH locates the last non-blank character in a string of
C arbitrary length. If no characters are found, ILASCH is set = 0.
C When used with the companion routine IFIRCH, the length of a string
C can be determined, and/or a concatenated substring containing the
C significant characters produced.
C Note that the FORTRAN intrinsic function LEN returns the length
C of a character string as declared, rather than as filled. The
C declared length includes leading and trailing blanks, and thus is
C not useful in generating ';significant'; substrings.
C-----------------------------------------------------------------------
C
C REFERENCES (NONE)
C ROUTINES CALLED (NONE)
C END PROLOGUE IFIRCH
IMPLICIT NONE
C
INTEGER I,ILASCH,NLOOP
CHARACTER*(*) STRING
C
C FIRST EXECUTABLE STATEMENT ILASCH
NLOOP = LEN(STRING)
IF (NLOOP.EQ.0 .OR. STRING.EQ.'; ';) THEN
ILASCH = 0
RETURN
ENDIF
C
DO 100 I = NLOOP, 1, -1
ILASCH = I
IF (STRING(I:I) .NE. '; ';) RETURN
100 CONTINUE
C
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKNCF (MDIM, ICKWRK, RCKWRK, NCF)
C
C START PROLOGUE
C
C SUBROUTINE CKNCF (MDIM, ICKWRK, RCKWRK, NCF)
C Returns the elemental composition of the species
C
C INPUT
C MDIM - First dimension of the two-dimensional array NCF;
C MDIM must be equal to or greater than the number of
C elements, MM.
C Data type - integer scalar
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C NCF - Matrix of the elemental composition of the species;
C NCF(M,K) is the number of atoms of the Mth element
C in the Kth species.
C Data type - integer array
C Dimension NCF(MDIM,*) exactly MDIM (at least MM,
C the total number of elements in the problem) for
C the first dimension and at least KK, the total
C number of species, for the second.
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,J,K,M,MDIM,NCF
DOUBLE PRECISION
1 RCKWRK
DIMENSION ICKWRK(*), RCKWRK(*), NCF(MDIM,*)
C
DO 150 K = 1, NKK
J = IcNC + (K-1)*NMM
DO 150 M = 1, NMM
NCF(M,K) = ICKWRK(J + M - 1)
150 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKXNUM (LINE, NEXP, LOUT, NVAL, RVAL, KERR)
C
C START PROLOGUE
C
C SUBROUTINE CKXNUM (LINE, NEXP, LOUT, NVAL, RVAL, KERR)
C This subroutine is called to parse a character string, LINE,
C that is composed of several blank-delimited substrings.
C Each substring is expected to represent a number, which
C is converted to entries in the array of real numbers, RVAL(*).
C NEXP is the number of values expected, and NVAL is the
C number of values found. This allows format-free input of
C numerical data. For example:
C
C input: LINE = " 0.170E+14 0 47780.0"
C NEXP = 3, the number of values requested
C LOUT = 6, a logical unit number on which to write
C diagnostic messages.
C output: NVAL = 3, the number of values found
C RVAL(*) = 1.700E+13, 0.000E+00, 4.778E+04
C KERR = .FALSE.
C
C INPUT
C LINE - A character string.
C Data type - CHARACTER*80
C NEXP - Number of real values to be found in character string.
C If NEXP is negative, then IABS(NEXP) values are
C expected. However, it is not an error condition,
C if less values are found.
C Data type - integer scalar
C LOUT - Output unit for printed diagnostics.
C Data type - integer scalar
C
C OUTPUT
C NVAL - Number of real values found in character string.
C Data type - integer scalar
C RVAL - Array of real values found.
C Data type - real array
C Dimension RVAL(*) at least NEXP
C KERR - Error flag; syntax or dimensioning error results
C in KERR = .TRUE.
C Data type - logical
C
C END PROLOGUE
C
C A ';!'; will comment out a line, or remainder of the line.
C
IMPLICIT NONE
C
INTEGER IERR,ILEN,IPPLEN,LOUT,N,NEXP,NVAL
DOUBLE PRECISION
1 RTEMP,RVAL
CHARACTER LINE*(*), ITEMP*80
DIMENSION RVAL(*), RTEMP(80)
LOGICAL KERR
C
C----------Find Comment String (! signifies comment)
C
ILEN = IPPLEN(LINE)
NVAL = 0
KERR = .FALSE.
C
IF (ILEN .LE. 0) RETURN
IF (ILEN .GT. 80) THEN
WRITE (LOUT,*) '; Error in CKXNUM...line length > 80 ';
WRITE (LOUT,';(A)';) LINE
KERR = .TRUE.
RETURN
ENDIF
C
ITEMP = LINE(:ILEN)
IF (NEXP .LT. 0) THEN
CALL IPPARR (ITEMP, -1, NEXP, RTEMP, NVAL, IERR, LOUT)
ELSE
CALL IPPARR (ITEMP, -1, -NEXP, RTEMP, NVAL, IERR, LOUT)
IF (IERR .EQ. 1) THEN
WRITE (LOUT, *) '; Syntax errors in CKXNUM...';
WRITE (LOUT,';(A)';) LINE
KERR = .TRUE.
ELSEIF (NVAL .NE. NEXP) THEN
WRITE (LOUT,*) '; Error in CKXNUM...';
WRITE (LOUT,';(A)';) LINE
KERR = .TRUE.
WRITE (LOUT,*) NEXP,'; values expected, ';,
1 NVAL,'; values found.';
ENDIF
ENDIF
IF (NVAL .LE. IABS(NEXP)) THEN
DO 20 N = 1, NVAL
RVAL(N) = RTEMP(N)
20 CONTINUE
ENDIF
C
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE IPPARR (STRING,ICARD,NEXPEC,RVAL,NFOUND,IERR,LOUT)
C BEGIN PROLOGUE IPPARR
C REFER TO IPGETR
C DATE WRITTEN 850625 (YYMMDD)
C REVISION DATE 851625 (YYMMDD)
C CATEGORY NO. J3.,J4.,M2.
C KEYWORDS PARSE
C AUTHOR CLARK,G.L.,GROUP C-3 LOS ALAMOS NAT';L LAB
C PURPOSE Parses real variables from a character variable. Called
C by IPGETR, the IOPAK routine used for interactive input.
C DESCRIPTION
C
C-----------------------------------------------------------------------
C IPPARR may be used for parsing an input record that contains real
C values, but was read into a character variable instead of directly
C into real variables.
C The following benefits are gained by this approach:
C - specification of only certain elements of the array is allowed,
C thus letting the others retain default values
C - variable numbers of values may be input in a record, up to a
C specified maximum
C - control remains with the calling program in case of an input
C error
C - diagnostics may be printed by IPPARR to indicate the nature
C of input errors
C
C The contents of STRING on input indicate which elements of RVAL
C are to be changed from their entry values, and values to which
C they should be changed on exit. Commas and blanks serve as
C delimiters, but multiple blanks are treated as a single delimeter.
C Thus, an input record such as:
C '; 1., 2,,4.e-5 , ,6.e-6';
C is interpreted as the following set of instructions by IPGETR:
C
C (1) set RVAL(1) = 1.0
C (2) set RVAL(2) = 2.0
C (3) leave RVAL(3) unchanged
C (4) set RVAL(4) = 4.0E-05
C (5) leave RVAL(5) unchanged
C (6) set RVAL(6) = 6.0E-06
C
C IPPARR will print diagnostics on the default output device, if
C desired.
C
C IPPARR is part of IOPAK, and is written in ANSI FORTRAN 77
C
C Examples:
C
C Assume RVAL = (0., 0., 0.) and NEXPEC = 3 on entry:
C
C input string RVAL on exit IERR NFOUND
C ------------- ---------------------- ---- ------
C '; 2.34e-3, 3 45.1'; (2.34E-03, 3.0, 45.1) 0 3
C ';2,,3.-5'; (2.0, 0.0, 3.0E-05) 0 3
C ';,1.4,0.028E4'; (0.0, 1.4, 280.0) 0 3
C ';1.0, 2.a4, 3.0'; (1.0, 0.0, 0.0) 1 1
C ';1.0'; (1.0, 0.0, 0.0) 2 1
C
C Assume RVAL = (0.,0.,0.,0.) and NEXPEC = -4 on entry:
C
C input string RVAL on exit IERR NFOUND
C ------------- ---------------------- ---- ------
C ';1.,2.'; (1.0, 2.0) 0 2
C ';,,3 4.0'; (0.0, 0.0, 3.0, 4.0) 0 4
C ';1,,3,,5.0'; (0.0, 0.0, 3.0, 0.0) 3 4
C
C arguments: (I=input,O=output)
C -----------------------------
C STRING (I) - the character string to be parsed.
C
C ICARD (I) - data statement number, and error processing flag
C < 0 : no error messages printed
C = 0 : print error messages, but not ICARD
C > 0 : print error messages, and ICARD
C
C NEXPEC (I) - number of real variables expected to be input. If
C < 0, the number is unknown, and any number of values
C between 0 and abs(nexpec) may be input. (see NFOUND)
C
C PROMPT (I) - prompting string, character type. A question
C mark will be added to form the prompt at the screen.
C
C RVAL (I,O) - the real value or values to be modified. On entry,
C the values are printed as defaults. The formal parameter
C corresponding to RVAL must be dimensioned at least NEXPEC
C in the calling program if NEXPEC > 1.
C
C NFOUND (O) - the number of real values represented in STRING,
C only in the case that there were as many or less than
C NEXPEC.
C
C IERR (O) - error flag:
C = 0 if no errors found
C = 1 syntax errors or illegal values found
C = 2 for too few values found (NFOUND < NEXPEC)
C = 3 for too many values found (NFOUND > NEXPEC)
C-----------------------------------------------------------------------
C
C REFERENCES (NONE)
C ROUTINES CALLED IFIRCH,ILASCH
C END PROLOGUE IPPARR
IMPLICIT NONE
C
INTEGER IBS,ICARD,IE,IERR,IES,ILASCH,LOUT,NC,NCH,NEXP,NEXPEC,
1 NFOUND
DOUBLE PRECISION
1 RVAL
CHARACTER STRING*(*), ITEMP*80
DIMENSION RVAL(*)
CHARACTER *8 FMT(16)
LOGICAL OKINCR
C
C FIRST EXECUTABLE STATEMENT IPPARR
IERR = 0
NFOUND = 0
NEXP = IABS(NEXPEC)
IE = ILASCH(STRING)
IF (IE .EQ. 0) GO TO 500
NC = 1
C
C--- OKINCR is a flag that indicates it';s OK to increment
C--- NFOUND, the index of the array into which the value
C--- should be read. It is set negative when a space follows
C--- a real value substring, to keep incrementing from
C--- occurring if a comma should be encountered before the
C--- next value.
C
OKINCR = .TRUE.
C
C--- begin overall loop on characters in string
C
100 CONTINUE
C
IF (STRING(NC:NC) .EQ. ';,';) THEN
IF (OKINCR) THEN
NFOUND = NFOUND + 1
ELSE
OKINCR = .TRUE.
ENDIF
C
GO TO 450
ENDIF
IF (STRING(NC:NC) .EQ. '; ';) GO TO 450
C
C--- first good character (non-delimeter) found - now find
C--- last good character
C
IBS = NC
160 CONTINUE
NC = NC + 1
IF (NC .GT. IE) GO TO 180
IF (STRING(NC:NC) .EQ. '; ';)THEN
OKINCR = .FALSE.
ELSEIF (STRING(NC:NC) .EQ. ';,';)THEN
OKINCR = .TRUE.
ELSE
GO TO 160
ENDIF
C
C--- end of substring found - read value into real array
C
180 CONTINUE
NFOUND = NFOUND + 1
IF (NFOUND .GT. NEXP) THEN
IERR = 3
GO TO 500
ENDIF
C
DATA FMT/ '; (E1.0)';, '; (E2.0)';, '; (E3.0)';, '; (E4.0)';,
1 '; (E5.0)';, '; (E6.0)';, '; (E7.0)';, '; (E8.0)';, '; (E9.0)';,
2 ';(E10.0)';, ';(E11.0)';, ';(E12.0)';, ';(E13.0)';, ';(E14.0)';,
3 ';(E15.0)';, ';(E16.0)';/
IES = NC - 1
NCH = IES - IBS + 1
ITEMP = '; ';
ITEMP = STRING(IBS:IES)
READ (ITEMP(:NCH), FMT(NCH), ERR = 400) RVAL(NFOUND)
GO TO 450
400 CONTINUE
IERR = 1
GO TO 510
450 CONTINUE
NC = NC + 1
IF (NC .LE. IE) GO TO 100
C
500 CONTINUE
IF (NEXPEC .GT. 0 .AND. NFOUND .LT. NEXP) IERR = 2
510 CONTINUE
C
IF (IERR .EQ. 0 .OR. ICARD .LT. 0) RETURN
IF (ICARD .NE. 0) WRITE(LOUT,';(A,I3)';)
1 ';!! ERROR IN DATA STATEMENT NUMBER';, ICARD
IF (IERR .EQ. 1) WRITE(LOUT,';(A)';)';SYNTAX ERROR, OR ILLEGAL VALUE';
IF (IERR .EQ. 2) WRITE(LOUT,';(A,I2, A, I2)';)
1 '; TOO FEW DATA ITEMS. NUMBER FOUND = '; , NFOUND,
2 '; NUMBER EXPECTED = ';, NEXPEC
IF (IERR .EQ. 3) WRITE(LOUT,';(A,I2)';)
1 '; TOO MANY DATA ITEMS. NUMBER EXPECTED = ';, NEXPEC
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKCPOR (T, ICKWRK, RCKWRK, CPOR)
C
C START PROLOGUE
C
C SUBROUTINE CKCPOR (T, ICKWRK, RCKWRK, CPOR)
C Returns the nondimensional specific heats at constant pressure;
C see Eq. (19).
C
C INPUT
C T - Temperature.
C cgs units - K
C Data type - real scalar
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C CPOR - Nondimensional specific heats at constant pressure
C for the species.
C cgs units - none
C Data type - real array
C Dimension CPOR(*) at least KK, the total number of
C species.
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K,L,N,NA1
DOUBLE PRECISION
1 CPOR,RCKWRK,T,TEMP,TN
DIMENSION ICKWRK(*), RCKWRK(*), TN(10), CPOR(*)
C
TN(1) = 1.D00
DO 150 N = 2, NCP
TN(N) = T**(N-1)
150 CONTINUE
C
DO 250 K = 1, NKK
L = 1
DO 220 N = 2, ICKWRK(IcNT + K - 1)-1
TEMP = RCKWRK(NcTT + (K-1)*MXTP + N - 1)
IF (T .GT. TEMP) L = L+1
220 CONTINUE
C
NA1 = NcAA + (L-1)*NCP2 + (K-1)*NCP2T
CPOR(K) = 0.D00
DO 250 N = 1, NCP
CPOR(K) = CPOR(K) + TN(N)*RCKWRK(NA1 + N - 1)
250 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKPY (RHO, T, Y, ICKWRK, RCKWRK, P)
C
C START PROLOGUE
C
C SUBROUTINE CKPY (RHO, T, Y, ICKWRK, RCKWRK, P)
C Returns the pressure of the gas mixture given the mass density,
C temperature and mass fractions; see Eq. (*).
C
C INPUT
C RHO - Mass density.
C cgs units - gm/cm**3
C Data type - real scalar
C T - Temperature.
C cgs units - K
C Data type - real scalar
C Y - Mass fractions of the species.
C cgs units - none
C Data type - real array
C Dimension Y(*) at least KK, the total number of
C species.
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C P - Pressure.
C cgs units - dynes/cm**2
C Data type - real scalar
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K
DOUBLE PRECISION
1 Y,RCKWRK,SUMYOW,P,RHO,T
DIMENSION Y(*), ICKWRK(*), RCKWRK(*)
C
SUMYOW = 0.0
DO 150 K = 1, NKK
SUMYOW = SUMYOW + Y(K)/RCKWRK(NcWT + K - 1)
150 CONTINUE
P = RHO * RCKWRK(NcRU) * T * SUMYOW
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKUMS (T, ICKWRK, RCKWRK, UMS)
C
C START PROLOGUE
C
C SUBROUTINE CKUMS (T, ICKWRK, RCKWRK, UMS)
C Returns the internal energies in mass units; see Eq. (30).
C
C INPUT
C T - Temperature.
C cgs units - K
C Data type - real scalar
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C UMS - Internal energies in mass units for the species.
C cgs units - ergs/gm
C Data type - real array
C Dimension UMS(*) at least KK, the total number of
C species.
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK(*),K
DOUBLE PRECISION
1 T, RCKWRK(*), UMS(*), RUT
C
CALL CKHMS (T, ICKWRK, RCKWRK, UMS)
RUT = T*RCKWRK(NcRU)
DO 150 K = 1, NKK
UMS(K) = UMS(K) - RUT/RCKWRK(NcWT+K-1)
150 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKCVBS (T, Y, ICKWRK, RCKWRK, CVBMS)
C
C START PROLOGUE
C
C SUBROUTINE CKCVBS (T, Y, ICKWRK, RCKWRK, CVBMS)
C Returns the mean specific heat at constant volume in mass units;
C see Eq. (36).
C
C INPUT
C T - Temperature.
C cgs units - K
C Data type - real scalar
C Y - Mass fractions of the species.
C cgs units - none
C Data type - real array
C Dimension Y(*) at least KK, the total number of
C species.
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C CVBMS - Mean specific heat at constant volume in mass units.
C cgs units - ergs/(gm*K)
C Data type - real scalar
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK(*),K
DOUBLE PRECISION
1 T,RCKWRK(*), Y(*), CVBMS
C
CALL CKCVMS (T, ICKWRK, RCKWRK, RCKWRK(NcK1))
C
CVBMS = 0.0
DO 100 K = 1, NKK
CVBMS = CVBMS + Y(K)*RCKWRK(NcK1 + K - 1)
100 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKCVMS (T, ICKWRK, RCKWRK, CVMS)
C
C START PROLOGUE
C
C SUBROUTINE CKCVMS (T, ICKWRK, RCKWRK, CVMS)
C Returns the specific heats at constant volume in mass units;
C see Eq. (29).
C
C INPUT
C T - Temperature.
C cgs units - K
C Data type - real scalar
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C CVMS - Specific heats at constant volume in mass units
C for the species.
C cgs units - ergs/(gm*K)
C Data type - real array
C Dimension CVMS(*) at least KK, the total number of
C species.
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK(*),K
DOUBLE PRECISION
1 T,RCKWRK(*), CVMS(*)
C
CALL CKCPMS (T, ICKWRK, RCKWRK, CVMS)
C
DO 150 K = 1, NKK
CVMS(K) = CVMS(K) - RCKWRK(NcRU) / RCKWRK(NcWT + K - 1)
150 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKHORT (T, ICKWRK, RCKWRK, HORT)
C
C START PROLOGUE
C
C SUBROUTINE CKHORT (T, ICKWRK, RCKWRK, HORT)
C Returns the nondimensional enthalpies; see Eq. (20).
C
C INPUT
C T - Temperature.
C cgs units - K
C Data type - real scalar
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C HORT - Nondimensional enthalpies for the species.
C cgs units - none
C Data type - real array
C Dimension HORT(*) at least KK, the total number of
C species.
C
C END PROLOGUE
C
IMPLICIT NONE
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK,K,L,N,NA1
DOUBLE PRECISION
x HORT,RCKWRK,SUM,T,TEMP,TN
DIMENSION TN(10), HORT(*), ICKWRK(*), RCKWRK(*)
C
DO 150 N = 1, NCP
TN(N) = T**(N-1)/N
150 CONTINUE
C
DO 250 K = 1, NKK
L = 1
DO 220 N = 2, ICKWRK(IcNT + K - 1)-1
TEMP = RCKWRK(NcTT + (K-1)*MXTP + N - 1)
IF (T .GT. TEMP) L = L+1
220 CONTINUE
C
NA1 = NcAA + (L-1)*NCP2 + (K-1)*NCP2T
SUM = 0.0
DO 225 N = 1, NCP
SUM = SUM + TN(N)*RCKWRK(NA1 + N - 1)
225 CONTINUE
HORT(K) = SUM + RCKWRK(NA1 + NCP1 - 1)/T
250 CONTINUE
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKAWT (ICKWRK, RCKWRK, AWT)
C
C START PROLOGUE
C
C SUBROUTINE CKAWT (ICKWRK, RCKWRK, AWT)
C Returns the atomic weights of the elements
C
C INPUT
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C AWT - Atomic weights of the elements.
C cgs units - gm/mole
C Data type - real array
C Dimension AWT(*) at least MM, the total number of
C elements in the problem.
C
C END PROLOGUE
C
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK(*)
DOUBLE PRECISION RCKWRK(*), AWT(*)
C
INTEGER M
C
DO 100 M = 1, NMM
AWT(M) = RCKWRK(NcAW + M - 1)
100 CONTINUE
C
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKWT (ICKWRK, RCKWRK, WT)
C
C START PROLOGUE
C
C SUBROUTINE CKWT (ICKWRK, RCKWRK, WT)
C Returns the molecular weights of the species
C
C INPUT
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C WT - Molecular weights of the species.
C cgs units - gm/mole
C Data type - real array
C Dimension WT(*) at least KK, the total number of
C species.
C
C END PROLOGUE
C
INCLUDE ';ckstrt.h';
C
INTEGER ICKWRK(*)
DOUBLE PRECISION RCKWRK(*), WT(*)
C
INTEGER K
C
DO 100 K = 1, NKK
WT(K) = RCKWRK(NcWT + K - 1)
100 CONTINUE
C
RETURN
END
C
C----------------------------------------------------------------------C
C C
SUBROUTINE CKATHM (NDIM1, NDIM2, ICKWRK, RCKWRK, MAXTP, NT, TMP,
1 A)
C
C START PROLOGUE
C
C SUBROUTINE CKATHM (NDIM1, NDIM2, ICKWRK, RCKWRK, MAXTP, NT, TMP,
C A)
C Returns the coefficients of the fits for thermodynamic properties
C of the species; see Eqns. (19)-(21).
C
C INPUT
C NDIM1 - First dimension of the three-dimensional array of
C thermodynamic fit coefficients, A; NDIM1 must be at
C least NCP2, the total number of coefficients for one
C temperature range.
C NDIM2 - Second dimension of the three-dimensionalarray of
C thermodynamic fit coefficients, A; NDIM2 must be at
C least MXPT-1, the total number of temperature ranges.
C ICKWRK - Array of integer workspace.
C Data type - integer array
C Dimension ICKWRK(*) at least LENIWK.
C RCKWRK - Array of real work space.
C Data type - real array
C Dimension RCKWRK(*) at least LENRWK.
C
C OUTPUT
C NT - Number of temperatures used for fitting coefficients of
C thermodynamic properties for the species.
C Data type - integer array
C Dimension NT(*) at least KK, the total number of
C species.
C TMP - Common temperatures dividing the thermodynamic fits for
C the species.
C cgs units - K
C Data type - real array
C Dimension TMP(MAXT,*) exactly MAXT for the first
C dimension (the maximum number of temperatures
C allowed for a species) , and at least KK for the
C second dimension (the total number of species)
C A - Three dimensional array of fit coefficients to the
C thermodynamic data for the species.
C The indicies in A(N,L,K) mean-
C N = 1,NN are polynomial coefficients in CP/R
C CP/R(K)=A(1,L,K) + A(2,L,K)*T + A(3,L,K)*T**2 + ...
C N = NN+1 is a6 in Eq. (20)
C N = NN+2 is a7 in Eq. (21)
C L = 1..MXTP-1 is for each temperature range.
C K is the species index
C Data type - real array
C Dimension A(NPCP2,NDIM2,*) exactly NPCP2 and MXTP-1
C for the first and second dimensions and at least
C KK for the third.
C
C END PROLOGUE
C
INCLUDE ';ckstrt.h';
C
INTEGER NDIM1, NDIM2, MAXTP, ICKWRK(*), NT(*)
DOUBLE PRECISION RCKWRK(*), TMP(MAXTP,*), A(NDIM1,NDIM2,*)
C
INTEGER K,L,NA1,M
DOUBLE PRECISION SUM
C
DO 150 K = 1, NKK
NT(K) = ICKWRK(IcNT + K - 1)
TMP(MXTP,K) = RCKWRK(NcTT + K*MXTP - 1)
C
DO 150 L = 1, MXTP-1
TMP(L,K) = RCKWRK(NcTT + (K-1)*MXTP + L - 1)
NA1 = NcAA + (L-1)*NCP2 + (K-1)*NCP2T
DO 150 M = 1, NCP2
A(M, L, K) = RCKWRK(NA1 + M - 1)
150 CONTINUE
C
RETURN
END
C
|
|
IP卡
狗仔卡
发表于 2006-5-10 22:30:11
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
楼主
