cgxbfgr.for
c
c
SUBROUTINE GXCHDA(MAXKK,MAXII,KK,II,SYMS,WT,ENTHCO,NDEC,
1 MAXPAT,NENTHP,NU,RGAS,PATM)
C-----------------------------------------------------------------------
C SUBROUTINE GXCHDA is activated by setting TMP1=GRND10 in the Q1
C file.
C If the variable CSG4 is blank then the data are read from the
C Q1 file in which case the data must be in columns numbered 3 or
C higher. If CSG4 is not blank, then the data is read from a file
C with name CSG4//'CHEM', eg. if CSG4=SCRS the file-name is
C SCRSCHEM.
C Data to be read must be placed between an upper line consisting of
C word CHEMBEGIN and a lower line consisting of the word CHEMEND.
C The keywords SIUNIT, SPECIES, THERM, REACT, RADIAT and TRANSP
C may be used to switch to SI units, to declare chemical species,
C to enter thermodynamic data, to enter reaction data, to enter
C radiation data (not active), and to enter transport-property
C data (not active).
C An example of a data-file is:
C
C ***********************************************************
C *
C * PHOENICS Chemistry-Data Input File for the Radian
C * Test-case (J.K. Worrell 02/09/94)
C *
C * Only the secondary reactions of the generalised SCRS
C * are active, the thermodynamic data is that given by
C * Adnani of Radian Corp. in his GROUND coding, and the
C * heats of formation are taken from Kanury and converted
C * from kcal/gm-mole to Joules/kg.
C *
C ***********************************************************
C *
C CHEMBEGIN
C *
C * Use SI units
C *
C SIUNIT
C *
C * Set up the species required
C *
C SPECIES
C O2 31.99880
C H2 2.01594
C CO 28.01055
C H2O 18.01534
C CO2 44.00995
C N2 28.01340
C END
C *
C * Set up the thermodynamic data for this case; all species
C * except CO2 and H2 have 2 "patch" approximations to the
C * temperature dependence of the specific heat. For CO2 there
C * is a single "patch" , and for H2 there are 3 "patches".
C *
C THERM
C O2 200. 700. 0.0
C 892. 1.41E-2 1.53E-4
C O2 600. 5000. 0.0
C 1423. -2.71E-2 8.36E-6 -9.96E3
C H2 200. 600. 0.0
C 1.19E4 10.9 -1.2E-2
C H2 500. 1450. 0.0
C 1.44E4 -0.345 9.55E-4
C H2 1350. 5000. 0.0
C 1.19E4 3.39 -4.23E-4
C CO 200. 700. -3.9466E6
C 1061. -0.177 3.65E-4
C CO 600. 5000. -3.9466E6
C 1157. 0.229 -5.28E-5 -4689.
C H2O 200. 1000. -13.4245E6
C 1786. 0.183 3.24E-4
C H2O 900. 5000. -13.4245E6
C 1371. 1.11 -1.43E-7
C CO2 200. 5000. -8.9417E6
C 1373. 0.24 -5.99E-5 -1.04E4
C N2 200. 700. 0.0
C 1051. -.123 2.77E-4
C N2 600. 5000. 0.0
C 918. 0.33 -7.E-5 -387.
C END
C *
C * Finally, set up the reactions;
C * 2CO + O2 --> CO2
C * and
C * 2H2 + O2 --> H2O
C *
C REACT CO -2 O2 -1 CO2 2
C REACT H2 -2 O2 -1 H2O 2
C CHEMEND
C
C-----------------------------------------------------------------------
C
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'pltcfile'
CHARACTER*(*) SYMS(MAXKK)
REAL WT(MAXKK),ENTHCO(NDEC,MAXPAT,MAXKK)
INTEGER NU(MAXKK,MAXII),NENTHP(MAXKK)
CHARACTER*256 NMSAVE
C
LOGICAL WORDIS,RDWERR
INTEGER H
COMMON/WORDI1/NWDS,NCHARS(20),NSEMI,H,NLINES
COMMON/WORDC1/WD(20),INLINE
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
CHARACTER WD*20,INLINE*120
C
COMMON/CHDATI/ NKK,NII,MAXP1,NDEC1
C
NAMSUB = 'GXCHDA'
C
IF(CSG4.NE.' ') THEN
NMSAVE=NMFIL(27)
NMFIL(27)=CSG4//'CHEM'
ENDIF
C
CALL SUB4(KK,0,II,0,MAXP1,MAXPAT,NDEC1,NDEC)
LSYM=LEN(SYMS(1))
DO 1 K=1,MAXKK
NENTHP(K)=0
WT(K)=0.0
SYMS(K)=' '
DO 1 N=1,MAXPAT
DO 1 J=1,NDEC
ENTHCO(J,N,K)=0.0
1 CONTINUE
CALL SUB2R(RGAS,8.314E7,PATM,1.01325E6)
C
CALL OPENQ1('CHEMBEGIN',IERR)
IF(IERR.EQ.0) THEN
CALL WRYT40('File opened for reading of chemical-data')
CALL WRITBL
IF(CSG4.EQ.' ') THEN
CALL WRIT40('>>>> Data read in from Q1 file <<<< ')
ELSE
CALL WRIT40('>>> Data read in from file '//CSG4//' <<<')
ENDIF
CALL WRITBL
4 CALL RDLNQ1('CHEMEND',IERR)
IF(IERR.EQ.0) THEN
C... If line is not blank, or CHEMEND or other problems echo it,
C weed out comments (ie. lines beginning with "*") and proceed
WRITE(LUPR1,*) INLINE
WRITE(LUPR3,*) INLINE
IF(WD(1)(1:1).EQ.'*') GO TO 4
IF(WORDIS(1,'SIUNIT')) THEN
C... SIUNIT: simply modify RGAS and PATM (more to be done if transport
C or radiation data is active?). The default is to CGS units,
C as is the case with CHEMKIN. This keyword may appear
C anywhere.
RGAS=RGAS*1.E-4
PATM=PATM*1.E-1
ELSE IF(WORDIS(1,'SPECIES')) THEN
C
C... SPECIES: this signals the start of the SPECIES block for declaring
C the species data. It MUST preceed all other blocks, and
C it must be terminated with an END. For each species, the
C name and molecular weight must be given. Any number of
C species may appear on a line so long as the name and
C molecular weight appear on the SAME line. Sample lines
C are;
C
C SPECIES
C O2 31.99880
C H2 2.01594
C CO 28.01055
C END
C
1004 CALL RDLNQ1('CHEMEND',IERR)
WRITE(LUPR1,*) INLINE
WRITE(LUPR3,*) INLINE
IF(WD(1)(1:1).EQ.'*') GO TO 1004
IF(IERR.EQ.0.AND..NOT.WORDIS(1,'END')) THEN
C... If not END or other exclusions enter block to decode line
IF(MOD(NWDS,2).NE.0) THEN
IERR=1
CALL WRIT40(' GXCHDA: SPECIES DATA FORMAT IS WRONG ')
ELSE IF(KK+NWDS/2.GT.MAXKK) THEN
IERR=2
CALL WRIT40(' GXCHDA: TRIED TO READ TOO MANY SPECIES ')
ELSE
C... No problems, so put the name (odd WorDs) into SYMS, and the mol.
C weight (even WorDs) into WT.
DO 1001 I=1,NWDS,2
KK=KK+1
SYMS(KK)(1:LSYM)=WD(I)
WT(KK)=RRDZZZ(I+1)
1001 CONTINUE
ENDIF
IF(IERR.NE.0.OR.RDWERR()) THEN
CALL WRYT40(' GXCHDA: failed to read species data ')
CALL SET_ERR(248,
* 'Error. GXCHDA: failed to read species data.',1)
C CALL EAROUT(1)
ENDIF
GO TO 1004
ENDIF
ELSE IF(WORDIS(1,'THERM')) THEN
C
C... THERM: this signals the start of a block for setting up
C THERModynamics data. This MUST appear after the SPECIES
C block for reasons which will be obvious. Provision is
C made for the specific heat which is a function of temperature
C to be approximated by up to MAXPAT "patches", ie. to have
C different coefficients for different temperature ranges.
C If temperature ranges overlap, then linear interpolation
C between the two different values is used. (See GXHMSK and
C GXCMSK). The data is given on 2 lines:
C
C Line 1 - 'species-name' 'Temp1' 'Temp2' 'Heat-of-form.'
C Line 2 - 'Coeff1' 'Coeff2' 'Coeff3' 'Coeff4'
C
C In this, Temp2 > Temp1 is essential, the Heat-of-form(ation)
C need only be given for the 1st "patch", and "trailing"
C Coeffs not specified are taken as 0. Sample lines are;
C
C THERM
C O2 200. 700. 0.0
C 892. 1.41E-2 1.53E-4
C O2 600. 5000. 0.0
C 1423. -2.71E-2 8.36E-6 -9.96E3
C END
C
2004 CALL RDLNQ1('CHEMEND',IERR)
WRITE(LUPR1,*) INLINE
WRITE(LUPR3,*) INLINE
IF(WD(1)(1:1).EQ.'*') GO TO 2004
IF(IERR.EQ.0.AND..NOT.WORDIS(1,'END')) THEN
IF(NWDS.LT.4-MIN(NENTHP(K),1)) THEN
CALL WRIT40(' GXCHDA: THERMO. DATA FORMAT IS WRONG ')
IERR=10
ELSE IF(KK.GT.0) THEN
C... Call GXMATC to find the current species (WD(1)) in SYMS and get K
C which indicates position.
CALL GXMATC(WD(1),SYMS,KK,K,IERR)
IF(IERR.NE.0) THEN
CALL WRIT40(
1 ' GXCHDA: (THERM) SPECIES UNRECOGNISED ')
ELSE
C... Now set N to the "patch" number which will be stored in NENTHP(K)
N=NENTHP(K)+1
IF(N.LE.MAXPAT) THEN
C... Having checked that we do not have too many "patches", put the
C incremented "patch" number in NENTHP(K), and read the T1, T2 and
C heat-of-form.
NENTHP(K)=N
DO 2010 J=1,4-MIN(N,2)
ENTHCO(J,N,K)=RRDZZZ(J+1)
2010 CONTINUE
IF(ENTHCO(2,N,K).LE.ENTHCO(1,N,K)) THEN
C... If T1 & T2 are incorrectly ordered; complain and fail
CALL WRIT40(
1 ' GXCHDA: (THERM) T2 < T1 IS DISALLOWED ')
IERR=16
ENDIF
C... Ensure consistency of heat-of-formation
IF(N.GT.1) ENTHCO(3,N,K)=ENTHCO(3,1,K)
ELSE
CALL WRIT40(
1 ' GXCHDA: (THERM) TOO MANY PATCHES ')
IERR=12
ENDIF
ENDIF
ELSE
CALL WRIT40(' GXCHDA: (THERM) NO SPECIES!! ')
IERR=13
ENDIF
IF(IERR.NE.0.OR.RDWERR()) THEN
CALL WRYT40(' GXCHDA: failed to read thermodyn. data ')
CALL SET_ERR(249,
* 'Error. GXCHDA: failed to read thermodyn. data.',1)
C CALL EAROUT(1)
ENDIF
C... Now read the second line for the coefficients
2014 CALL RDLNQ1('CHEMEND',IERR)
WRITE(LUPR1,*) INLINE
WRITE(LUPR3,*) INLINE
IF(WD(1)(1:1).EQ.'*') GO TO 2014
IF(IERR.EQ.0.AND..NOT.WORDIS(1,'END')) THEN
IF(NWDS.LE.NDEC-3) THEN
C... If there is enough space, read all the coefficients found
DO 2020 J=1,NWDS
ENTHCO(J+3,N,K)=RRDZZZ(J)
2020 CONTINUE
ELSE
CALL WRIT40(
1 ' GXCHDA: (THERM) NO SPACE FOR CP COEFF.S')
IERR=15
ENDIF
IF(IERR.NE.0.OR.RDWERR()) THEN
CALL WRYT40(
1 ' GXCHDA: failed to read thermodyn. data ')
CALL SET_ERR(250,
* 'Error. GXCHDA: failed to read thermodyn. data.',1)
C CALL EAROUT(1)
ENDIF
GO TO 2004
ENDIF
ENDIF
ELSE IF(WORDIS(1,'REACT')) THEN
C
C... REACT: this signals a line of REACTion data. Currently the
C format is that pairs of species name and the change in number
C moles are input for the species involved in each reaction.
C These lines MUST come after the SPECIES block. Sample lines
C are;
C
C REACT CH4 -2 O2 -1 CO 2 H2 4
C REACT CO -2 O2 -1 CO2 2
C REACT H2 -2 O2 -1 H2O 2
C
IF(IERR.EQ.0) THEN
IF(MOD(NWDS-1,2).NE.0) THEN
CALL WRIT40(
1 ' GXCHDA: (REACT) INCORRECT DATA FORMAT ')
IERR=21
ELSE IF(II+1.GT.MAXII) THEN
CALL WRIT40(
1 ' GXCHDA: (REACT) TOO MANY REACTIONS ')
IERR=22
ELSE
II=II+1
DO 3010 I=2,NWDS,2
CALL GXMATC(WD(I),SYMS,KK,K,IERR)
IF(IERR.NE.0) THEN
CALL WRIT40(
1 ' GXCHDA: (REACT) SPECIES UNRECOGNISED ')
GO TO 3020
ELSE
NU(K,II)=IRDZZZ(I+1)
ENDIF
3010 CONTINUE
3020 CONTINUE
ENDIF
ENDIF
ELSE IF(WORDIS(1,'RADIAT')) THEN
C
C... RADIAT: this block is for the entry of RADIATion data
C
CALL WRIT40(' GXCHDA: (RADIAT) NOT IMPLEMENTED ')
ELSE IF(WORDIS(1,'TRANSP')) THEN
C
C... TRANSP: this block is for the entry of TRANSPort data
C
CALL WRIT40(' GXCHDA: (TRANSP) NOT IMPLEMENTED ')
ELSE
CALL WRITBL
CALL WRITST
CALL WRIT40(' GXCHDA: INCORRECT DATA-TYPE KEYWORD ')
CALL WRIT40(' (ONLY SPECIES, THERM, REACT, RADIAT, ')
CALL WRIT40(' SIUNIT & END ARE ALLOWED) ')
CALL WRYT40(' GXCHDA: INCORRECT DATA-TYPE KEYWORD ')
CALL WRYT40(' (ONLY SPECIES, THERM, REACT, RADIAT, ')
CALL WRYT40(' SIUNIT & END ARE ALLOWED) ')
CALL WRITST
CALL WRITBL
CALL SET_ERR(251, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
IF(RDWERR())
1 CALL WRYT40(' GXCHDA: failed to read data from file ')
GO TO 4
ENDIF
ENDIF
C
C... Store NKK and NII in COMMON
C
CALL SUB2(NKK,KK,NII,II)
C
C... Now fiddle around with COEFF(3,,), ie. with dHform, the enthalpy
C of formation to ensure consistency. note, that at this point we
C also check that thermodynamic data has been set for ALL species.
C
DO 4010 K=1,KK
IF(NENTHP(K).GT.0) THEN
CALL GXFIXE(298.0,NENTHP(K),ENTHCO(1,1,K),ent)
ELSE
CALL WRIT40(' GXCHDA: NO DATA THERMODYNAMIC DATA FOR ')
CALL WRIT1I(' SPECIES',K)
ENDIF
4010 CONTINUE
C
IF(CSG4.NE.' ') THEN
NMFIL(27)=NMSAVE
CALL WRIT40('>>>> End of data input from '//CSG4//' <<<< ')
ELSE
CALL WRIT40('>>>>>> End of data input from Q1 <<<<< ')
ENDIF
C
NAMSUB = 'gxchda'
END
C
C-----------------------------------------------------------------------
c
SUBROUTINE GXMATC(STRING,LIST,KK,K,IERR)
CHARACTER*(*) STRING,LIST(*)
C
C-----------------------------------------------------------------------
C GXMATC: Match a string (STRING) against the list (LIST) of length
C KK and return an index to it (K)
C-----------------------------------------------------------------------
C
IERR=0
DO 1 K1=1,KK
IF(STRING.EQ.LIST(K1)) THEN
K=K1
RETURN
ENDIF
1 CONTINUE
IERR=11
END
C
C-----------------------------------------------------------------------
c
SUBROUTINE GXHBMS(TEMP,Y,NP,COEFF,ENTH)
C
C-----------------------------------------------------------------------
C GXHBMS: this calculates the enthalpy per unit mass for a
C gas mixture with temperature TEMP and composition in mass-frac.s
C Y. The NP and COEFF arrays contain the thermodynamic data to be
C passed down to GXHMSK. The resultannt enthalpy is returned in
C ENTH.
C-----------------------------------------------------------------------
C
COMMON/CHDATI/ NKK,NII,MAXP1,NDEC1
REAL Y(*),COEFF(NDEC1,MAXP1,*)
INTEGER NP(*)
ENTH=0.0
DO 10 K=1,NKK
CALL GXHMSK(TEMP,NP(K),COEFF(1,1,K),ENTHK)
ENTH=ENTH+Y(K)*ENTHK
10 CONTINUE
END
C
C-----------------------------------------------------------------------
c
SUBROUTINE GXHMSK(TEMP,NP,COEFF,ENTH)
C
C-----------------------------------------------------------------------
C GXHMSK: this calculates the enthalpy per unit mass for a single
C species. The NP and COEFF arrays contain the thermodynamic data to
C be used. So far only the the model used by Radian has been
C implemented; ie.
C
C Cp = Cp0 + Cp1*TEMP + Cp2*TEMP**2 + Cp3/SQRT(TEMP)
C so
C H = dHform' + Cp0*TEMP + 0.5*Cp1*TEMP**2 + 0.3333*Cp2*TEMP**3
C - 2.0*Cp3*SQRT(TEMP)
C
C Note the optional use of multiple overlapping "patch"es (as in
C the Radian model) and the interpolation within the overlap regions
C between the 2 active "patch"es (TEMP lying in more than 2 "patch"es
C is disallowed and causes failure). The resultant enthalpy is
C returned in ENTH.
C-----------------------------------------------------------------------
C
COMMON/CHDATI/ NKK,NII,MAXP1,NDEC1
REAL COEFF(NDEC1,MAXP1)
C... Note the limitted use of D.P.
DOUBLE PRECISION ENT(2)
INTEGER IUSE(2)
LOGICAL QNE
C
ISUM=0
DO 100 N=1,NP
T1=COEFF(1,N)
T2=COEFF(2,N)
IF(TEMP.GE.T1.AND.TEMP.LE.T2) THEN
ISUM=ISUM+1
IF(ISUM.GT.2) GOTO 9001
IUSE(ISUM)=N
ENT(ISUM)=COEFF(6,N)/3.
DO 20 I=5,3,-1
ENT(ISUM)=TEMP*ENT(ISUM)+COEFF(I,N)/MAX(I-3,1)
20 CONTINUE
CO7=COEFF(7,N)
IF(QNE(CO7,0.0)) ENT(ISUM)=ENT(ISUM)-2.*CO7*SQRT(TEMP)
ENDIF
100 CONTINUE
IF(ISUM.EQ.1) THEN
ENTH=ENT(1)
ELSE IF(ISUM.EQ.2) THEN
TA=MAX(COEFF(1,IUSE(1)),COEFF(1,IUSE(2)))
TB=MIN(COEFF(2,IUSE(1)),COEFF(2,IUSE(2)))
A=(TEMP-TA)/(TB-TA)
ENTH=A*ENT(1)+(1.0-A)*ENT(2)
ENDIF
RETURN
9001 CALL WRIT40(' GXHMSK: TEMP. FALLS IN TOO MANY PATCHES')
END
C
C-----------------------------------------------------------------------
c
SUBROUTINE GXCPBS(TEMP,Y,NP,COEFF,CP)
C
C-----------------------------------------------------------------------
C GXCPBS: this calculates the specific heat-capacity per unit mass for
C a gas mixture with temperature TEMP and composition in mass-frac.s
C Y. The NP and COEFF arrays contain the thermodynamic data to be
C passed down to GXCMSK. The resultannt specific-heat is returned in
C CP.
C-----------------------------------------------------------------------
C
COMMON/CHDATI/ NKK,NII,MAXP1,NDEC1
REAL Y(*),COEFF(NDEC1,MAXP1,*)
INTEGER NP(*)
CP=0.0
DO 10 K=1,NKK
CALL GXCMSK(TEMP,NP(K),COEFF(1,1,K),CPK)
CP=CP+Y(K)*CPK
10 CONTINUE
END
C
C-----------------------------------------------------------------------
c
SUBROUTINE GXCMSK(TEMP,NP,COEFF,CP)
C
C-----------------------------------------------------------------------
C GXCMSK: this calculates the specific heat-capacity per unit mass for
C a single species. The NP and COEFF arrays contain the thermodynamic
C data to be used. So far only the the model used by Radian has been
C implemented; ie.
C
C Cp = Cp0 + Cp1*TEMP + Cp2*TEMP**2 + Cp3/SQRT(TEMP)
C
C Note the optional use of multiple overlapping "patch"es (as in
C the Radian model) and the interpolation within the overlap regions
C between the 2 active "patch"es (TEMP lying in more than 2 "patch"es
C is disallowed and causes failure). The resultant enthalpy is
C returned in CP.
C-----------------------------------------------------------------------
C
COMMON/CHDATI/ NKK,NII,MAXP1,NDEC1
REAL COEFF(NDEC1,MAXP1)
C... Note the limitted use of D.P.
DOUBLE PRECISION CPA(2)
INTEGER IUSE(2)
LOGICAL QNE
C
ISUM=0
DO 100 N=1,NP
T1=COEFF(1,N)
T2=COEFF(2,N)
IF(TEMP.GE.T1.AND.TEMP.LE.T2) THEN
ISUM=ISUM+1
IF(ISUM.GT.2) GOTO 9001
IUSE(ISUM)=N
CPA(ISUM)=COEFF(6,N)
DO 20 I=5,4,-1
CPA(ISUM)=TEMP*CPA(ISUM)+COEFF(I,N)
20 CONTINUE
CO7=COEFF(7,N)
IF(QNE(CO7,0.0)) CPA(ISUM)=CPA(ISUM)+CO7/SQRT(TEMP)
ENDIF
100 CONTINUE
IF(ISUM.EQ.1) THEN
CP=CPA(1)
ELSE IF(ISUM.EQ.2) THEN
TA=MAX(COEFF(1,IUSE(1)),COEFF(1,IUSE(2)))
TB=MIN(COEFF(2,IUSE(1)),COEFF(2,IUSE(2)))
A=(TEMP-TA)/(TB-TA)
CP=A*CPA(1)+(1.0-A)*CPA(2)
ENDIF
RETURN
9001 CALL WRIT40(' GXCMSK: TEMP. FALLS IN TOO MANY PATCHES')
END
C
C-----------------------------------------------------------------------
c
SUBROUTINE GXFIXE(TEMP,NP,COEFF,ENTH)
C
C-----------------------------------------------------------------------
C GXFIXE: this calculates the enthalpy per unit mass for a single
C species and then modifies the value of COEFF(3) for each "patch" in
C by subtracting off the value of H at temperature TEMP (ie. 298.K)
C in order that at TEMP the enthalpy is equal to the specified enthalpy
C of formation. Note that AFTER this call the COEFF(3) will no longer
C be the same in each "patch". The NP and COEFF arrays contain the
C thermodynamic data to be used. So far only the the model used by
C Radian has been implemented; ie.
C
C Cp = Cp0 + Cp1*TEMP + Cp2*TEMP**2 + Cp3/SQRT(TEMP)
C so
C H = dHform + Cp0*TEMP + 0.5*Cp1*TEMP**2 + 0.3333*Cp2*TEMP**3
C - 2.0*Cp3*SQRT(TEMP)
C
C We need;
C
C H(T) = dHform + integral[Cp(T).dT]
C
C with limits T and T0 on the integral, where T0 is TEMP in this
C routine. Also,
C
C H(T) = dHform
C
C So, COEFF(3) must be dHform' where
C
C dHform' = dHform - integral(Cp(T).dT)T=T0
C
C = 2*dHform - H(T0)
C
C This subroutine converts dHform to dHform' and puts it back in
C COEFF(3).
C
C Note the optional use of multiple overlapping "patch"es (as in
C the Radian model) and the interpolation within the overlap regions
C between the 2 active "patch"es (TEMP lying in more than 2 "patch"es
C is disallowed and causes failure). The resultant enthalpy is
C returned in ENTH.
C-----------------------------------------------------------------------
C
COMMON/CHDATI/ NKK,NII,MAXP1,NDEC1
REAL COEFF(NDEC1,MAXP1)
C... Note the limitted use of D.P.
DOUBLE PRECISION ENT
LOGICAL QNE
C
DO 100 N=1,NP
T1=COEFF(1,N)
T2=COEFF(2,N)
ENT=COEFF(6,N)/3.
DO 20 I=5,3,-1
ENT=TEMP*ENT+COEFF(I,N)/MAX(I-3,1)
20 CONTINUE
CO7=COEFF(7,N)
IF(QNE(CO7,0.0)) ENT=ENT-2.0*CO7*SQRT(TEMP)
COEFF(3,N)=2.*COEFF(3,N)-ENT
100 CONTINUE
END
c