cgxbfgr.for
cFile-name GXSCRS.htm
c
C**** SUBROUTINE GXSCRI is called from the following sections of
C GREX3: Group 1 Section 1 to create storage and aqquire data
C Group 9 Section 1;
C " " 10
C " " 14
C " " 13
C
C It makes extensive use of data transmitted from the Q1 file by
C means of the SPEDAT command.
C
C Its functions, and its activation by the PIL command SCRS, are
c described in the PHOENICS Encyclopaedia, under the heading:
c Extended Simple Chemically-Reacting System.
C
SUBROUTINE GXSCRI
C INCLUDE 'chmkin'
INCLUDE 'patcmn'
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
COMMON/CMOBTP/OBJTYP(20) /CMNOTP/NMOTYP,NMCHAM
CHARACTER*14 OBJTYP,COBJTP,OBNAME
PARAMETER (NISCRS=150, NRSCRS=250, NCSCRS=20)
C Dynamic storage for SCRS model via the arrays RGSC, IGSC & CGSC,
C which at a later stage will be absorbed into the PHOENICS F-array
C via the use of GXMAKE.
COMMON/SCRSR/RGSC(NRSCRS)/SCRSI/IGSC(NISCRS)/SCRSC/CGSC(NCSCRS)
LOGICAL GRN,NEZ,ERROR,INLOBJ
SAVE CFUEL,COXID,CFP1,CFP2,JYFU,JYOX,JYMIX,JYWT,JYMIXP,JYMIXM,
1 JYECO,JYELFU,JYELOX,JL0Y,JNU,JNEP,JL0K,JAWT,JNCF,
2 LBMIXF,LBMMWT,NREAC,NRSYS,NSPEC,NELEM,RGAS,
3 TEMPFU,TEMPOX
CHARACTER*4 CFUEL,COXID,CFP1,CFP2,CGSC,CTEMP*1
C
C... Initialisation in Group 1, Section 1
IF(IGR.EQ.1) THEN
IF(ISC.NE.1) RETURN
c
C.... get values of NSPEC, NREAC, NRSYS, NELEM
CALL GETSDI('SCRS','NSPEC',NSPEC)
CALL GETSDI('SCRS','NREAC',NREAC)
CALL GETSDI('SCRS','NRSYS',NRSYS)
CALL GETSDI('SCRS','NELEM',NELEM)
call writ4i('nspec ',nspec,'nreac ',nreac,'nrsys ',nrsys,
1 'nelem ',nelem)
CALL WRYT40('special data read by GXSCRI')
C.... allocate storage
c.... the purpose of the following sequence is....
JYFU = 1
JYOX = JYFU + NSPEC
JYMIX = JYOX + NSPEC
JYWT = JYMIX + NSPEC
JYMIXP= JYWT + NSPEC
JYMIXM= JYMIXP+ NSPEC
JYECO = JYMIXM+ NSPEC
JAWT = JYECO + 21*NSPEC
JYELFU= JAWT + NELEM
JYELOX= JYELFU+ NELEM
C
JL0Y=1
JNU=JL0Y + NSPEC
JNEP=JNU + NSPEC*NREAC
JL0K=JNEP+ NSPEC
JNCF=JL0K+ NSPEC
C
NRW=6*NSPEC+ 21* NSPEC + 3*NELEM
NIW=3*NSPEC+ NSPEC*NREAC + NELEM*NSPEC
NCW=NSPEC + NELEM
IF(GRNDNO(10,TMP1)) THEN
NRW=NRW + NSPEC*3*7
NIW=NIW + NSPEC
ENDIF
IF(NRSCRS.LT.NRW) THEN
CALL WRIT40(' NOT ENOUGH RGSC STORAGE IN GXSCRS: ')
CALL WRIT40(' INCREASE RGSC')
CALL WRIT2I(' CURRENT',NRSCRS,'REQUIRED',NRW)
CALL SET_ERR(272, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
IF(NISCRS.LT.NIW) THEN
CALL WRIT40(' NOT ENOUGH IGSC STORAGE IN GXSCRS: ')
CALL WRIT40(' INCREASE IGSC')
CALL WRIT2I(' CURRENT',NISCRS,'REQUIRED',NIW)
CALL SET_ERR(273, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
IF(NCSCRS.LT.NCW) THEN
CALL WRIT40(' NOT ENOUGH CGSC STORAGE IN GXSCRS: ')
CALL WRIT40(' INCREASE CGSC')
CALL WRIT2I(' CURRENT',NCSCRS,'REQUIRED',NCW)
CALL SET_ERR(274, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
C
LBMIXF=LBNAME('F')
LBMMWT=LBNAME('MMWT')
c
c.... the purpose of the following sequence is....
CALL GETSDC('SCRS','CFUEL' ,CFUEL)
CALL GETSDC('SCRS','COXID' ,COXID)
CALL GETSDC('SCRS','CFP1 ' ,CFP1 )
CALL GETSDC('SCRS','CFP2 ' ,CFP2 )
c
CALL GETSDR('SCRS','TEMPFU',TEMPFU)
CALL GETSDR('SCRS','TEMPOX',TEMPOX)
c
DO 10121 I=1,NSPEC
WRITE(CTEMP,'(I1)') I
CALL GETSDR('SCRS','YFU'//CTEMP,RGSC(I) )
CALL GETSDR('SCRS','YOX'//CTEMP,RGSC(I+NSPEC))
10121 CONTINUE
c.... could not the content of these arguments be less clumsily be
c conveyed via common?
CALL GXSCRS(0,NRSYS,NSPEC,NREAC,NELEM,TMP1,CGSC,LBMIXF,
1 CFUEL,COXID,CFP1,CFP2,RGAS,RGSC(JYWT),RGSC(JYMIX),
2 TMPD,RGSC(JYFU),RGSC(JYOX),RGSC(JYMIXP),
3 RGSC(JYMIXM),TEMPFU,TEMPOX,IGSC(JL0Y),IGSC(JNU),
4 IGSC(JNEP),IGSC(JL0K),IGSC(JNCF),RGSC(JYECO),
5 RGSC(JAWT),RGSC(JYELFU),RGSC(JYELOX),
6 LBMMWT,TEMP1,H1,DEN1,0)
C
ELSEIF(IGR.EQ.9) THEN
C... Calculate density
IF(ISC.EQ.1) THEN
CALL GXSCRS(-1,NRSYS,NSPEC,NREAC,NELEM,TMP1,CGSC,LBMIXF,
1 CFUEL,COXID,CFP1,CFP2,RGAS,RGSC(JYWT),RGSC(JYMIX),
2 TMPD,RGSC(JYFU),RGSC(JYOX),RGSC(JYMIXP),
3 RGSC(JYMIXM),TEMPFU,TEMPOX,IGSC(JL0Y),IGSC(JNU),
4 IGSC(JNEP),IGSC(JL0K),IGSC(JNCF),RGSC(JYECO),
5 RGSC(JAWT),RGSC(JYELFU),RGSC(JYELOX),
6 LBMMWT,TEMP1,H1,DEN1,0)
C... Calculate temperature
ELSEIF(ISC.EQ.10) THEN
CALL GXSCRS(1,NRSYS,NSPEC,NREAC,NELEM,TMP1,CGSC,LBMIXF,
1 CFUEL,COXID,CFP1,CFP2,RGAS,RGSC(JYWT),RGSC(JYMIX),
2 TMPD,RGSC(JYFU),RGSC(JYOX),RGSC(JYMIXP),
3 RGSC(JYMIXM),TEMPFU,TEMPOX,IGSC(JL0Y),IGSC(JNU),
4 IGSC(JNEP),IGSC(JL0K),IGSC(JNCF),RGSC(JYECO),
5 RGSC(JAWT),RGSC(JYELFU),RGSC(JYELOX),
6 LBMMWT,TEMP1,H1,DEN1,0)
C... Calculate specific heat
ELSEIF(ISC.EQ.14) THEN
CALL GXSCRS(-2,NRSYS,NSPEC,NREAC,NELEM,TMP1,CGSC,LBMIXF,
1 CFUEL,COXID,CFP1,CFP2,RGAS,RGSC(JYWT),RGSC(JYMIX),
2 TMPD,RGSC(JYFU),RGSC(JYOX),RGSC(JYMIXP),
3 RGSC(JYMIXM),TEMPFU,TEMPOX,IGSC(JL0Y),IGSC(JNU),
4 IGSC(JNEP),IGSC(JL0K),IGSC(JNCF),RGSC(JYECO),
5 RGSC(JAWT),RGSC(JYELFU),RGSC(JYELOX),
6 LBMMWT,TEMP1,H1,DEN1,0)
ENDIF
C
ELSEIF(IGR.EQ.13) THEN
INLOBJ=.FALSE.
IF(NPATCH(1:1).EQ.'^') THEN
OBNAME=COBJTP(IPAT)
INLOBJ=(OBNAME(1:5).EQ.'INLET'.OR.
1 OBNAME(1:6).EQ.'BURNER')
ENDIF
C... Calculate reaction-rate source terms
IF(NPATCH(6:8).EQ.'RSO') THEN
CALL GXSCRS(3,NRSYS,NSPEC,NREAC,NELEM,TMP1,CGSC,LBMIXF,
1 CFUEL,COXID,CFP1,CFP2,RGAS,RGSC(JYWT),RGSC(JYMIX),
2 TMPD,RGSC(JYFU),RGSC(JYOX),RGSC(JYMIXP),
3 RGSC(JYMIXM),TEMPFU,TEMPOX,IGSC(JL0Y),IGSC(JNU),
4 IGSC(JNEP),IGSC(JL0K),IGSC(JNCF),RGSC(JYECO),
5 RGSC(JAWT),RGSC(JYELFU),RGSC(JYELOX),
6 LBMMWT,TEMP1,H1,DEN1,0)
C.... Calculate inlet enthalpy & density for inlet source terms
ELSEIF(NPATCH(5:5).EQ.'F'.OR.
1 (INLOBJ.AND.(OBNAME(6:7).EQ.'_F'.OR.OBNAME(7:8).EQ.'_F')))THEN
IADD=0
TEMPIN=TEMPFU
ELSEIF(NPATCH(5:5).EQ.'O'.OR.
1 (INLOBJ.AND.(OBNAME(6:7).EQ.'_O'.OR.OBNAME(7:8).EQ.'_O')))THEN
IADD=NSPEC
TEMPIN=TEMPOX
ENDIF
ENTH=0.0
IF(INDVAR.EQ.P1) THEN
IF(GRN(RHO1)) THEN
WTMIX=0.0
DO 13001 K=1,NSPEC
WTMIX=WTMIX+RGSC(K+IADD)/RGSC(JYWT+K-1)
13001 CONTINUE
WTMIX=1./WTMIX
RHOMIX=PRESS0*WTMIX/RGAS/(TEMPIN+TEMP0)
ELSE
RHOMIX=RHO1
ENDIF
CALL FN1(VAL,RHOMIX)
ENDIF
IF(BFC.AND.ISC.EQ.15) THEN
IF(INDVAR.EQ.P1.OR.(INDVAR.GE.U1.AND.INDVAR.LE.W1))
1 CALL GXBFC
ELSE
IF(INDVAR.EQ.P1) THEN
C... Cell type
IF(ITYPE.EQ.1) THEN
CALL FN1(VAL,0.0)
DO 13002 IVEL=U1,W1,2
CALL GETCOV(NPATCH,IVEL,VCO,VVALA)
IF(NEZ(VVALA)) CALL FN33(VAL,VVALA*VVALA)
13002 CONTINUE
CALL FN30(VAL)
CALL FN25(VAL,RHOMIX)
C... East (2) and West (3) type
ELSE IF(ITYPE.EQ.2.OR.ITYPE.EQ.3) THEN
CALL GETCOV(NPATCH,U1,VCO,VVALA)
C... North(4) and South(5) type
ELSE IF(ITYPE.EQ.4.OR.ITYPE.EQ.5) THEN
CALL GETCOV(NPATCH,V1,VCO,VVALA)
C... High (6) and Low (7) type
ELSE IF(ITYPE.EQ.6.OR.ITYPE.EQ.7) THEN
CALL GETCOV(NPATCH,W1,VCO,VVALA)
ENDIF
VVALA=ABS(VVALA)
CALL FN25(VAL,VVALA)
ENDIF
ENDIF
IF(INDVAR.EQ.H1) THEN
IF(GRNDNO(9,TMP1)) THEN
NTS=-1
CALL GXCKTC(RGSC(1+IADD),TEMPIN,TEMP0,ENTH,NTS,CP,ERROR)
C... Convert chemkin h1 from erg/g to J/kg
ENTH=ENTH*1.E-4
ELSE
CALL SET_ERR(275, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
CALL FN1(VAL,ENTH)
ENDIF
C
C... Calculate elemental composition & mole fractions for output
ELSEIF(IGR.EQ.19.AND.ISC.EQ.5) THEN
CALL GXSCRS(4,NRSYS,NSPEC,NREAC,NELEM,TMP1,CGSC,LBMIXF,
1 CFUEL,COXID,CFP1,CFP2,RGAS,RGSC(JYWT),RGSC(JYMIX),
2 TMPD,RGSC(JYFU),RGSC(JYOX),RGSC(JYMIXP),
3 RGSC(JYMIXM),TEMPFU,TEMPOX,IGSC(JL0Y),IGSC(JNU),
4 IGSC(JNEP),IGSC(JL0K),IGSC(JNCF),RGSC(JYECO),
5 RGSC(JAWT),RGSC(JYELFU),RGSC(JYELOX),
6 LBMMWT,TEMP1,H1,DEN1,0)
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXSCRS(ITASK,JSYS,KK,II,NE,TCP,SYMS,IMIXF,CFUEL,
2 COXID,CFP1,CFP2,RGAS,WT,YMIX,TMP,YFU,YOX,YMIXP,
3 YMIXM,TMPFU,TMPOX,L0Y,NU,NEP,L0KIND,NCF,ECO,
4 AWT,YELFU,YELOX,IMMWT,ITMP,IENT,IDEN,I1A)
C-----------------------------------------------------------------------
C This subroutine is called as follows:
C ITASK = 0 : initialisation via Group 1 Section 1 of GREX3.
C = 1 : calculation of composition, density and specific
C heat from GXTEMP when TMP1=GRND9 or GRND10.
C = -1 : set density from GXRHO when RHO1=GRND9;RHO1A=GRND9
C = -2 : set specific heat from GXSPHT when CP1=GRND7.
C = 3 : calculate and set chemical source terms of fuel and
C intermediate species from Group 13 of GREX3 ( only
C when using SCRS with finite-rate chemistry )
C = 4 : calculation and whf storage of elemental composition
C and mole fractions from Group 19 Section 5 of GREX3
C ( elements only when ELH, ELN, ELO or ELC are STOREd
C in the Q1 file; mole fractions only when MO2, MH2O
C etc are STOREd in the Q1 file ).
C Various 'systems' may be specified
C MOD(ISYS,2) = 0 : Primary step is active
C MOD(ISYS,3) = 0 : 1 or 2 secondary steps are active
C Both, either or none of the latter 2 conditions may be true.
C-----------------------------------------------------------------
C INCLUDE 'chmkin'
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
C
COMMON/GENI/NXNY,IFIL5(50),ISPH1,IFIL1(8)
COMMON/LMAIN2/LFIL1(2),DBGTZ,LFIL2(8)
LOGICAL STOCMP,STEP1,STEP2,TWOFPS,SOLENT,SOLTEM,STODEN,ERROR,
1 LFIL1,DBGTZ,LFIL2,DBGLOC,DBSCRS
INTEGER L0Y(KK),NU(KK,II),IR(3),NEP(KK),NCF(NE,KK),L0KIND(KK)
REAL YMIX(KK),YFU(KK),YOX(KK),WT(KK),ECO(21*KK),AWT(NE),
1 YMIXP(KK),YMIXM(KK),YELFU(NE),YELOX(NE)
CHARACTER*4 SYMS(KK+NE)
CHARACTER*4 CFUEL,COXID,CFP1,CFP2,LINE*40
LOGICAL FINRAT,SOLVFU,SOLVP1,SOLVP2,RFUMIN,SLD,DENCMP
INTEGER COGRN,VALGRN
C
SAVE OM1,FPM1,FPM2,OM12,OM22,PM1,PM2,IR,KFUEL,KOXID,KFP1,KFP2,
1 KPROD1,KPROD2,STOCMP,ENTHFU,ENTHOX,RGASL,L0DEN,IEQTSK,
2 L0CPM,IFULIN,LBYSUM,IREACT,IPRMOD,ISRMOD,FINRAT,LBFSQ,
3 IPCHEM,ISCHEM,PREXP,FUEXP,OXEXP,EACTP,PREXS1,FUEXS1,
4 OXEXS1,EACTS1,PREXS2,FUEXS2,OXEXS2,EACTS2
SAVE L0P1RS,L0S1RS,L0S2RS,L0PVAL,L0SVAL,SOLVFU,SOLVP1,MM,SOLVP2,
1 INDFU,INDOX,INDP1,INDP2,INDPR1,INDPR2,LBP1RS,LBS1RS,LBS2RS,
2 STEP1,STEP2,NREACT,CEBP,CEBS1,CEBS2,IPDF,LBABET,LBBBET,
3 DENCMP
C
IF(ITASK.EQ.1.or.itask.eq.-2) THEN
dbscrs=dbt
JZADD =ITWO(0,(IZ-1)*NXNY,PARAB)
cc ELSEIF(ITASK.EQ.-2) THEN
cc dbscrs=dbt
cc JZADD =ITWO(0,(IZ-1)*NXNY,PARAB)
ELSEIF(ITASK.EQ.-1) THEN
dbscrs=dbrho
JZADD =ITWO(0,(IZ-1)*NXNY,PARAB)
ELSE
dbscrs=.false.
ENDIF
dbgloc=dbscrs
ISYS=ABS(JSYS)
STEP1=MOD(ISYS,2).EQ.0
STEP2=MOD(ISYS,3).EQ.0
STODEN=IDEN.NE.0
L0ENT=0
IF(IENT.GT.0) THEN
SOLENT=SOLVE(IENT)
ELSE
SOLENT=.FALSE.
L0ENT=-IENT
ENDIF
if(dbgloc) call banner(1,'gxscrs',221096)
C--------------------------------------------------------------------
CHAPTER 1 : INITIALISATION (Group 1 Section 1) ITASK=0
C--------------------------------------------------------------------
IF(ITASK.EQ.0) THEN
C.... Retrieve satellite input data
c.... Reals
CALL GETSDR('SCRS','CEBP' ,CEBP )
CALL GETSDR('SCRS','CEBS1' ,CEBS1 )
CALL GETSDR('SCRS','CEBS2' ,CEBS2 )
CALL GETSDR('SCRS','PREXP' ,PREXP )
CALL GETSDR('SCRS','FUEXP' ,FUEXP )
CALL GETSDR('SCRS','OXEXP' ,OXEXP )
CALL GETSDR('SCRS','EACTP' ,EACTP )
CALL GETSDR('SCRS','PREXS1',PREXS1)
CALL GETSDR('SCRS','FUEXS1',FUEXS1)
CALL GETSDR('SCRS','OXEXS1',OXEXS1)
CALL GETSDR('SCRS','EACTS1',EACTS1)
CALL GETSDR('SCRS','PREXS2',PREXS2)
CALL GETSDR('SCRS','FUEXS2',FUEXS2)
CALL GETSDR('SCRS','OXEXS2',OXEXS2)
CALL GETSDR('SCRS','EACTS2',EACTS2)
c.... integers
CALL GETSDI('SCRS','IFULIN',IFULIN)
CALL GETSDI('SCRS','IPDF' ,IPDF )
CALL GETSDI('SCRS','IREACT',IREACT)
CALL GETSDI('SCRS','IPRMOD',IPRMOD)
CALL GETSDI('SCRS','ISRMOD',ISRMOD)
CALL GETSDI('SCRS','IPCHEM',IPCHEM)
CALL GETSDI('SCRS','ISCHEM',ISCHEM)
c.... logicals
DENCMP=.FALSE.
CALL GETSDL('SCRS','DENCMP',DENCMP)
C
C... Sort out the reactions
C
CALL WRITBL
CALL WRITST
CALL WRIT40('GREX3 file is GXSCRS.F of: 080497 ')
CALL WRIT40(' PHOENICS EXTENDED SCRS - Version 1.0 ')
CALL WRITST
CALL SUB3R(OM1,0.0,FPM1,0.0,FPM2,0.0)
CALL SUB4R(OM12,0.0,OM22,0.0,PM1,0.0,PM2,0.0)
C 3d stores although only 2*nx*ny required.
C MZ=ITWO(1,2,NZ.GT.1.AND..NOT.PARAB)
IF(STODEN) CALL GXMAKE(L0DEN,NXNY*NZ,'SDEN')
CALL GXMAKE(L0CPM,NXNY*NZ,'SCPM')
C
IF(GRNDNO(9,TCP)) THEN
C
C... Initialise various data using the CHEMKIN Interface
C LSG33 used to indicate equilibrium solution!!!!!
C
CALL GXCKIN(KK,KK1,II,MM,WT,AWT,RGAS,SYMS,NU,NCF,L0Y,LSG33)
call writ3i('Elements',MM,'Species ',KK,'Reaction',II)
IF(KK1.NE.KK) THEN
CALL WRIT40('GXSCRS: NO. SPECIES IN .CKM AND Q1 DO ')
CALL WRIT40('NOT MATCH ')
CALL SET_ERR(276, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
IEQTSK=0
ELSEIF(GRNDNO(10,TCP)) THEN
C
C... Initialise various data using the PHOENICS Chemistry system
C Note that it is ASSUMED that no more than 3 reactions are used
C (arg. 2), thay the thermodynamic data requires no. more than
C 7 elements per "patch" (arg. 8) and that no more than 3 "patch"es
C will be used for the thermodynamic data (arg. 9). RGAS and PRESS0
C are set by GXCHDA.
C
CALL GXCHDA(KK,3,KK1,II,SYMS,WT,ECO,7,3,NEP,NU,RGAS,PRESS0)
IF(KK1.NE.KK) THEN
CALL WRIT40('GXSCRS: NO. SPECIES IN .CKM AND Q1 DO ')
CALL WRIT40('NOT MATCH ')
CALL SET_ERR(277, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
ENDIF
RGASL=RGAS
IF(II.GT.3) THEN
CALL WRIT40('GXSCRS: NO MORE THAN 3 REACTIONS ARE ')
CALL WRIT40('ALLOWED ')
CALL SET_ERR(278, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
C
IF(II.EQ.0) THEN
IF(ISYS.NE.1) THEN
CALL WRIT40(' NO REACTIONS, BUT REACTION CODING IS ')
CALL WRIT40(' ACTIVATED (IE. ISYS .NE. 1) ')
CALL SET_ERR(279, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
ELSEIF(II.EQ.1) THEN
IF(STEP1) THEN
IF(STEP2) THEN
CALL WRIT40(' 1 REACTION, PRIMARY REACTION IS ACTIVE,')
CALL WRIT40(' SO NO SECONDARY REACTION IS PERMITTED ')
CALL SET_ERR(280, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
ELSEIF(STEP2) THEN
IF(STEP1) THEN
CALL WRIT40(' 1 REACTION, SECONDARY REACTION IS ')
CALL WRIT40(' ACTIVE SO NO PRIMARY REACTION IS ')
CALL WRIT40(' PERMITTED ')
CALL SET_ERR(281, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
ENDIF
ENDIF
IF(II.GT.0) THEN
CALL WRITST
CALL WRIT40(' REACTION STOICHOMETRIC COEFFICIENTS ')
DO 1011 I=1,II
CALL WRITBL
CALL WRIT1I('Reaction',I)
DO 1011 K=1,KK
LINE=SYMS(K)
WRITE(LINE(15:40),'(I3)') NU(K,I)
CALL WRIT40(LINE)
1011 CONTINUE
ENDIF
C
C... The array IR says which which reaction is which; IR(1) is
C index to the primary reaction, IR(2) to the 1st secondary
C reaction and IR(3) to the 2nd optional secondary reaction
C
DO 20 I=1,3
IR(I)=0
20 CONTINUE
CALL SUB2(KFUEL,0,KOXID,0)
IF(STEP1) THEN
IF(CFUEL.EQ.' ') THEN
CALL SET_ERR(302, 'Error. See result file.',1)
C CALL EAROUT(1)
ELSE
NUKI=0
DO 100 I=1,II
DO 110 K=1,KK
NUKI=NU(K,I)
IF(NUKI.GE.0) GO TO 110
IF(CFUEL.EQ.SYMS(K)) THEN
KFUEL=K
IMFU=NUKI
IF(IR(1).EQ.0) THEN
IR(1)=I
ELSEIF(IR(1).NE.I) THEN
GO TO 115
ENDIF
ELSEIF(COXID.EQ.SYMS(K)) THEN
KOXID=K
IMOX=NUKI
IR(1)=I
ENDIF
IF(KFUEL.GT.0.AND.KOXID.GT.0) GO TO 120
110 CONTINUE
100 CONTINUE
115 CALL WRIT40(' PRIMARY REACTION SPECIFIED, BUT FUEL ')
CALL WRIT40(' WAS NOT IDENTIFIED ')
CALL SET_ERR(282, 'Error. See result file.',1)
C CALL EAROUT(1)
120 I=IR(1)
CALL SUB2(KFP1,0,KFP2,0)
DO 130 K=1,KK
NUKI=NU(K,I)
IF(NUKI.GT.0) THEN
IF(CFP1.EQ.SYMS(K)) THEN
KFP1=K
IMFP1=NUKI
ELSEIF(CFP2.EQ.SYMS(K)) THEN
KFP2=K
IMFP2=NUKI
ENDIF
ENDIF
130 CONTINUE
ENDIF
C
C... Calculate the mass of OXIDant (OM1) consumed by burning
C unit mass of FUEL, and the corresponding masses of the
C intermediates created (FPM1 & FPM2)
C
OM1 =FLOAT(IMOX) * WT(KOXID) / (FLOAT(IMFU) * WT(KFUEL))
FPM1=FLOAT(IMFP1) * WT(KFP1) / (FLOAT(-IMFU)* WT(KFUEL))
FPM2=FLOAT(IMFP2) * WT(KFP2) / (FLOAT(-IMFU)* WT(KFUEL))
C define KPROD1 & KPROD2 as zero
CALL SUB2(KPROD1,0,KPROD2,0)
ENDIF
IF(STEP2) THEN
IF(II.EQ.2.AND.(CFP1.EQ.' '.OR.CFP2.EQ.' ')) THEN
CALL SET_ERR(283, 'Error. See result file.',1)
C CALL EAROUT(1)
ELSE
IF(.NOT.(STEP1.AND.II.EQ.2)) KFP2=0
CALL SUB2(KFP1,0,NUKI,0)
DO 200 I=1,II
IF(I.EQ.IR(1)) GO TO 200
DO 210 K=1,KK
NUKI=NU(K,I)
c!!!! what curious coding
IF(NUKI.EQ.0) GO TO 210
IF(NUKI.LT.0.AND.CFP1.EQ.SYMS(K)) THEN
IF(KFP1.EQ.0) THEN
KFP1=K
IMFP1=NUKI
ELSEIF(KFP1.NE.K) THEN
GO TO 205
ENDIF
IR(2)=I
IF(I.EQ.IR(3)) GO TO 205
ELSEIF(NUKI.LT.0.AND.CFP2.EQ.SYMS(K)) THEN
IF(KFP2.EQ.0) THEN
KFP2=K
IMFP2=NUKI
ELSEIF(KFP2.NE.K) THEN
GO TO 205
ENDIF
IR(3)=I
IF(I.EQ.IR(2)) GO TO 205
ENDIF
IF(IR(2).GT.0.AND.IR(3).GT.0) GO TO 220
210 CONTINUE
200 CONTINUE
205 IF(STEP1.AND.II.EQ.2.AND.IR(2).GT.0) GO TO 220
IF(.NOT.STEP1.AND.II.EQ.1.AND.IR(2).GT.0) GO TO 220
CALL WRIT40(' SECONDARY REACTION SPECIFIED, BUT ')
CALL WRIT40(' FUELS WERE NOT CORRECTLY IDENTIFIED ')
CALL SET_ERR(284, 'Error. See result file.',1)
C CALL EAROUT(1)
220 CALL SUB3(KPROD1,0,KPROD2,0,IMOX2,0)
IF(.NOT.STEP1) KOXID=0
MLAST=II
IF(.NOT.STEP1.AND.II.EQ.1) MLAST=2
DO 230 M=2,MLAST
I=IR(M)
DO 240 K=1,KK
NUKI=NU(K,I)
IF(NUKI.EQ.0) GO TO 240
IF(NUKI.GT.0.AND.I.EQ.IR(2)) THEN
IF(KPROD1.EQ.0) THEN
KPROD1=K
IMP1 =NUKI
ELSE
GO TO 235
ENDIF
ELSEIF(NUKI.GT.0.AND.I.EQ.IR(3)) THEN
IF(KPROD2.EQ.0) THEN
KPROD2=K
IMP2 =NUKI
ELSE
GO TO 235
ENDIF
ELSEIF(NUKI.LT.0.AND.I.EQ.IR(2).AND.K.NE.KFP1) THEN
IF(KOXID.EQ.0) THEN
KOXID=K
ELSEIF(KOXID.NE.K) THEN
K1=K
GO TO 237
ENDIF
IMOX1=NUKI
ELSEIF(NUKI.LT.0.AND.I.EQ.IR(3).AND.K.NE.KFP2) THEN
IF(KOXID.EQ.0) THEN
KOXID=K
ELSEIF(KOXID.NE.K) THEN
K1=K
GO TO 237
ENDIF
IMOX2=NUKI
ENDIF
IF(KPROD1.GT.0.AND.KPROD2.GT.0.AND.
1 IMOX1.NE.0.AND.(IMOX2.NE.0.OR.CFP2.EQ.' ')) GO TO 250
240 CONTINUE
230 CONTINUE
235 IF(STEP1.AND.II.EQ.2.AND.IMOX1.NE.0.AND
1 .KPROD1.GT.0) GO TO 250
IF(.NOT.STEP1.AND.II.EQ.1.AND.IMOX1.NE.0.AND
1 .KPROD1.GT.0) GO TO 250
CALL WRIT40(' TOO MANY PRODUCTS SPECIFIED FOR ')
CALL WRIT40(' SECONDARY REACTIONS ')
CALL SET_ERR(285, 'Error. See result file.',1)
C CALL EAROUT(1)
237 CALL WRIT40(' INCONSISTENT SPECIFICATION OF OXIDANT ')
CALL WRIT40(' SPECIES: CHECK THE REACTIONS ')
CALL WRIT2I(' KOXID1',KOXID,' KOXID2',K1)
CALL SET_ERR(286, 'Error. See result file.',1)
C CALL EAROUT(1)
250 CONTINUE
ENDIF
C
C... Calculate the mass of OXIDant (OM12 & OM22) consumed by
C burning unit mass of intermediates, and the corresponding
C masses of the products created (PM1 & PM2)
C
OM12=FLOAT(IMOX1)*WT(KOXID)/FLOAT(IMFP1)/WT(KFP1)
PM1 =FLOAT(IMP1)*WT(KPROD1)/FLOAT(-IMFP1)/WT(KFP1)
IF(II.GT.1) THEN
OM22=FLOAT(IMOX2)*WT(KOXID)/FLOAT(IMFP2)/WT(KFP2)
IF(KPROD2.GT.0) THEN
PM2=FLOAT(IMP2)*WT(KPROD2)/FLOAT(-IMFP2)/WT(KFP2)
ELSE
PM2=0
ENDIF
ELSE
OM22=0
PM2=0
ENDIF
ENDIF
NREACT=II
CALL WRITST
CALL WRIT40('SPECIES INDICES ')
CALL WRIT2I('KFUEL ',KFUEL,'KOXID ',KOXID)
CALL WRIT2I('KFP1 ',KFP1, 'KFP2 ',KFP2)
CALL WRIT2I('KPROD1 ',KPROD1,'KPROD2 ',KPROD2)
CALL WRIT40('STOICHOMETRIC MASS COEFFICIENTS ')
CALL WRIT2R('OM1 ',OM1,'FPM1 ',FPM1)
CALL WRIT2R('FPM2 ',FPM2,'OM12 ',OM12)
CALL WRIT1R('OM22 ',OM22)
CALL WRIT2R('PM1 ',PM1,'PM2 ',PM2)
C
C... If the first non-zero entry of YFU < 0, then inlet composition is
C in mole-frac.s so convert to mass-fracs
C
DO 252 K=1,KK
K1=K
IF(YFU(K).LT.0.0) GO TO 255
252 CONTINUE
C... If none are <0 then skip to 277 for mass-frac. input
GO TO 277
255 YFU(K1)=-YFU(K1)
DO 257 K=K1+1,KK
IF(YFU(K).LT.0.0) THEN
CALL WRIT40('MOLE-FRAC INPUT, BUT SOME FRAC.S ARE < 0')
CALL SET_ERR(287, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
257 CONTINUE
C
C... Print the initial compositions in mole-frac.s (if entry is in
C mole-frac.s)
C
CALL INCOND(-1,KK,SYMS,YFU,YOX,0.,0.,0.,0.,0.,0.)
WTMIX=0.0
DO 260 K=1,KK
WTMIX=WTMIX+YFU(K)*WT(K)
260 CONTINUE
DO 265 K=1,KK
YFU(K)=YFU(K)*WT(K)/WTMIX
265 CONTINUE
WTMIX=0.0
DO 270 K=1,KK
WTMIX=WTMIX+YOX(K)*WT(K)
270 CONTINUE
DO 275 K=1,KK
YOX(K)=YOX(K)*WT(K)/WTMIX
275 CONTINUE
C... Print the specified or derived mass-frac. inlet compositions
277 CONTINUE
ITIC=1
c IF(.NOT.SOLENT) THEN
C... Calculate the enthalpies in the inlet streams for later
C derivation of enthalpies and temperatures when enthalpy is
C not solved
IF(GRNDNO(9,TCP)) THEN
NTS=-1
CALL GXCKTC(YFU,TMPFU,TEMP0,ENTHFU,NTS,CP,ERROR)
CALL GXCKTC(YOX,TMPOX,TEMP0,ENTHOX,NTS,CP,ERROR)
C.... Convert chemkin enthalpy from erg/g to J/kg
ENTHFU=ENTHFU*1.E-4
ENTHOX=ENTHOX*1.E-4
ITIC=2
ELSEIF(GRNDNO(10,TCP)) THEN
CALL GXPHTC(YFU,TMPFU,TEMP0,ECO,NEP,ENTHFU,-1,CP,ERROR)
CALL GXPHTC(YOX,TMPOX,TEMP0,ECO,NEP,ENTHOX,-1,CP,ERROR)
ITIC=2
ENDIF
C... Calculate inlet densities for printout only
WTFU=0.0
WTOX=0.0
DO 2771 K=1,KK
WTFU=WTFU+YFU(K)/WT(K)
WTOX=WTOX+YOX(K)/WT(K)
2771 CONTINUE
WTFU=1./WTFU
WTOX=1./WTOX
RHOFU=PRESS0*WTFU/(RGASL*(TMPFU+TEMP0))
RHOOX=PRESS0*WTOX/(RGASL*(TMPOX+TEMP0))
IF(ERROR) THEN
CALL WRIT1I(' CELL',I)
CALL WRIT4R(' TMPFU',TMPFU ,' TMPOX',TMPOX,
1 ' ENTHFU',ENTHFU,' ENTHOX',ENTHOX)
CALL WRIT1R(' TEMP0',TEMP0)
CALL SET_ERR(288, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
c ENDIF
C
C... Print the final compositions in mass-frac.s with inlet
C temperatures and enthalpies
C
CALL INCOND(ITIC,KK,SYMS,YFU,YOX,TMPFU,TMPOX,ENTHFU,
1 ENTHOX,RHOFU,RHOOX)
LBYSUM=LBNAME('YSUM')
IF(IPDF.GT.0) LBFSQ=LBNAME('FSQ')
C.... Finite-Rate Chemistry Initialisation Sequence
LBCKF=NINT(CHSOB)
DO 281 K=1,KK
IF(STORE(LBCKF+K-1)) THEN
L0KIND(K)=LBCKF+K-1
ELSE
L0KIND(K)=0
ENDIF
281 CONTINUE
SOLVFU=SOLVE(L0KIND(KFUEL))
SOLVP1=SOLVE(L0KIND(KFP1))
SOLVP2=SOLVE(L0KIND(KFP2))
FINRAT=SOLVFU.OR.SOLVP1.OR.SOLVP2
IF(IREACT.EQ.2.OR.FINRAT) THEN
C... fuel and oxidiser
INDFU =L0KIND(KFUEL)
INDOX =L0KIND(KOXID)
C... intermediates
INDP1 =L0KIND(KFP1)
INDP2 =L0KIND(KFP2)
C... products
IF(KPROD1.NE.0) THEN
INDPR1=L0KIND(KPROD1)
ELSE
INDPR1=0
ENDIF
IF(KPROD2.NE.0) THEN
INDPR2=L0KIND(KPROD2)
ELSE
INDPR2=0
ENDIF
LBP1RS=LBNAME('P1RS')
LBS1RS=LBNAME('S1RS')
LBS2RS=LBNAME('S2RS')
IF(LBP1RS.EQ.0) CALL GXMAKE(L0P1RS,NXNY,'P1RS')
IF(LBS1RS.EQ.0) CALL GXMAKE(L0S1RS,NXNY,'S1RS')
IF(LBS2RS.EQ.0) CALL GXMAKE(L0S2RS,NXNY,'S2RS')
CALL GXMAKE(L0PVAL,NXNY,'PVAL')
CALL GXMAKE(L0SVAL,NXNY,'SVAL')
C... Calculate inlet elemental composition
CALL WRIT40(' INLET ELEMENTAL COMPOSITION BY MASS IS:')
CALL WRIT40(' ELEMENT FUEL OXIDANT ')
CALL WRIT40('MIXF 1 0 ')
DO 283 M=1,MM
YELFU(M)=0.0
LINE=SYMS(KK+M)
DO 282 K=1,KK
YELFU(M)=YELFU(M)+NCF(M,K)*AWT(M)*YFU(K)/WT(K)
YELOX(M)=YELOX(M)+NCF(M,K)*AWT(M)*YOX(K)/WT(K)
282 CONTINUE
WRITE(LINE(17:40),'(1P2E10.3)') YELFU(M),YELOX(M)
CALL WRIT40(LINE)
283 CONTINUE
CALL WRITST
ENDIF
C--------------------------------------------------------------------
CHAPTER 2 : DENSITY (Group 9 Section ?) ITASK=-1
C--------------------------------------------------------------------
ELSEIF(ITASK.EQ.-1) THEN
C
C... Copy the locally stored density into the PHOENICS F-array
C segment
C
if(dbgloc) then
call writ2i('igr ',igr,'isc ',isc)
call writ1i('iz ',iz)
call prn('3den',iden)
call prn('*den',-(l0den+jzadd))
endif
IF(L0DEN.GT.0) CALL FN0(IDEN,-(L0DEN+JZADD))
if(dbgloc) call prn('3den',iden)
C--------------------------------------------------------------------
CHAPTER 3 : COMPOSITION AND TEMPERATURE (Group 9 Section ?) ITASK=1
C--------------------------------------------------------------------
ELSEIF(ITASK.EQ.1) THEN
C
C... Calculate the various stages of the composition
C
TWOFPS=KFP2.GT.0
L0MIXF=L0F(IMIXF)
L0MMWT=0
IF(IMMWT.GT.0) L0MMWT=L0F(IMMWT)
IF(IENT.GT.0) L0ENT=L0F(IENT)
IF(IPDF.GT.0) L0FSQ=L0F(LBFSQ)
IF(LBABET.NE.0) L0ABET=L0F(LBABET)
IF(LBBBET.NE.0) L0BBET=L0F(LBBBET)
SOLTEM=ITMP.NE.0
IF(SOLTEM) L0TMP=L0F(ITMP)
IF(STODEN) L0P1=L0F(P1)
J0DEN=L0DEN+JZADD
J0CPM=L0CPM+JZADD
C IF(STOCMP) CALL GXCKIN(-KK,KK1,II,MM,WT,AWT,RGAS,SYMS,NU,
C 1 NCF,L0Y)
STOCMP=CHSOB.GT.0
IF(STOCMP) THEN
LBCKF=NINT(CHSOB)
STOCMP=.FALSE.
DO 300 K=1,KK
IF(STORE(LBCKF+K-1)) THEN
L0Y(K)=L0F(LBCKF+K-1)
STOCMP=.TRUE.
ELSE
L0Y(K)=0
ENDIF
300 CONTINUE
ENDIF
C
C... Determine whether only a single cell is to be calculated
C (for the CHEMKIN properties), or a whole-slab is to be
C calculated for all other calls
C
IF(I1A.LE.0) THEN
CALL SUB2(I1,1,I2,NXNY)
ELSE
CALL SUB2(I1,I1A,I2,I1A)
SOLTEM=.FALSE.
ENDIF
c IF(VPOR.NE.0) L0VPOR=L0F(VPOR)
cc IF(IPRPS.GT.0) L0PRPS = L0F(IPRPS)
IF(LBYSUM.NE.0) L0YSUM=L0F(LBYSUM)
C... START OF THE MAIN DO LOOP
DO 500 I=I1,I2
c IF(VPOR.NE.0) THEN
c IF(F(L0VPOR+I).LT.1.E-10) THEN
cc IF(IPRPS.GT.0) THEN
cc IF(F(L0PRPS+I).GE.SOLPRP) THEN
IF(SLD(I)) THEN
IF(SOLTEM) F(L0TMP+I)=1.E-10
IF(STODEN) F(J0DEN+I)=1.E-10
GO TO 500
ENDIF
cc ENDIF
FMIX=F(L0MIXF+I)
FMIX1=1.-FMIX
WTMIX=0.0
C... IREACT=0 Unburnt composition calculation
C =1 Fast-chemistry calculation
C =2 Finite-rate chemistry calculation
IF(IREACT.EQ.0) THEN
DO 410 K=1,KK
YMIX(K)=FMIX*YFU(K)+FMIX1*YOX(K)
410 CONTINUE
ELSE
IF(IREACT.EQ.1.AND.IPDF.EQ.1) THEN
C DOUBLE-DELTA PDF
GFSQ=AMAX1(F(L0FSQ+I),TINY)
CALL GXDPDF(GFSQ,FMIX,FPLUS,FMINUS,FALFA)
CALL GXCOMP(IREACT,KK,YFU,YOX,L0Y,FPLUS,STEP1,STEP2,
1 TWOFPS,SOLVP2,KFUEL,KOXID,KFP1,KFP2,KPROD1,KPROD2,
2 OM1,OM12,OM22,FPM1,FPM2,PM1,PM2,I,YMIXP)
CALL GXCOMP(IREACT,KK,YFU,YOX,L0Y,FMINUS,STEP1,STEP2,
1 TWOFPS,SOLVP2,KFUEL,KOXID,KFP1,KFP2,KPROD1,KPROD2,
2 OM1,OM12,OM22,FPM1,FPM2,PM1,PM2,I,YMIXM)
FALFA1=1.-FALFA
DO 411 K=1,KK
YMIX(K)=FALFA*YMIXP(K)+FALFA1*YMIXM(K)
411 CONTINUE
GFSQ=AMIN1((F(L0FSQ+I)+TINY),FMIX*FMIX1)
ELSE
C... Composition for finite-rate chemistry and no-pdf fast chemistry
CALL GXCOMP(IREACT,KK,YFU,YOX,L0Y,FMIX,STEP1,STEP2,
1 TWOFPS,SOLVP2,KFUEL,KOXID,KFP1,KFP2,KPROD1,KPROD2,
2 OM1,OM12,OM22,FPM1,FPM2,PM1,PM2,I,YMIX)
ENDIF
ENDIF
C... Ensure composition sums to unity
GYSUM=0.0
DO 4790 K=1,KK
4790 GYSUM=GYSUM+YMIX(K)
IF(LBYSUM.NE.0) F(L0YSUM+I)=GYSUM
DO 4791 K=1,KK
4791 YMIX(K)=YMIX(K)/GYSUM
C... Calculate mean molecular weight of the mixture
WTMIX=0.0
DO 430 K=1,KK
WTMIX=WTMIX+YMIX(K)/WT(K)
430 CONTINUE
WTMIX=1./WTMIX
C... Transfer composition in storage
IF(STOCMP) THEN
DO 450 K=1,KK
IF(L0Y(K).GT.0) F(L0Y(K)+I)=YMIX(K)
450 CONTINUE
ENDIF
IF(SOLTEM.OR.SOLENT.OR.L0ENT.GT.0.OR.I1.EQ.I2) THEN
C IF(SOLTEM) THEN
C TMP= F(L0TMP+I)
C ELSE
TMP=FMIX*TMPFU+FMIX1*TMPOX
C ENDIF
C
C... Either get the enthalpy from the F-array, if it is SOLVEd,
C or calculate it from the "adiabatic" SCRS model. If it is
C STOREd, then put the value into the F-array.
C
IF(SOLENT) THEN
ENTH=F(L0ENT+I)
ELSE
ENTH=FMIX*ENTHFU+FMIX1*ENTHOX
IF(L0ENT.GT.0) F(L0ENT+I)=ENTH
ENDIF
C
C... Calculate the temperature from the enthalpy using either
C the CHEMKIN Interface routine (GRND9) or the PHOENICS
C chemistry system (GRND10). Note the iteration limit, in
C the arg. list, is hard-wired to 5.
C
IF(GRNDNO(9,TCP)) THEN
C... Convert PHOENICS enthalpy from J/kg to egr/g for CHEMKIN
ENTH=ENTH*1.E4
NTS=5
CALL GXCKTC(YMIX,TMP,TEMP0,ENTH,NTS,CP,ERROR)
C... Convert chemkin cp from erg/g/K to J/kg/K
CP=CP*1.E-4
F(J0CPM+I)=CP
IF(LSG33) THEN
C
C... If LSG33 is true then do the equilibrium solution (fixed pressure
C and enthalpy). This returns the temperature, composition and
C molecular weight WTMIX.
C
CALL GXCKEQ(LU,IEQTSK,YMIX,TMP,TEMP0,PRESS0,WTMIX,
1 ERROR)
ENDIF
ELSEIF(GRNDNO(10,TCP)) THEN
CALL GXPHTC(YMIX,TMP,TEMP0,ECO,NEP,ENTH,5,CP,ERROR)
F(J0CPM+I)=CP
ENDIF
C IF(ERROR) THEN
C CALL WRIT40(' Temperature iteration not converged at:')
C CALL WRIT3I(' ISWEEP',ISWEEP,' IGR',IGR,
C 1 ' ISC',ISC)
C CALL WRIT2I(' CELL',I, ' IZ',IZSTEP)
C CALL WRIT2R(' TMP',TMP,' TEMP0',TEMP0)
C CALL WRIT3R(' ENTH',ENTH,' CP',CP,' FMIX',FMIX)
C DO 4301 K=1,KK
C WRITE(14,*) 'Y(K) = ',YMIX(K),' WT(K) = ',WT(K)
C 4301 CONTINUE
C CALL EAROUT(1)
C ENDIF
C
C... If temperature is STOREd (check the logic!) put the temperature
C into the F-array
C
IF(SOLTEM) F(L0TMP+I)=TMP
ENDIF
IF(L0MMWT.GT.0) F(L0MMWT+I)=WTMIX
C
C... If the density is stored, then put this into a temporary buffer
C for later use (this has NOT been tested for 2-phase cases!).
C Note that the solved pressure is omitted and only PRESS0 is used
C which is intended to increase stability but which is WRONG for
C compressible cases (especially detonations!!)
C
IF(STODEN) THEN
IF(DENCMP) THEN
RHOL=(F(L0P1+I)+PRESS0)*WTMIX/(RGASL*(TMP+TEMP0))
ELSE
RHOL=PRESS0*WTMIX/(RGASL*(TMP+TEMP0))
ENDIF
F(J0DEN+I)=RHOL
ENDIF
500 CONTINUE
if(dbgloc) then
call writ2i('igr ',igr,'isc ',isc)
call writ2i('isweep ',isweep,'iz ',iz)
if(soltem) call prn('*tmp',-l0tmp)
if(stoden) call prn('*den',-J0den)
call prn('*cpm',-J0cpm)
endif
C--------------------------------------------------------------------
CHAPTER 4 : MIXTURE SPECIFIC HEAT (Group 9 Section 14) ITASK=-2
C--------------------------------------------------------------------
ELSEIF(ITASK.EQ.-2) THEN
if(dbgloc) then
call writ2i('igr ',igr,'isc ',isc)
call writ1i('iz ',iz)
call prn('3cpm',lcp1)
call prn('*cpm',-(L0CPM+JZADD))
endif
CALL FN0(LCP1,-(L0CPM+JZADD))
IF(ISPH1.GT.0) CALL FN0(ISPH1,-(L0CPM+JZADD))
if(dbgloc) call prn('3cpm',lcp1)
C--------------------------------------------------------------------
CHAPTER 5 : FINITE-RATE REACTION SOURCES (Group 13 ) ITASK=2 & 3
C--------------------------------------------------------------------
ELSEIF(ITASK.EQ.3) THEN
C.... CHEMICAL SOURCE TERMS -----------------------GRND1
CALL SUB2(COGRN,ISC-1,VALGRN,ISC-12)
CALL SUB2(L0KE,L0F(KE),L0EP,L0F(EP))
C.... Primary reaction chemical source terms
C CO=GRND2
C IPRMOD=0 RFUi = -CR*MIN(MFU,MO/s)*RHO*VOL*EP/K
C =1 = -CR*MIN(MFU,MO/s,B*MPR/[1+s])*RHO*VOL*EP/K
C =2 = -CR*MFU*R*(1-R)**RHO*VOL*EP/K
C where R=MFU/MFUU & MFUU=F*MFU,A+(1-F)*MFU,B
C =3 = -PREXP*(FUMC**B)*(OXMC**C)*EXP(-TA/T)*VOL
IF(NPATCH.EQ.'SCRSPRSO') THEN
C.....KFUEL
IF(INDVAR.EQ.INDFU) THEN
IF(LBP1RS.NE.0) L0P1RS=L0F(LBP1RS)
IF(COGRN.EQ.2) THEN
CALL SUB3(L0MFU,L0F(INDVAR),L0MOX,L0F(INDOX),L0CO,L0F(CO))
cc IF(IPRPS.NE.0) L0PRPS=L0F(IPRPS)
RSTOIC=1./OM1
IF(IPRMOD.EQ.1) THEN
BEBU=0.5
RPROD=BEBU/(1.+OM1)
CALL SUB2(L0FP1,L0F(INDP1),L0FP2,L0F(INDP2))
ENDIF
IF(IPRMOD.EQ.3.OR.IPCHEM.EQ.1) THEN
TACTP=EACTP/1.987
PRECNP=1.E3*PREXP*WT(KFUEL)/1.E3**(FUEXP+OXEXP)
ITEMP=TEMP1
IF(TEMP0.GT.0.0) THEN
CALL FN2(EASP1,ITEMP,TEMP0,1.0)
ITEMP=EASP1
ENDIF
CALL SUB2(L0RHO,L0F(DEN1),L0TEM,L0F(ITEMP))
ENDIF
IF(IPRMOD.EQ.2) L0MIXF=L0F(IMIXF)
IF(IPRMOD.EQ.3) THEN
CALL FN1(CO,PRECNP)
ELSE
CALL FN15(CO,EP,KE,0.0,CEBP)
ENDIF
RFUMIN=IPRMOD.LE.2.AND.IPCHEM.EQ.1
C
DO 515 IX = IXF,IXL
IPLUS=(IX-1)*NY
CDIR$ IVDEP
DO 510 IY = IYF,IYL
I = IY + IPLUS
FU = F(L0MFU+I)
OX = F(L0MOX+I)
cc IF(IPRPS.NE.0) THEN
cc IF(F(L0PRPS+I).GE.SOLPRP) THEN
IF(SLD(I)) THEN
F(L0CO+I) = 0.0
F(L0P1RS+I) = 0.0
F(L0PVAL+I) = 0.0
GO TO 510
ENDIF
cc ENDIF
IF(IPRMOD.EQ.0) THEN
F(L0CO+I)=F(L0CO+I)*AMIN1(FU,OX*RSTOIC)/(FU+TINY)
ELSEIF(IPRMOD.EQ.1) THEN
PR=F(L0FP1+I)+F(L0FP2+I)
F(L0CO+I)=F(L0CO+I)*AMIN1(FU,OX*RSTOIC,RPROD*PR)
1 /(FU+TINY)
ELSEIF(IPRMOD.EQ.2) THEN
FMIX=F(L0MIXF+I)
FUUB=FMIX*YFU(KFUEL)+(1-FMIX)*YOX(KFUEL)
RPV1=FU/(FUUB+TINY)
F(L0CO+I)=F(L0CO+I)*RPV1*(1.-RPV1)
ELSEIF(IPRMOD.EQ.3) THEN
DENGAS=F(L0DEN+I)
FUMC=(DENGAS*FU/WT(KFUEL))**FUEXP
OXMC=(DENGAS*OX/WT(KOXID))**OXEXP
EXPTRM=EXP(-TACTP/F(L0TEM+I))
F(L0CO+I)=F(L0CO+I)*FUMC*OXMC*EXPTRM/
1 (DENGAS*(FU+TINY))
ENDIF
F(L0P1RS+I) = - F(L0CO+I)*FU
GSOR=F(L0P1RS+I)
IF(RFUMIN) THEN
DENGAS=F(L0DEN+I)
FUMC=(DENGAS*FU/WT(KFUEL))**FUEXP
OXMC=(DENGAS*OX/WT(KOXID))**OXEXP
EXPTRM=EXP(-TACTP/F(L0TEM+I))
GSORAR=-PRECNP*FUMC*OXMC*EXPTRM/DENGAS
RFURAT=GSORAR/(GSOR+TINY)
IF(ABS(RFURAT).LT.1.0) THEN
F(L0P1RS+I)=GSORAR
F(L0CO+I)=-GSORAR/(FU+TINY)
ENDIF
ENDIF
IF(IFULIN.EQ.1) THEN
FUB = AMAX1(0.0, (AMIN1(FU,FU-OX*RSTOIC)) )
FUMFUB=AMAX1(TINY,FU-FUB)
F(L0CO+I) = -F(L0P1RS+I)/FUMFUB
F(L0PVAL+I) = FUB
GCO=F(L0CO+I)
GVAL=F(L0PVAL+I)
ELSE
F(L0PVAL+I) = 0.0
ENDIF
510 CONTINUE
515 CONTINUE
ELSEIF(VALGRN.GE.0.AND.VALGRN.LE.2) THEN
CALL FN0(VAL,-L0PVAL)
ENDIF
C KFP1
ELSEIF (INDVAR.EQ.INDP1) THEN
CALL FN2(VAL,-L0P1RS,0.0,-FPM1)
C KFP2
ELSEIF (INDVAR.EQ.INDP2) THEN
CALL FN2(VAL,-L0P1RS,0.0,-FPM2)
ENDIF
ENDIF
C
C.... Secondary reaction chemical source terms per unit mass
C COGRN =2 One secondary reaction OXLIM=MOX/s
C =5 Two secondary reactions OXLIM=MFU*MOX/(MFU*s+MFU2*s2)
C ISRMOD=0 RFU = -CR*MIN(MFU,OXLIM)*RHO*VOL*EP/K
C =1 = -CR*MIN(MFU,OXLIM,B*MPR[1+s])*RHO*VOL*EP/K
C =2 = -CR*MFU*R*(1-R)*RHO*VOL*EP/K
C where R=MFU/MFUU & MFUU=F*MFU,A+(1-F)*MFU,B
C =3 = -PREXP*(FUMC**B)*(OXMC**C)*EXP(-TA/T)*VOL
IF(NPATCH.EQ.'SCRSSRSO') THEN
C KFP1 & KFP2
IF(COGRN.LE.5) THEN
CALL SUB3(L0VAR,L0F(INDVAR),L0MOX,L0F(INDOX),L0CO,L0F(CO))
IF(INDVAR.EQ.INDP1) THEN
STOIC = OM12
RSTOIC=1./STOIC
IF(LBS1RS.NE.0) L0S1RS=L0F(LBS1RS)
L0SSOR = L0S1RS
CREACT = CEBS1
KWTFP = KFP1
IF(COGRN.EQ.5) THEN
STOIC2=OM22
L0MFU2=L0F(INDP2)
ENDIF
IF(ISRMOD.EQ.1) THEN
BEBU=0.5
RPROD=BEBU/(1.+OM12)
L0FPRD=L0F(INDPR1)
ENDIF
IF(ISRMOD.EQ.3.OR.ISCHEM.EQ.1) THEN
CALL SUB2R(PREXS,PREXS1,FUEXS,FUEXS1)
CALL SUB2R(OXEXS,OXEXS1,EACTS,EACTS1)
ENDIF
ELSE
STOIC = OM22
RSTOIC=1./STOIC
IF(LBS2RS.NE.0) L0S2RS=L0F(LBS2RS)
L0SSOR = L0S2RS
CREACT = CEBS2
KWTFP = KFP2
IF(COGRN.EQ.5) THEN
STOIC2=OM12
L0MFU2=L0F(INDP1)
ENDIF
IF(ISRMOD.EQ.1) THEN
BEBU=0.5
RPROD=BEBU/(1.+OM22)
L0FPRD=L0F(INDPR2)
ENDIF
IF(ISRMOD.EQ.3.OR.ISCHEM.EQ.1) THEN
CALL SUB2R(PREXS,PREXS2,FUEXS,FUEXS2)
CALL SUB2R(OXEXS,OXEXS2,EACTS,EACTS2)
ENDIF
ENDIF
IF(ISRMOD.EQ.2) L0MFU=L0F(INDFU)
IF(ISRMOD.EQ.3.OR.ISCHEM.EQ.1) THEN
TACTS=EACTS/1.987
PRECNS=1.E3*PREXS*WT(KWTFP)/1.E3**(FUEXS+OXEXS)
ITEMP=TEMP1
IF(TEMP0.GT.0.0) THEN
CALL FN2(EASP1,ITEMP,TEMP0,1.0)
ITEMP=EASP1
ENDIF
CALL SUB2(L0RHO,L0F(DEN1),L0TEM,L0F(ITEMP))
ENDIF
IF(ISRMOD.EQ.3) THEN
CALL FN1(CO,PRECNS)
ELSE
CALL FN15(CO,EP,KE,0.0,CREACT)
ENDIF
RFUMIN=IPRMOD.LE.2.AND.IPCHEM.EQ.1
cc IF(IPRPS.NE.0) L0PRPS=L0F(IPRPS)
C
DO 535 IX = IXF,IXL
IPLUS=(IX-1)*NY
CDIR$ IVDEP
DO 530 IY = IYF,IYL
I = IY + IPLUS
FU = F(L0VAR+I)
OX = F(L0MOX+I)
cc IF(IPRPS.NE.0) THEN
cc IF(F(L0PRPS+I).GE.SOLPRP) THEN
IF(SLD(I)) THEN
F(L0CO+I) = 0.0
F(L0SSOR+I) = 0.0
F(L0SVAL+I) = 0.0
GO TO 530
ENDIF
cc ENDIF
IF(ISRMOD.EQ.0) THEN
IF(COGRN.EQ.2) THEN
OXLIM = OX*RSTOIC
ELSEIF(COGRN.EQ.5) THEN
FU2 = F(L0MFU2+I)
OXLIM = FU*OX/((FU*STOIC+FU2*STOIC2)+TINY)
ENDIF
F(L0CO+I) = F(L0CO+I)*AMIN1(FU,OXLIM)/(FU+TINY)
ELSEIF(ISRMOD.EQ.1) THEN
PR=F(L0FPRD+I)
IF(COGRN.EQ.2) THEN
OXLIM = OX*RSTOIC
ELSEIF(COGRN.EQ.5) THEN
FU2 = F(L0MFU2+I)
OXLIM = FU*OX/((FU*STOIC+FU2*STOIC2)+TINY)
ENDIF
F(L0CO+I)=F(L0CO+I)*AMIN1(FU,OXLIM,RPROD*PR)
1 /(FU+TINY)
ELSEIF(ISRMOD.EQ.2) THEN
F(L0CO+I) = -F(L0P1RS+I)/(F(L0MFU+I)+TINY)
ELSEIF(ISRMOD.EQ.3) THEN
DENGAS=F(L0DEN+I)
FUMC=(DENGAS*FU/WT(KWTFP))**FUEXS
OXMC=(DENGAS*OX/WT(KOXID))**OXEXS
EXPTRM=EXP(-EACTS/F(L0TEM+I))
F(L0CO+I)=F(L0CO+I)*FUMC*OXMC*EXPTRM/
1 (DENGAS*(FU+TINY))
ENDIF
F(L0SSOR+I) = - F(L0CO+I)*FU
IF(RFUMIN) THEN
DENGAS=F(L0DEN+I)
FUMC=(DENGAS*FU/WT(KWTFP))**FUEXS
OXMC=(DENGAS*OX/WT(KOXID))**OXEXS
EXPTRM=EXP(-EACTP/F(L0TEM+I))
GSORAR=-PRECNS*FUMC*OXMC*EXPTRM/DENGAS
RFURAT=GSORAR/(GSOR+TINY)
IF(ABS(RFURAT).LT.1.0) THEN
F(L0SSOR+I)=GSORAR
F(L0CO+I)=-GSORAR/(FU+TINY)
ENDIF
ENDIF
IF(IFULIN.EQ.1) THEN
FUB = AMAX1(0.0, AMIN1(FU,FU-OXLIM) )
C FUMFUB=FU-FUB
C F(L0CO+I) = -F(L0SSOR+I)/(FUMFUB+TINY)
FUMFUB=AMAX1(TINY,FU-FUB)
F(L0CO+I) = -F(L0SSOR+I)/FUMFUB
F(L0SVAL+I) = FUB
ELSE
F(L0SVAL+I) = 0.0
ENDIF
530 CONTINUE
535 CONTINUE
ELSEIF(VALGRN.GE.0.AND.VALGRN.LE.5) THEN
CALL FN0(VAL,-L0SVAL)
ENDIF
ENDIF
C--------------------------------------------------------------------
CHAPTER 6 : ELEMENTAL COMPOSITION (Group 19 Section 5) ITASK=4
C--------------------------------------------------------------------
ELSEIF(ITASK.EQ.4) THEN
C--- Calculate the elemental composition and mole fractions for
C output purposes if those EL & M variables are stored for the
C whole field.
IF(STOCMP) THEN
LBCKF=NINT(CHSOB)
STOCMP=.FALSE.
DO 1951 K=1,KK
IF(STORE(LBCKF+K-1)) THEN
L0Y(K)=L0F(LBCKF+K-1)
STOCMP=.TRUE.
ELSE
L0Y(K)=0
ENDIF
1951 CONTINUE
ENDIF
C... Elemental mass composition
DO 1955 M=1,MM
IEL=LBNAME('EL'//SYMS(KK+M)(1:2))
IF(STORE(IEL)) THEN
L0EL=L0F(IEL)
DO 1954 J=1,NXNY
SUMEL=0.0
DO 1952 K=1,KK
YK=F(L0Y(K)+J)
SUMEL=SUMEL+NCF(M,K)*AWT(M)*YK/WT(K)
1952 CONTINUE
F(L0EL+J)=SUMEL
1954 CONTINUE
ENDIF
1955 CONTINUE
C... Species mole fractions
DO 1959 K=1,KK
IMOL=LBNAME('M'//SYMS(K)(1:3))
IF(STORE(IMOL)) THEN
KMOL=K
L0MOL=L0F(IMOL)
WTMIX=0.0
DO 1958 J=1,NXNY
WTMIX=0.0
DO 1956 KS=1,KK
YK=F(L0Y(KS)+J)
WTMIX=WTMIX+YK/WT(KS)
1956 CONTINUE
WTMIX=1./WTMIX
F(L0MOL+J)=F(L0Y(KMOL)+J)*WTMIX/WT(KMOL)
1958 CONTINUE
ENDIF
1959 CONTINUE
ENDIF
if(debug) call banner(2,'gxscrs',0)
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXCKTC(YMIX,TMP,TEMP0,ENTH,NTS,CP,ERROR)
C-----------------------------------------------------------------------
C GXCKTC: interfaces with the CHEMKIN Interface in order to
C calculate the TeMPerature from an input ENTHalpy and composition
C (YMIX-in mass-frac.s) when NTS>0, or to calculate the specific heat
C (CP) and ENTHalpy from the input TeMPerature and composition. If
C convergence fails in NTS iterations in the 1st case, then ERROR is
C set true.
C-----------------------------------------------------------------------
C--- The following section declares the REAL CHEMKIN interface
C buffer to be optionally double- or single-precision, and
C must be configured (using the CHANGE program) to match the
C precision chosen for CHEMKIN
C
C PARAMETER (NDCMR=2500,NDCMI=1000,NDCMC=100,NDCML=100)
C PARAMETER (NCKMI0=84)
C
C*****double precision
INCLUDE 'chmkin'
DOUBLE PRECISION
C*****END double precision
C*****single precision
C REAL
C*****END single precision
C 1 RCKMC,RU,RUC,PRES,CPMIX,RHOMIX,VISMIX,
C 2 CONMIX,VOLUME,ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,
C 3 UFAC,DFAC,U0,U,RTMP,QDOT,DTTP,DTMIN,DTMAX,DTEMP0
1 RU,RUC,CPMIX,
2 CONMIX,ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,
3 UFAC,DFAC,RTMP,DTMIN,DTMAX,DTEMP0
CHARACTER*16 CCKMC
COMMON /CKMCI0/LBCKF,KMCR,KMCI,K0MWT,K0BUF,LINK,LOUTCM,MM,KK,
1 II,L0HK,L0HKH,L0CP,L0CPH,L0H,L0HH,L0DK,L0DKH,
2 L0DJK,L0DJKH,L0VDE,L0VDN,L0VDH,L0CDK,L0TDK,L0HCO,
3 L0HVAL,L0DTF,K0L0,K0SYM,K0ELT,K0SYME,NTDIF,
4 K0ABV,K0BLO,K0TPB,K0TWP,K0IPV,K0SCR,K0VAR,
5 K0VARN,K0RES,K0RESN,K0A,K0SU,K0AP,L0SUK,L0APK,
6 K0H0K,LEVEL(2),LENTWP,ITIND,K0TD,
7 L0VHI,L0DTK,L0HSU,L0HAP,K0WDT,K0ROR,IGRND,
8 IQDOT,IENT,ICON,IMMWT,ICKOP1,IMIXF,IFUEL,IOXYG,
9 IPRO1,IPRO2,IDILU,L0APHI,L0SUHI,ICKOPT,IMCOPT,
1 ICALL1,LEVEL1,LEVEL2,NJAC,ITJAC,IRETIR,NUMDT,
2 NINIT,NRTOT,NITOT,NCTOT,NLTOT,K0EQRW,K0EQIW,K0AWT
C 1 /CKMCI/ICKMC(NDCMI)
2 /CKMCR0/ATIM,RTIM,ATOL,RTOL,UFAC,DFAC,ABSOL,RELAT,
2 RU,RUC,DTMIN,DTMAX,DTEMP0
C 2 /CKMCR/RCKMC(NDCMR)
C 3 /CKMCC/CCKMC(NDCMC)
4 /CKMCL0/VARPRS,CONTEM
C 4 /CKMCL/LCKMC(NDCML)
C
LOGICAL ERROR
REAL YMIX(KK)
C
ENTH0=0.0
DO 10 K=1,KK
10 RCKMC(K0BUF+K)=YMIX(K)
C
ERROR=.FALSE.
C
RTMP=TMP
NTS=MAX(-1,NTS)
IF(NTS.LT.0) ENTH=0.0
DO 100 ITS=1,ABS(NTS)
C... Determine mean mixture enthalpy (CONMIX)in mass units
CALL CKHBMS((RTMP+TEMP0),RCKMC(K0BUF+1),ICKMC,RCKMC,CONMIX)
C... Determine mean mixture Cp in ergs/(gm*K)
CALL CKCPBS((RTMP+TEMP0),RCKMC(K0BUF+1),ICKMC,RCKMC,CPMIX)
CONMIX=CONMIX-ENTH0
DTMP=-(CONMIX-ENTH)/CPMIX
RTMP=RTMP+DTMP
IF(ABS(DTMP/RTMP).LT.1.E-5.AND.
1 ABS(1.0-ENTH/CONMIX).LT.1.E-5.OR.NTS.LT.0) GO TO 200
100 CONTINUE
C
ERROR=NTS.GT.0
RETURN
C
200 IF(NTS.GT.0) THEN
TMP=RTMP
ELSE
ENTH=CONMIX
ENDIF
CP=CPMIX
C
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXCKIN(KK0,KKK,III,MMM,WT,AWT,RGAS,SYMS,NU,NCF,L0Y,
1 EQUSOL)
C-----------------------------------------------------------------------
C GXCKIN: interfaces with the CHEMKIN Interface in order to
C return various bits of data to the SCRS code (GXSCRS) for
C initialisation purposes.
C-----------------------------------------------------------------------
C--- The following section declares the REAL CHEMKIN interface
C buffer to be optionally double- or single-precision, and
C must be configured (using the CHANGE program) to match the
C precision chosen for CHEMKIN
C
INCLUDE 'chmkin'
PARAMETER (NDCMR=2500,NDCMI=1000,NDCMC=100,NDCML=100)
C PARAMETER (NCKMI0=84)
C
C*****double precision
DOUBLE PRECISION
C*****END double precision
C*****single precision
C REAL
C*****END single precision
C 1 RCKMC,RU,RUC,PRES,CPMIX,RHOMIX,VISMIX,
C 2 CONMIX,VOLUME,ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,
C 3 UFAC,DFAC,U0,U,RTMP,QDOT,DTTP,DTMIN,DTMAX,DTEMP0
1 RU,RUC,
2 ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,
3 UFAC,DFAC,DTMIN,DTMAX,DTEMP0
CHARACTER*16 CCKMC
COMMON /CKMCI0/LBCKF,KMCR,KMCI,K0MWT,K0BUF,LINK,LOUTCM,MM,KK,
1 II,L0HK,L0HKH,L0CP,L0CPH,L0H,L0HH,L0DK,L0DKH,
2 L0DJK,L0DJKH,L0VDE,L0VDN,L0VDH,L0CDK,L0TDK,L0HCO,
3 L0HVAL,L0DTF,K0L0,K0SYM,K0ELT,K0SYME,NTDIF,
4 K0ABV,K0BLO,K0TPB,K0TWP,K0IPV,K0SCR,K0VAR,
5 K0VARN,K0RES,K0RESN,K0A,K0SU,K0AP,L0SUK,L0APK,
6 K0H0K,LEVEL(2),LENTWP,ITIND,K0TD,
7 L0VHI,L0DTK,L0HSU,L0HAP,K0WDT,K0ROR,IGRND,
8 IQDOT,IENT,ICON,IMMWT,ICKOP1,IMIXF,IFUEL,IOXYG,
9 IPRO1,IPRO2,IDILU,L0APHI,L0SUHI,ICKOPT,IMCOPT,
1 ICALL1,LEVEL1,LEVEL2,NJAC,ITJAC,IRETIR,NUMDT,
2 NINIT,NRTOT,NITOT,NCTOT,NLTOT,K0EQRW,K0EQIW,K0AWT
C 1 /CKMCI/ICKMC(NDCMI)
2 /CKMCR0/ATIM,RTIM,ATOL,RTOL,UFAC,DFAC,ABSOL,RELAT,
2 RU,RUC,DTMIN,DTMAX,DTEMP0
C 2 /CKMCR/RCKMC(NDCMR)
C 3 /CKMCC/CCKMC(NDCMC)
4 /CKMCL0/VARPRS,CONTEM
C 4 /CKMCL/LCKMC(NDCML)
1 /EQWR/REQWR
C
REAL WT(*),AWT(*)
CHARACTER*(*) SYMS(*)
INTEGER NU(KK,*),L0Y(*),NCF(MM,KK)
LOGICAL EQUSOL
C
IF(KK0.GT.0) THEN
IF(KK.NE.KK0) THEN
CALL WRIT40('GXCKIN: NO. OF SPECIES IS INCORRECT ')
CALL SET_ERR(289,
* 'Error. GXCKIN: NO. OF SPECIES IS INCORRECT.',1)
C CALL EAROUT(1)
ENDIF
C
LSYM=LEN(SYMS(1))
IF(LSYM.LT.16)
1 CALL WRIT40('GXCKIN: SYMBOL-NAMES WILL BE TRUNCATED ')
RGAS=RU
CALL SUB3(MMM,MM,KKK,KK,III,II)
C
C... Copy the species-names and molecular weights from the CHEMKIN
C data and get by means of CKNU the reaction data, ie. number of
C moles created or destroyed of each species for each reaction.
C
DO 10 K=1,KK
SYMS(K)=CCKMC2(K0SYM+K)
10 WT(K)=RCKMC(K0MWT+K)
CALL CKNU(KK,ICKMC,RCKMC,NU)
ENDIF
DO 12 M=1,MM
SYMS(KK+M)=CCKMC2(K0SYME+M)
AWT(M)=RCKMC(K0AWT+M)
KEM=K0ELT+M
DO 11 K=1,KK
NCF(M,K)=ICKMC(KEM)
KEM=KEM+MM
11 CONTINUE
12 CONTINUE
C
IF(EQUSOL) THEN
NCON = 0
CALL EQLEN (MM, KK, NCON, LENIEQ, LENREQ)
IF(NRTOT+LENREQ.GT.NDCMR) THEN
CALL WRIT40(' NOT ENOUGH REAL SPACE')
CALL SET_ERR(290, 'Error. NOT ENOUGH REAL SPACE.',1)
C CALL EAROUT(1)
ENDIF
IF(NITOT+LENIEQ.GT.NDCMI) THEN
CALL WRIT40(' NOT ENOUGH INTEGER SPACE')
CALL SET_ERR(291, 'Error. NOT ENOUGH INTEGER SPACE.',1)
C CALL EAROUT(1)
ENDIF
K0EQRW=NRTOT
K0EQIW=NITOT
NRTOT=NRTOT+LENREQ
NITOT=NITOT+LENIEQ
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE INCOND(ITASK,KK,SYMS,YFU,YOX,TMPFU,TMPOX,ENTHFU,
1 ENTHOX,RHOFU,RHOOX)
C-----------------------------------------------------------------------
C INCOND: this routine simply prints out inlet conditions for the
C fuel and oxidant streams. It is fairly self explanatory.
C-----------------------------------------------------------------------
REAL YFU(KK),YOX(KK)
CHARACTER*(*) SYMS(KK)
CHARACTER*40 LINE
C
LSYM=LEN(SYMS(1))
C
CALL WRITBL
CALL WRITST
CALL WRITBL
C
IF(ITASK.LT.0) THEN
CALL WRIT40(' INLET COMPOSITION IN MOLE-FRACTIONS IS:')
ELSE
CALL WRIT40(' INLET COMPOSITION IN MASS-FRACTIONS IS:')
ENDIF
CALL WRIT40(' SPECIES FUEL OXIDANT ')
CALL WRIT40('MIXF 1 0 ')
DO 100 K=1,KK
LINE=SYMS(K)
WRITE(LINE(17:40),'(1P2E10.3)') YFU(K),YOX(K)
CALL WRIT40(LINE)
100 CONTINUE
C
IF(ITASK.GT.0) THEN
LINE='TEMPERATURE'
WRITE(LINE(17:40),'(1P2E10.3)') TMPFU,TMPOX
CALL WRIT40(LINE)
c IF(ITASK.GT.1) THEN
LINE='ENTHALPY'
WRITE(LINE(17:40),'(1P2E10.3)') ENTHFU,ENTHOX
CALL WRIT40(LINE)
LINE='DENSITY'
WRITE(LINE(17:40),'(1P2E10.3)') RHOFU,RHOOX
CALL WRIT40(LINE)
c ENDIF
ENDIF
C
CALL WRITBL
CALL WRITST
CALL WRITBL
C
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXPHTC(YMIX,TMP,TEMP0,ECO,NP,ENTH,NTS,CP,ERROR)
C-----------------------------------------------------------------------
C GXPHTC: interfaces with the PHOENICS Chemistry system in order to
C calculate the TeMPerature from an input ENTHalpy and composition
C (YMIX-in mass-frac.s) when NTS>0, or to calculate the specific heat
C (CP) and ENTHalpy from the input TeMPerature and composition. If
C convergence fails in NTS iterations in the 1st case, then ERROR is
C set true. The same as GXCKTC, but note the extra arguments which
C carry the data for the thermodynamic functions, and the different
C way in which CKHMS & GXHMSK and CKCPMS & GXCMSK work. Note also
C the difference in use of double-precision and te slightly more
C complex iteration method (which MAY well not be needed) in this
C routine.
C-----------------------------------------------------------------------
C
REAL YMIX(*),ECO(*)
INTEGER NP(*)
LOGICAL ERROR
C
ERROR=.FALSE.
C
RTMP=0.0
NTS=MAX(-1,NTS)
IF(NTS.LT.0) ENTH=0.0
DO 100 ITS=1,ABS(NTS)
CALL GXHBMS((TMP+RTMP+TEMP0),YMIX,NP,ECO,ENTHA)
CALL GXCPBS((TMP+RTMP+TEMP0),YMIX,NP,ECO,CP)
DTMP=-(ENTHA-ENTH)/CP
RTMP=RTMP+DTMP
IF(ABS(DTMP/(TMP+RTMP)).LT.1.E-5.AND.
1 ABS(1.0-ENTH/ENTHA).LT.1.E-5.OR.NTS.LT.0) GO TO 200
100 CONTINUE
C
ERROR=NTS.GT.0
RETURN
C
200 IF(NTS.LE.0) THEN
ENTH=ENTHA
ELSE
TMP=TMP+RTMP
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXCOMP(IREACT,KK,YFU,YOX,L0Y,FMIX,STEP1,STEP2,TWOFPS,
1 SOLVP2,KFUEL,KOXID,KFP1,KFP2,KPROD1,KPROD2,
2 OM1,OM12,OM22,FPM1,FPM2,PM1,PM2,I,YMIX)
INCLUDE 'farray'
REAL YMIX(KK),YFU(KK),YOX(KK)
INTEGER L0Y(KK)
LOGICAL STEP1,STEP2,TWOFPS,SOLVP2
C... Unburnt composition
FMIX1=1.-FMIX
DO 1 K=1,KK
YMIX(K)=FMIX*YFU(K)+FMIX1*YOX(K)
1 CONTINUE
C
C... Fast-chemistry composition calculation
IF(IREACT.EQ.1) THEN
IF(STEP1) THEN
C... If the first step is active, modify the composition.
DOXID=MIN(OM1*YMIX(KFUEL),YMIX(KOXID))
IF(DOXID.GT.0.0) THEN
YMIX(KOXID)=YMIX(KOXID)-DOXID
DFUEL=DOXID/OM1
YMIX(KFUEL)=YMIX(KFUEL)-DFUEL
YMIX(KFP1) =YMIX(KFP1) +DFUEL*FPM1
IF(TWOFPS) YMIX(KFP2) =YMIX(KFP2)+DFUEL*FPM2
ENDIF
ENDIF
IF(STEP2.AND.YMIX(KOXID).GT.0.0) THEN
C... If the second step is active, modify the composition.
OFPSUM=OM12*YMIX(KFP1)
IF(TWOFPS) OFPSUM=OFPSUM+OM22*YMIX(KFP2)
DOXID=MIN(OFPSUM,YMIX(KOXID))
IF(DOXID.GT.0.0) THEN
YMIX(KOXID) =YMIX(KOXID) -DOXID
DFP1=DOXID*(YMIX(KFP1)/OFPSUM)
YMIX(KFP1) =YMIX(KFP1) -DFP1
YMIX(KPROD1)=YMIX(KPROD1)+DFP1*PM1
IF(TWOFPS) THEN
DFP2=DOXID*(YMIX(KFP2)/OFPSUM)
YMIX(KFP2) =YMIX(KFP2) -DFP2
YMIX(KPROD2)=YMIX(KPROD2)+DFP2*PM2
ENDIF
ENDIF
ENDIF
C
C... Finite-Rate Chemistry Composition Calculation
ELSEIF(IREACT.EQ.2) THEN
IF(STEP1) THEN
YKFUEL=F(L0Y(KFUEL)+I)
DFUEL=YKFUEL-YMIX(KFUEL)
YMIX(KFUEL)=YKFUEL
DOXID=MAX(OM1*DFUEL,-YMIX(KOXID))
DFUEL=DOXID/OM1
IF(DOXID.LT.0.0) THEN
YMIX(KOXID)= YMIX(KOXID) + DOXID
YMIX(KFP1) = YMIX(KFP1) - DFUEL*FPM1
IF(TWOFPS) YMIX(KFP2) = YMIX(KFP2) - DFUEL*FPM2
ENDIF
ENDIF
c... is the following 'and.ymix(koxid).gt.0.0' now valid
IF(STEP2.AND.YMIX(KOXID).GT.0.0) THEN
C... If the second step is active, modify the composition
YKFP1 =F(L0Y(KFP1)+I)
DFP1=YKFP1-YMIX(KFP1)
DOXID=MAX(OM12*DFP1,-YMIX(KOXID))
DFP1=DOXID/OM12
IF(DOXID.LT.0.0) THEN
YMIX(KOXID) = YMIX(KOXID) + DOXID
YMIX(KFP1) = YMIX(KFP1) + DFP1
YMIX(KPROD1) = YMIX(KPROD1) - DFP1*PM1
IF(TWOFPS.AND.SOLVP2) THEN
YKFP2 =F(L0Y(KFP2)+I)
DFP2 = YKFP2 - YMIX(KFP2)
DOXID = MAX(OM22*DFP2,-YMIX(KOXID))
DFP2 = DOXID/OM22
YMIX(KFP2) = YMIX(KFP2) + DFP2
YMIX(KOXID) = YMIX(KOXID) + DOXID
C... mrm 31/05/96 prevent negative products
YMIX(KPROD2)= AMAX1(YMIX(KPROD2),0.)
YMIX(KPROD2)= YMIX(KPROD2) - DFP2*PM2
ENDIF
ENDIF
ENDIF
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXDPDF(G,F,FP,FM,FA)
FPRIME=SQRT(G)
FPP=F+FPRIME
FPM=F-FPRIME
FP =MIN(FPP,1.0)
FM =MAX(FPM,0.0)
IF(FPP.GE.1.0.AND.FM.GT.0.0) THEN
FM=F-G/((1.-F)+1.E-20)
ELSEIF(FPM.LE.0.0.AND.FP.LT.1.0) THEN
FP=F+G/(F+1.E-20)
ENDIF
FPMFM=AMAX1(1.E-20,FP-FM)
FA=(F-FM)/FPMFM
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
SUBROUTINE GXCKEQ(LU,ITASK,YMIX,TMP,TEMP0,PRESS0,WTMIX,ERROR)
C-----------------------------------------------------------------------
C GXCKEQ: interfaces with the CHEMKIN Interface
C-----------------------------------------------------------------------
C--- The following section declares the REAL CHEMKIN interface
C buffer to be optionally double- or single-precision, and
C must be configured (using the CHANGE program) to match the
C precision chosen for CHEMKIN
C
C PARAMETER (NDCMR=2500,NDCMI=1000,NDCMC=100,NDCML=100)
C PARAMETER (NCKMI0=84)
C
INCLUDE 'chmkin'
C*****double precision
DOUBLE PRECISION
C*****END double precision
C*****single precision
C REAL
C*****END single precision
C 1 RCKMC,RU,RUC,PRES,CPMIX,RHOMIX,VISMIX,
C 2 CONMIX,VOLUME,ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,
C 3 UFAC,DFAC,U0,U,RTMP,QDOT,DTTP,DTMIN,DTMAX,DTEMP0
1 RU,RUC,PRES,CPMIX,
2 CONMIX,ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,
3 UFAC,DFAC,RTMP,DTMIN,DTMAX,DTEMP0
CHARACTER*16 CCKMC
COMMON /CKMCI0/LBCKF,KMCR,KMCI,K0MWT,K0BUF,LINK,LOUTCM,MM,KK,
1 II,L0HK,L0HKH,L0CP,L0CPH,L0H,L0HH,L0DK,L0DKH,
2 L0DJK,L0DJKH,L0VDE,L0VDN,L0VDH,L0CDK,L0TDK,L0HCO,
3 L0HVAL,L0DTF,K0L0,K0SYM,K0ELT,K0SYME,NTDIF,
4 K0ABV,K0BLO,K0TPB,K0TWP,K0IPV,K0SCR,K0VAR,
5 K0VARN,K0RES,K0RESN,K0A,K0SU,K0AP,L0SUK,L0APK,
6 K0H0K,LEVEL(2),LENTWP,ITIND,K0TD,
7 L0VHI,L0DTK,L0HSU,L0HAP,K0WDT,K0ROR,IGRND,
8 IQDOT,IENT,ICON,IMMWT,ICKOP1,IMIXF,IFUEL,IOXYG,
9 IPRO1,IPRO2,IDILU,L0APHI,L0SUHI,ICKOPT,IMCOPT,
1 ICALL1,LEVEL1,LEVEL2,NJAC,ITJAC,IRETIR,NUMDT,
2 NINIT,NRTOT,NITOT,NCTOT,NLTOT,K0EQRW,K0EQIW
C 1 /CKMCI/ICKMC(NDCMI)
2 /CKMCR0/ATIM,RTIM,ATOL,RTOL,UFAC,DFAC,ABSOL,RELAT,
2 RU,RUC,DTMIN,DTMAX,DTEMP0
C 2 /CKMCR/RCKMC(NDCMR)
C 3 /CKMCC/CCKMC(NDCMC)
4 /CKMCL0/VARPRS,CONTEM
C 4 /CKMCL/LCKMC(NDCML)
C
LOGICAL ERROR,LEQUI,LCNTN
REAL YMIX(KK)
C
DO 100 K=1,KK
RCKMC(K0BUF+K)=YMIX(K)
100 CONTINUE
C
IF(ITASK.EQ.0) THEN
LEQUI=.TRUE.
LCNTN=.FALSE.
ELSE
LEQUI=.FALSE.
LCNTN=.TRUE.
ENDIF
IOPT=5
KMON=1
RTMP=TMP+TEMP0
PRES=PRESS0
NCONST=0
C
c all arguments must be declared in the real/double precision
c section at the top of this routine
C
CALL EQUIL(LU,.FALSE.,0,LEQINI,LCNTN,ICKMC,RCKMC,
1 LENIEQ,ICKMC(K0EQIW+1),LENREQ,RCKMC(K0EQRW+1),
1 MM,KK,CCKMC1(K0ELT+1),CCKMC1(K0SYM+1),IOPT,KMON,
2 RCKMC(K0BUF+1),RTMP,RTMP,PRES,PRES,NCONST,
3 RCKMC(K0VAR+1),RCKMC(K0VAR+1))
C
C... Next gets solution using CPMIX as a dummy variable and CONMIX
C to pick up the mixture molecular weight
C
CALL EQSOL (KK, RCKMC(K0EQRW+1), RCKMC(K0BUF+1), RCKMC(K0BUF+1),
1 TMP, CPMIX, CPMIX, CPMIX, CPMIX, CONMIX, CPMIX,
2 CPMIX)
C
DO 200 K=1,KK
YMIX(K)=RCKMC(K0BUF+K)
200 CONTINUE
WTMIX=CONMIX
TMP=RTMP-TEMP0
END
c