cgxbfgr.for

c

File-name GXSCRS.htm 20803

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 '/phoenics/d_includ/chmkin' INCLUDE '/phoenics/d_includ/patcmn' INCLUDE '/phoenics/d_includ/farray' INCLUDE '/phoenics/d_includ/satear' INCLUDE '/phoenics/d_includ/grdloc' INCLUDE '/phoenics/d_includ/satgrd' INCLUDE '/phoenics/d_includ/grdear' INCLUDE '/phoenics/d_includ/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 '/phoenics/d_includ/chmkin' INCLUDE '/phoenics/d_includ/farray' INCLUDE '/phoenics/d_includ/satear' INCLUDE '/phoenics/d_includ/grdloc' INCLUDE '/phoenics/d_includ/satgrd' INCLUDE '/phoenics/d_includ/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)) ) C FUMFUB=FU-FUB C F(L0CO+I) = -F(L0P1RS+I)/(FUMFUB+TINY) 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 '/phoenics/d_includ/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 '/phoenics/d_includ/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 '/phoenics/d_includ/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 C FA=(F-FM)/((FP-FM)+1.E-20) 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 '/phoenics/d_includ/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