cgxbfgr.for

c
c
c
c
c
      SUBROUTINE GXNOX
C---------------------------------------------------------------
      INCLUDE 'farray'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE 'grdbfc'
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
C================================================================
C 1         CREK COMMON BLOCKS
C          ====================
      LOGICAL SLD
      COMMON/CINDEX/JHCPS,JCK1,JITER,JCK2(15)
      COMMON/CKINET/GCPSUM,GCK1(8)
      LOGICAL LADIAB,LCONVG,LDEBUG,LEQUIL,LREACT,LNRG
      COMMON/CHMLGC/LADIAB,LCONVG,LDEBUG,LEQUIL,LREACT,LNRG
      COMMON/CSPECE/GH0(20),GCK2(20),GSMW(20),GCK3(300)
      COMMON/CPARAM/GEMV,GHSUB0,JNS,GPA,GQ0,GQ1,GQ2,
     1              GQ3,GQ4,GRHOP,GSM,GS1(20),GS2(20),GTK
C================================================================
      SAVE GSNAME
      SAVE TEN1, TEN2, GBX1, GBX2, TACT1, TACT2, GFSN,
     1     L0SPR2, L0SPR3, L0SPR4, RSPR, RSOUR, ID
      SAVE JXDBCK, JYDBCK, JZDBCK, JSDBCK, NS, NR
      SAVE CRKCAL, CRKBUG, DBGNOX
      SAVE GN2PR, GN2OX, GCO2PR, GH2OPR, GHM, GOM, GNM,
     1     GNOM, GN2M, GOHM, GO2M, GCO2M, GH2OM, GASCON,
     1     GCOM, GH2M, GFUM, SETMIN
      SAVE JXNO, JXN, JXO, JXH, JXOH, JNOSR, JYOXID,
     1     JXO2, JYO2, JYCO2, JYH2, JYH2O, JYCO, JEQUI,
     1     JDEGF, JYN2, JYPR, LUREAD, IERR, JTMP
      SAVE JCRKT, L0VOL, GFLOT, GFNOT, GFNOVT, GFLNT,
     1     GEXC, GEXH, GEXN, GEXO
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION STARTS:
      PARAMETER (NSM=10,NRM=3)
      CHARACTER*4 GSNAME
      DIMENSION GSNAME(8)
      DIMENSION RSOUR(NSM),RSPR(NSM)
      DIMENSION GBX1(NRM),GBX2(NRM),TACT1(NRM),TACT2(NRM)
      DIMENSION TEN1(NRM),TEN2(NRM),ID(4,NRM),GFSN(NSM)
C
      LOGICAL CRKCAL,QEQ,CRKBUG,DBGNOX
C
      DATA GSNAME/'CO ','CO2','H  ','H2 ','H2O','O  ','OH ','O2'/
C Zeldovich rate constants
C   GBX1 & GBX2 are in units of kg**2/m**3/kg-mol/s when TEN1
C             and TEN2 are zero.
      DATA GBX1/7.6E10,6.4E6,6.3E8/
      DATA GBX2/1.6E10,1.5E6,2.5E9/
C   TACT1 & TACT2 are in units of K ( i.e. E/Ro ) where
C     is the activation energy & Ro is the gas constant.
      DATA TACT1/3.8E4,3150.0,0.0/
      DATA TACT2/0.0,1.95E4,2.446E4/
C   TEN1 & TEN2 are exponents for the temperature
      DATA TEN1/0.0,1.0,0.5/
      DATA TEN2/0.0,1.0,0.5/
      DATA ID/4,5,3,2,
     1        7,2,3,5,
     1        6,2,3,1/
C
      NAMSUB='GXNOX'
      IXL=IABS(IXL)
C     Molecular mass in Kg/(Kg mol)
C*****************************************************************
C--- GROUP 1. Run title and other preliminaries
C*****************************************************************
      IF(IGR.EQ.1.AND.ISC.EQ.1) THEN
C     User may here change message transmitted to the VDU screen or
C     batch-run log file.
        CALL WRYT40('ground file is gxnox.f of 25.04.99      ')
        CALL WRYT40('PHOENICS version number is 3.2          ')
C
C.... Read information from Q1
C     integers
        CALL RQ1I('NOXD','NS',NS)
        CALL RQ1I('NOXD','NR',NR)
        IF(LSG63) THEN
          CALL RQ1I('NOXD','JXDBCK',JXDBCK)
          CALL RQ1I('NOXD','JYDBCK',JYDBCK)
          CALL RQ1I('NOXD','JZDBCK',JZDBCK)
          CALL RQ1I('NOXD','JSDBCK',JSDBCK)
C           reals
          CALL RQ1R('NOXD','GN2OX',GN2OX)
          CALL RQ1R('NOXD','GN2PR',GN2PR)
          CALL RQ1R('NOXD','GCO2PR',GCO2PR)
          CALL RQ1R('NOXD','GH2OPR',GH2OPR)
C           logicals
          CALL RQ1L('NOXD','CRKBUG',CRKBUG)
          CALL RQ1L('NOXD','DBGNOX',DBGNOX)
        ENDIF
C
        CALL SUB3R(GHM,1.00797,GOM,15.9994,GNM,14.0067)
        CALL SUB2R(GNOM,30.0061,GN2M,28.0134)
        CALL SUB3R(GOHM,17.00737,GO2M,31.9988,GCO2M,44.00995)
        CALL SUB3R(GH2OM,18.01534,GASCON,8314.0,SETMIN,1.0E-10)
C
        IF(LSG62) THEN
          CALL SUB2R(GCOM,28.01055,GH2M,2.01594)
        ELSEIF(LSG63) THEN
          GFUM=16.04303
        ENDIF
C
        CALL SUB2(JXNO,LBNAME('XNO'),JXN,LBNAME('XN'))
        CALL SUB2(JXO,LBNAME('XO'),JXH,LBNAME('XH'))
        CALL SUB2(JXOH,LBNAME('XOH'),JNOSR,LBNAME('NOSR'))
        CALL SUB2(JCRKT,LBNAME('CRKT'),JXO2,LBNAME('XO2'))
        IF(LSG62) THEN
          CALL SUB2(JYCO2,LBNAME('YCO2'),JYH2,LBNAME('YH2'))
          CALL SUB2(JYH2O,LBNAME('YH2O'),JYCO,LBNAME('YCO '))
          CALL SUB2(JEQUI,LBNAME('EQUI'),JDEGF,LBNAME('DEGF'))
          CALL SUB2(JYN2,LBNAME('YN2'),JYO2,LBNAME('YO2'))
          JTMP=LBNAME('TMP1')
        ELSEIF(LSG63) THEN
          JTMP=LBNAME('TMP1')
          JYPR=LBNAME('PROD')
          JYOXID=LBNAME('OXID')
        ENDIF
C
        CRKCAL=.TRUE.
        LUREAD=20
        IERR=-1
        CALL OPENZZ(29,IERR)
        IF(IERR.NE.0) THEN
          CALL WRIT40('ERROR OPENING CH4DAT IN GXNOX         ')
          CALL WRIT1I('ERROR NO',IERR)
          CALL SET_ERR(252, 'ERROR OPENING CH4DAT IN GXNOX.',1)
C          CALL EAROUT(1)
        ENDIF
        CALL GXCRK0( LUREAD , LUPR1 )
        GEMV=0.0
        LEQUIL=.TRUE.
        CALL GXMAKE(L0SPR2,NX*NY,'SPR2')
        CALL GXMAKE(L0SPR3,NX*NY,'SPR3')
        CALL GXMAKE(L0SPR4,NX*NY,'SPR4')
C*****************************************************************
C--- GROUP 13. Boundary conditions and special sources
C*****************************************************************
      ELSEIF(IGR.EQ.13.AND.ISC.EQ.1) THEN
C--coefficient = GRND
        IF(NPATCH.EQ.'NOXSR') CALL FN0(CO,-L0SPR2)
      ELSEIF(IGR.EQ.13.AND.ISC.EQ.12) THEN
C--value = GRND
        IF(NPATCH.EQ.'NOXSR') THEN
          IF(INDVAR.EQ.JXNO) THEN
            CALL FN0(VAL,-L0SPR3)
C... evaluate and store local NO source per unit volume of gas
            CALL SUB4(L0CO,L0SPR2,L0VAL,L0SPR3,L0SOR,L0F(JNOSR),
     1                L0NO,L0F(JXNO))
            DO 1311 J=1,NX*NY
              F(L0SOR+J)=F(L0CO+J)*(F(L0VAL+J)-F(L0NO+J))
 1311       CONTINUE
          ELSEIF(INDVAR.EQ.JXN) THEN
            CALL FN0(VAL,-L0SPR4)
          ENDIF
        ENDIF
***************************************************************
C--- GROUP 19. SECTION 4 ---- START OF ITERATION.
C***************************************************************
      ELSEIF(IGR.EQ.19.AND.ISC.EQ.4) THEN
        IF(CRKCAL) THEN
          CALL SUB3(L0XH,L0F(JXH),L0XO,L0F(JXO),L0XOH,L0F(JXOH))
          CALL SUB3(L0XO2,L0F(JXO2),L0CRKT,L0F(JCRKT),L0VOL,L0F(VOL))
          IF(LSG62) THEN
            CALL SUB3(L0YCO,L0F(JYCO),L0YO2,L0F(JYO2),L0YN2,L0F(JYN2))
            CALL SUB4(L0YCO2,L0F(JYCO2),L0YH2O,L0F(JYH2O),L0YH2,
     1                 L0F(JYH2),L0FTMP,L0F(JTMP))
          ELSEIF(LSG63) THEN
            CALL SUB2(L0FYPR,L0F(JYPR),L0FTMP,L0F(JTMP))
            L0OXID=L0F(JYOXID)
          ENDIF
          IF(L0FTMP.EQ.0) THEN
            WRITE(14,*) 'Temperature is not stored. Execution ended.'
            CALL SET_ERR(253, 'Error. Temperature is not stored.',1)
C            CALL EAROUT(1)
          ENDIF
C**Equilibrium Calculation
          if(dbgrnd) then
            call writst
            call writbl
            call writ40('Crek calculation starts for current slab')
            call writ1i('IZ      ',IZ)
            call writbl
            call writst
          endif
          DO 1946 JX=1,NX
            DO 1946 JY=1,NY
              J = JY + NY*(JX-1)
C--Skip crek call for blocked cells
              IF(SLD(J)) GO TO 1946
          IF(L0VOL.NE.0) THEN
            IF(F(L0VOL+J).LE.TINY) GO TO 1946
          ENDIF
C
C--Supply all species except N2 for equilibrium calculation
C   H2O = H +OH     ; CO2 = CO +0.5O2 ; O2 = 2O
C   H2O = 0.5H2 + OH ; H2 = 2H
C
              IF(LSG62) THEN
            GYPR=1.0-F(L0YN2+J)
            GYN2=F(L0YN2+J)
C
              ELSEIF(LSG63) THEN
            GCPR = F(L0FYPR +J)
            GYPR = AMAX1(TINY,GCPR - GN2PR*GCPR)
            GYN2 = GN2PR*GCPR+GN2OX*(F(L0OXID+J))
              ENDIF
C
C--Skip crek call if only N2 present
              IF(GYN2.GT.0.98) GO TO 1946
C
C--Initialise CREK arrays S1 & S2 with mol fractions
C
              DO 1941 IS=1,JNS
C CO
        IF(IS.EQ.1.AND.LSG62) THEN
                  GS1(IS)=F(L0YCO+J)/(GCOM*GYPR)
                  GS2(IS)=GS1(IS)
C CO2(coal)
                ELSEIF(IS.EQ.2) THEN
                  IF(LSG62) THEN
                    GS1(IS)=F(L0YCO2+J)/(GCO2M*GYPR)
C CO2(gas)
                  ELSEIF(LSG63) THEN
                    GS1(IS)=GCPR*GCO2PR/(GCO2M*GYPR)
                  ENDIF
                  GS2(IS)=GS1(IS)
C H2
                ELSEIF(IS.EQ.4.AND.LSG62) THEN
                  GS1(IS)=F(L0YH2+J)/(GH2M*GYPR)
                  GS2(IS)=GS1(IS)
C H2O(coal)
                ELSEIF(IS.EQ.5) THEN
                  IF(LSG62) THEN
                    GS1(IS)=F(L0YH2O+J)/(GH2OM*GYPR)
C H2O(gas)
                  ELSEIF(LSG63) THEN
                    GS1(IS)=GCPR*GH2OPR/(GH2OM*GYPR)
                  ENDIF
                  GS2(IS)=GS1(IS)
C O2
                ELSEIF(IS.EQ.8.AND.LSG62) THEN
                  GS1(IS)=F(L0YO2+J)/(GO2M*GYPR)
                  GS2(IS)=GS1(IS)
C H,O & OH species
                ELSE
                  GS1(IS)=SETMIN
                  GS2(IS)=GS1(IS)
                ENDIF
1941          CONTINUE
C
C Temperature, pressure & other CREK input parameters
              GTK   = F(L0FTMP+J)
              GPA   = PRESS0
              JHCPS = 1
C Enter CREK routine to determine enthalpy
              CALL HCPS
C
              GHSUM=0.0
              DO 1942 IS=1,JNS
            GHSUM=GHSUM+GH0(IS)*GS2(IS)
1942      CONTINUE
C
              GHSUB0=GHSUM*GASCON*GTK
C
C  Enter CREK to perform equilibrium dissociation calculation
              CALL GXCREK
C  H mass fraction from CREK
              F(L0XH+J)=GS2(3)*GHM*GYPR
C  O mass fraction from CREK
              F(L0XO+J)=GS2(6)*GOM*GYPR
C  OH mass fraction from CREK
              F(L0XOH+J)=GS2(7)*GOHM*GYPR
C  O2 mass fraction from CREK
              F(L0XO2+J)=GS2(8)*GO2M*GYPR
C  T from CREK equilibrium calculation
              F(L0CRKT+J)=GTK
C--Check total mass of array GS2 = 1.0
              GFRAC=0.0
              DO 1944 IS=1,JNS
                GFRAC=GFRAC+GS2(IS)*GSMW(IS)
1944          CONTINUE
C
              IF(GFRAC.LT.0.995) THEN
                WRITE(LUPR1,*)' gs1 input to crek;',
     1                          ' gs2 output from crek'
                WRITE(LUPR1,*)' mol fractions of',
     1                           'equilibrium species'
                WRITE(LUPR1,*)'           GS1',
     1                           '        GS2       GH0 '
                DO 1945 IS=1,JNS
                  WRITE(LUPR1,951) IS,GSNAME(IS),
     1               GS1(IS),GS2(IS),GH0(IS)
1945            CONTINUE
                CALL WRITBL
                WRITE(LUPR1,*)'sum of mass fractions : ',GFRAC
                WRITE(LUPR1,*)'Temperature= ',GTK,' Hsub0 = ',GHSUB0
                CALL WRITBL
                IF(LSG62) THEN
                  WRITE(LUPR1,*) ' mass fractions',' including N2'
                  WRITE(LUPR1,*)'            input','    output    '
                  DO 1947 IS = 1 , JNS
                    WRITE(LUPR1,951) IS,GSNAME(IS),
     1                                  GS1(IS)*GSMW(IS)*GYPR,
     1                                  GS2(IS)*GSMW(IS)*GYPR
1947              CONTINUE
                  WRITE(LUPR1,953) F(L0YN2+J),F(L0YN2+J)
                ENDIF
 951            FORMAT(1X,I2,1X,A4,1P4E11.3)
 953            FORMAT(1X,' 9 N2  ',1P2E11.3)
              ENDIF
1946      CONTINUE
          CALL WRITST
          CALL WRITBL
          WRITE(LUPR1,*) 'Crek calculation ended at slab = ',IZ
          CALL WRITBL
          CALL WRITST
          IF(IZ.EQ.NZ) CRKCAL=.FALSE.
        ENDIF
C**End of equilibrium calculation*******************************
C**Calculation of NOX Zeldovich reaction sources.
        CALL SUB3(L0XO,L0F(JXO),L0XN,L0F(JXN),L0XOH,L0F(JXOH))
        CALL SUB3(L0XH,L0F(JXH),L0XO2,L0F(JXO2),L0XNO,L0F(JXNO))
        CALL SUB3(L0FEP2,L0SPR2,L0FEP3,L0SPR3,L0FEP4,L0SPR4)
        L0RHO1=L0F(DEN1)
        IF(LSG62) THEN
          L0YN2=L0F(JYN2)
          L0VOL=L0F(VOL)
        ELSEIF(LSG63) THEN
          L0FYPR=L0F(JYPR)
        ENDIF
C
       if(dbgrnd) then
         call writ40('NOX calculation starts for current slab ')
         call writ2i('isweep  ',isweep,'iz      ',iz)
         call writbl
       endif
       DO 1301 JX=1,NX
         DO 1301 JY=1,NY
           J=JY+NY*(JX-1)
c Initialize arrays
           F(L0FEP2+J) = 0.0
           F(L0FEP3+J) = 0.0
           F(L0FEP4+J) = 0.0
c-- Store mechanism values in GFSN array
c-- Convert mass fractions to mole fractions ( kg-mol/Kg )
c H mole fraction
          GFSN(1)=F(L0XH+J)/GHM
c N mole fraction
          GFSN(2)=F(L0XN+J)/GNM
c NO mole fraction
          GFSN(3)=F(L0XNO+J)/GNOM
c N2 mole fraction
          IF(LSG62) THEN
            GFSN(4)=F(L0YN2 +J)/GN2M
          ELSEIF(LSG63) THEN
            GFSN(4)=F(L0FYPR +J)*GN2PR/GN2M
          ENDIF
c O mole fraction
          GFSN(5)=F(L0XO +J)/GOM
c OH mole fraction
          GFSN(6)=F(L0XOH+J)/GOHM
c O2 mole fraction
          GFSN(7)=F(L0XO2+J)/GO2M
c
c-- For blocked cells, skip the following.
          IF(.NOT.SLD(J)) THEN
CRJ            IF(F(L0VOL+J).GT.TINY) THEN
            GTK=F((L0F(JCRKT))+J)+TINY
            TKINV=1.0/GTK
            GRHO2=F(L0RHO1+J)*F(L0RHO1+J)
c
c-- Zeldovich mechanism - source terms
            DO 1303 IS=1,NS
              GFSN(IS)=AMAX1( TINY , GFSN(IS) )
              RSOUR(IS) = 0.0
              RSPR(IS) = 0.0
 1303       CONTINUE
c
            DO 1304 IR=1,NR
c-- Forward & reverse reaction rate calculations
              TK1=AMAX1(AMIN1((TACT1(IR)*TKINV),100.),-100.)
              TK2=AMAX1(AMIN1((TACT2(IR)*TKINV),100.),-100.)
c-- Forward & backward rates in kg-mol/s/m**3
              FRATE=GBX1(IR)*EXP(-TK1)*GTK**TEN1(IR)*GRHO2
              RRATE=GBX2(IR)*EXP(-TK2)*GTK**TEN2(IR)*GRHO2
c-- Species selection for reaction, ir
              J1=ID(1,IR)
              J2=ID(2,IR)
              J3=ID(3,IR)
              J4=ID(4,IR)
              FF=-FRATE*GFSN(J1)*GFSN(J2)+RRATE*GFSN(J3)
     1                   *GFSN(J4)
              DSDF=-FRATE*(GFSN(J1)+GFSN(J2))-RRATE*(GFSN(J3)
     1                     +GFSN(J4))
c-- Store source for N & NO in array RSOUR
              RSOUR(J1)=RSOUR(J1)+FF
              RSOUR(J2)=RSOUR(J2)+FF
              RSOUR(J3)=RSOUR(J3)-FF
              RSOUR(J4)=RSOUR(J4)-FF
              RSPR(J1)=RSPR(J1)+DSDF
              RSPR(J2)=RSPR(J2)+DSDF
              RSPR(J3)=RSPR(J3)+DSDF
              RSPR(J4)=RSPR(J4)+DSDF
1304        CONTINUE
          ENDIF
c-- Coefficient for N
          F(L0FEP2+J)=-RSPR(2)
c-- Value for XN ( multiply by molecular weight so that final
c                 reaction rate source is in units of kg/m**3/s )
          F(L0FEP4+J)=-(RSOUR(2)/(RSPR(2)+TINY))*GNM+F(L0XN+J)
c-- Value for XNO ( multiply by molecular weight so that final
c                 reaction-rate source is in units of kg/m**3/s )
          F(L0FEP3+J)=RSOUR(3)*GNOM
1301    CONTINUE
C***************************************************************
C--- GROUP 19. SECTION 6 ---- Finish of IZ slab.
C***************************************************************
      ELSEIF(IGR.EQ.19.AND.ISC.EQ.6) THEN
C
        IF(LSG62) THEN
C-- Equivalence ratio calculation
          IF(ISWEEP.EQ.LSWEEP-1.OR.ENUFSW) THEN
            IF(JEQUI.NE.0) THEN
            CALL SUB4(L0YCO,L0F(JYCO),L0YO2,L0F(JYO2),L0YH2,
     1                L0F(JYH2),L0YCO2,L0F(JYCO2))
            CALL SUB2(L0YH2O,L0F(JYH2O),L0EQUI,L0F(JEQUI))
C-- Valence (C-4,O-2,H-1) / molecular mass; the division by mol
C   mass is to convert mass frac in the loops below into moles.
            CALL SUB4R(GE1,4./28.,GE2,4./44.,GE3,2./18.,GE4,2./2.)
            CALL SUB4R(GE5,2./28.,GE6,4./44.,GE7,2./18.,GE8,4./32.)
C-- Equivalence ratio = ((4*moles C)+(1*moles H)) / (2*moles O)
            DO 1961 J=1,NX*NY
              GNUM=GE1*F(L0YCO+J)+GE2*F(L0YCO2+J)+GE3*F(L0YH2O+J)+
     1                                            GE4*F(L0YH2+J)
              GDEN=GE5*F(L0YCO+J)+GE6*F(L0YCO2+J)+GE7*F(L0YH2O+J)+
     1                                            GE8*F(L0YO2+J)
              GEQUIV=GNUM/(GDEN+TINY)
              F(L0EQUI+J)=GEQUIV
 1961       CONTINUE
          ENDIF
          IF(JDEGF.NE.0) THEN
C-- Calculate temperatures in Fahrenheit
            L0DEGF=L0F(JDEGF)
C-- Convert crkt temperatures to fahrenheit
            DO 1962 J=1,NX*NY
 1962         F(L0DEGF+J)=(F((L0F(JCRKT))+J)-273.)*1.8+32.
          ENDIF
        ENDIF
C
C--  Calculate nox outlet data
        CALL GETPTC('OUTLET',GTYPE,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
        IF(IZ.LT.JZF .OR. IZ.GT.JZL) RETURN
        IF(IZ.EQ.JZF) THEN
          CALL SUB3R (GFLOT,0.0,GFLNT,0.0,GFNOT,0.0)
          CALL SUB4R(GCOM,28.01055,GO2M,31.9988,GN2M,28.0134,
     1               CO2M,44.00995)
          CALL SUB3R(GH2M,2.01594,GH2OM,18.01534,GNOM,30.0061)
          GFNOVT=0.
        ENDIF
        CALL SUB4(L0P1,L0F(P1),L0R1,L0F(R1),L0RHO,L0F(DEN1),
     1            L0XN,L0F(JXN))
        CALL SUB4(L0XNO,L0F(JXNO),L0YCO,L0F(JYCO),L0YO2,L0F(JYO2),
     1            L0YCO2,L0F(JYCO2))
        CALL SUB3(L0YN2,L0F(JYN2),L0YH2,L0F(JYH2),L0YH2O,L0F(JYH2O))
        IF(QEQ(GTYPE,6.0).OR.QEQ(GTYPE,7.0)) THEN
          L0AREA=L0F(AHIGH)
        ELSEIF (QEQ(GTYPE,4.0).OR.QEQ(GTYPE,5.0)) THEN
          L0AREA=L0F(ANORTH)
        ELSEIF (QEQ(GTYPE,2.0).OR.QEQ(GTYPE,3.0)) THEN
          L0AREA=L0F(AEAST)
        ELSE
          WRITE (LUPR1,*) 'Outlet boundary is not an AREA type.'
          WRITE (LUPR1,*) 'Outlet summary coding skipped.'
          RETURN
        ENDIF
        CALL GETCOV('OUTLET',P1,GCO,GVAL)
        IF(IZ.EQ.JZF) CALL SUB4R(GEXC,0.,GEXH,0.,GEXO,0.,GEXN,0.)
          DO 1963 JX=JXF,JXL
          DO 1963 JY=JYF,JYL
            J=JY+(JX-1)*NY
            GFLOW=-GCO*F(L0R1+J)*(GVAL-F(L0P1+J))*F(L0AREA+J)
            GFLONO=GFLOW*F(L0XNO+J)
            RMIXP=(F(L0YCO+J)/GCOM+F(L0YO2+J)/GO2M+F(L0YN2+J)/GN2M
     1                  +F(L0YCO2+J)/GCO2M+F(L0YH2+J)/GH2M
     1                  +F(L0YH2O+J)/GH2OM)*GASCON
            GFLNOV=GFLONO* (GASCON/(GNOM*RMIXP))
            GFLOT=GFLOT+GFLOW
            GFNOT= GFNOT + GFLONO
            GFNOVT= GFNOVT + GFLNOV
            GEXC=GEXC+(F(L0YCO+J)*12./28.+F(L0YCO2+J)*12./44.)*GFLOW
            GEXH=GEXH+(F(L0YH2+J)+F(L0YH2O+J)*2./18.)*GFLOW
            GEXO=GEXO+(F(L0YCO+J)*16./28.+F(L0YCO2+J)*32./44.
     1           +F(L0YH2O+J)*16./18.+F(L0YO2+J))*GFLOW
            GEXN=GEXN+(F(L0YN2+J))*GFLOW
 1963     CONTINUE
          IF(IZ.LT.JZL) RETURN
          GNOEX=1.E6*GFNOT/GFLOT
          GNOVX=1.E6*GFNOVT/GFLOT
          IF(ISWEEP.LT.LSWEEP.AND..NOT.ENUFSW) RETURN
          WRITE(LUPR1,*)
          WRITE(LUPR1,*)'      *******************************'
          WRITE(LUPR1,*)'      * NOX FORMATION : OUTLET DATA *'
          WRITE(LUPR1,*)'      *******************************'
          WRITE(LUPR1,1969) GFLOT,GFNOT,GNOEX,GNOVX,GEXC,GEXH,
     1                      GEXO,GEXN
          CALL WRITBL
1969      FORMAT(/6X,' Exit mass flow rate    =  ',1PE12.3,' kg/s',
     1       /6X,' ===================',
     1       /6X,' Exit NO mass flow rate =  ',1PE12.3,' kg/s'
     1       /6X,' ======================',
     1       /6X,' Exit NO mass fraction  =  ',1PE12.3,' ppm by mass'
     1       /6X,' ===================== ',
     1       /6X,' Exit NO vol fraction   =  ',1PE12.3,' ppm by vol',
     1       /6X,' ==================== ',
     1       /6X,' Exit C mass flow rate  =  ',1PE12.3,' kg/s'
     1       /6X,' =====================',
     1       /6X,' Exit H mass flow rate  =  ',1PE12.3,' kg/s'
     1       /6X,' =====================',
     1       /6X,' Exit O mass flow rate  =  ',1PE12.3,' kg/s'
     1       /6X,' =====================',
     1       /6X,' Exit N mass flow rate  =  ',1PE12.3,' kg/s'
     1       /6X,' =====================',//)
       ELSEIF(LSG63) THEN
C
         IF(IZ.EQ.NZ) THEN
           IF(ISWEEP.EQ.LSWEEP.OR.ENUFSW) THEN
C
             CALL SUB2(JFU,LBNAME('FUEL'),JMF,LBNAME('MIXF'))
             CALL SUB2(L0P1,L0F(P1),L0XN,L0F(JXN))
             CALL SUB2(L0XNO,L0F(JXNO),L0FU,L0F(JFU))
             L0MF=L0F(JMF)
C
             GCO = FIXP
             GVAL= 0.0
             CALL SUB3R(GFLOT,0.0,GFLOF,0.0,GFLOFU,0.0)
             CALL SUB2R(GFLON,0.0,GFLONO,0.0)
             DO 1965 J = 1 , NX*NY
               GFLOW   = -GCO*(GVAL-F(L0P1+J))
               GFLOT  = GFLOT + GFLOW
               GFLOF  = GFLOF  + GFLOW*F(L0MF+J)
               GFLOFU = GFLOFU + GFLOW*F(L0FU+J)
               GFLON  = GFLON  + GFLOW*F(L0XN+J)
               GFLONO = GFLONO + GFLOW*F(L0XNO+J)
1965         CONTINUE
             GFEX  = (GFLOF /(GFLOT+1.E-10))
             GFUEX = (GFLOFU/(GFLOT+1.E-10))
             GNEX = (GFLON /(GFLOT+1.E-10))*1.E6
             GNOEX = (GFLONO/(GFLOT+1E-10))*1.E6
             CALL WRITST
             WRITE(LUPR1,*)'* NOX FORMATION : OUTLET DATA *'
             WRITE(LUPR1,*)'*******************************'
             WRITE(LUPR1,*) '   Isweep = ',ISWEEP
             WRITE(LUPR1,1966) GFLOT,GFEX,GFUEX,GNEX,GNOEX
1966         FORMAT(/6X,' Exit mass flow rate =         ',
     1         1PE12.4,' kg/s',
     1         /6X,' ===================',/,
     1         /6X,' Exit mixture fraction =       ',1PE12.4,
     1         /6X,' =====================',/,
     1         /6X,' Exit fuel mass fraction =     ',1PE12.4,
     1         /6X,' =======================',/,
     1         /6X,' Exit  N mass fraction = ',1PE12.4,' ppm',
     1         /6X,' =======================',/,
     1         /6X,' Exit  NO mass fraction =      ',1PE12.4,' ppm',
     1         /6X,' ======================',//)
C
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      NAMSUB='gxnox'
      END
C----------------------------------------------------------------------
c
      SUBROUTINE GXCREK
C***********************************************************************
C.... This subroutine is called only from GXNOX in PHOENICS-3.2; but it
C     may be called from other subroutines if the user provides them.
C     The PHOENICS Encyclopaedia provides information about the function
C     of the CREK package generally.
C***********************************************************************
C
      COMMON/CEQUIL/AL(7,20),ATOM(2,7),B0(7),PI(7)
      COMMON/NCEQUL/CATOM(7)
      CHARACTER*1 CATOM
      COMMON/CINDEX/IHCPS,IMAT,ITER,JJ,N1,N2,N3,NA,NGLOB,NLM,NSM,NQ,
     1IDCO,IDCO2,IDH2,IDH2O,IDN2,IDO2
      COMMON/CKINET/CPSUM,ER,PPLN,FQ,RGAS,RGASIN,SMINV,TKINV,TLN
      COMMON/CPARAM/EMV,HSUB0,NS,PA,Q0,Q1,Q2,Q3,Q4,RHOP,SM,
     1S1(20),S2(20),TK
      LOGICAL LADIAB,LCONVG,LDEBUG,LEQUIL,LREACT,LNRG
      COMMON/CHMLGC/LADIAB,LCONVG,LDEBUG,LEQUIL,LREACT,LNRG
      COMMON/CASUB/ASUB(20,3)
      CHARACTER*4 ASUB
      COMMON/CSPECE/H0(20),S0(20),SMW(20),SSAVE(20),Z(7,2,20)
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      CHARACTER*1 CARB,HYDR
      LOGICAL EQZ,QEQ
      DATA FACTOR/5.0/,PATM/101325.0/,SMALL/1.0E-6/
      DATA CARB/'C'/,HYDR/'H'/
C
C     NORMAL SOLUTION
C
C.... Determine equivalence ratio; if outside interval (.1,10), assume
C     no reaction and 230694 adiabatic non-reacted mixture properties;
C     save given estimates or program generated estimates if tk is small.
C     if solution is successful, 230694 to calling program.  otherwise,
C     enter problem cell logic below
C
c      CALL WRITST
c      CALL WRITBL
c      CALL WRIT40('Entering subroutine GXCREK, from GXNOX  ')
c      CALL WRITBL
c      CALL WRITST
      NAMSUB='GXCREK'
      CALL ERATIO
      IF(.NOT.LREACT) GO TO 400
C.... Allow combustion for all mixtures....
      EMVSAV=EMV
      PPLN=ALOG(PA/PATM)
      TSAVE=TK
      DO 10 I=1,NS
   10 SSAVE(I)=S2(I)
      IF(TK.LT.SMALL) GO TO 100
C
      CALL SPECE
C
      IF(LCONVG) GO TO 900
C
C     PROBLEM CELL
C
C.... Solution logic is different for four types of problems as follows:
C     MODE 1 ... LEQUIL = T, LADIAB = T
C     MODE 2 ... LEQUIL = T, LADIAB = F
C     MODE 3 ... LEQUIL = F, LADIAB = T
C     MODE 4 ... LEQUIL = F, LADIAB = F
C
C.... Always try rt=0.0 (lnrg=f) after solution failure when lequil=f.
C     logic to find solution is controlled in chapters 1 and 2 below
C     where in each section, new estimates are determined either by
C     saved given ones, new assigned ones, or solution found not at
C     required conditions. The variables nextok and nextng are assigned
C     the statement numbers of where to go if the solution attempt is
C     successful or not, respectively.
C
      ASSIGN 900 TO NEXTOK
      ASSIGN 100 TO NEXTNG
      GO TO 520
C
  100 MODE=4
      IF(LEQUIL) MODE=MODE-2
      IF(LADIAB) MODE=MODE-1
C
C* * * * * * * * * * * * * * * * * * * * * * * * CHAPTER 1 * * * * * * *
C
C     EQUILIBRIUM
C.... This chapter makes equilibrium estimates amd initiates strategy for
C     cases in which convergence was not achieved on first call to spece.
C
      LEQUIL=.TRUE.
      LADIAB=.TRUE.
      IF(MODE.LT.3) GO TO 130
C
C.... First use given estimates for equil soln in mode 3 and 4 problems
C
      TK=TSAVE
      DO 120 I=1,NS
  120   S2(I)=SSAVE(I)
      ASSIGN 200 TO NEXTOK
      ASSIGN 130 TO NEXTNG
      GO TO 500
C
C.... Complete combustion estimate
C
  130 IF(IDCO2.EQ.0) GO TO 170
C.... If no carbon in fuel, skip to garbage estimates
      X=0.0
      Y=0.0
      DO 132 I=1,NS
        S2(I)=SMALL
        DO 131 L=1,NLM
          IF(EQZ(AL(L,I))) GO TO 131
          IF(CATOM(L).EQ.CARB) X=X+AL(L,I)*S1(I)
          IF(CATOM(L).EQ.HYDR) Y=Y+AL(L,I)*S1(I)
  131   CONTINUE
  132 CONTINUE
      S2(IDN2)=S1(IDN2)
C
      ER1=(4.0*X+Y)/(2.0*X+Y)
      ER2=2.0+Y/(2.0*X+1.e-20)
      IF(ER.LE.1.0) GO TO 133
      IF(ER.LT.ER1) GO TO 134
      IF(ER.LT.ER2) GO TO 135
C.... ER.gt.ER2 ---> all c and h in co, h2 and unburned cxhy
      S2(IDCO)=2.0*(X+Y/4.0)
      S2(IDH2)=Y*(1.0+Y/(4.0*X))
C.... Mole number for cxhy should be (er-2*(1+y/(4*x)))
      GO TO 136
C.... Fuel-lean combustion ---> only co2 and h2o formed
  133 S2(IDCO2)=X*ER
      S2(IDH2O)=Y*ER/2.0
      S2(IDO2)=(1.0-ER)*(X+Y/4.0)
      GO TO 136
C.... Sightly fuel-rich ---> co,co2, and h2o present
  134 S2(IDCO)=2.0*(X+Y/4.0)*(ER-1.0)
      S2(IDCO2)=Y*(1.0-ER)/2.0+X*(2.0-ER)
      S2(IDH2O)=Y*ER/2.0
      GO TO 136
C.... Fuel-rich ---> co,h2, and h2o present
  135 S2(IDCO)=ER*X
      S2(IDH2)=Y*(ER-1.0)/2.0+X*(ER-2.0)
      S2(IDH2O)=Y/2.0+X*(2.0-ER)
  136 SM=0.0
      DO 137 I=1,NS
  137 SM=SM+S2(I)
      TK=1500.0
      TKINV=6.6666666E-4
      IHCPS=1
      XLO=TK
      DO 139 K=1,10
        CALL HCPS
        HSUM=0.0
        DO 138 I=1,NS
  138   HSUM=HSUM+H0(I)*S2(I)
        TK=TK*(1.0+(HSUB0*RGASIN*TKINV-HSUM)/CPSUM)
        TKINV=1.0/TK
        XHI=ABS(TK-XLO)
        XLO=TK
        IF(XHI.LT.1.0) GO TO 141
  139 CONTINUE
      WRITE(14,140) K,XHI
  140 FORMAT(/10X,'MIXTURE TEMPERATURE FAILED TO CONVERGE IN  ',
     1    I2,' ITERATIONS     TEMP DIFFERENCES =',1PE12.3/)
  141 CONTINUE
C
      IF(LDEBUG) WRITE(14,142) (S2(K),K=1,NS),SM,TK
  142 FORMAT(/5X,'COMPLETE COMBUSTION ESTIMATE',/1P11E10.2/1P11E10.2/)
      IF(MODE.EQ.1) ASSIGN 900 TO NEXTOK
      IF(MODE.EQ.2) ASSIGN 300 TO NEXTOK
      IF(MODE.GE.3) ASSIGN 200 TO NEXTOK
      ASSIGN 170 TO NEXTNG
      GO TO 500
C
C.... Garbage estimates  (gordon and mcbride)
C
  170 TK=3800.0
      SM=0.1/FLOAT(NS)
      DO 171 I=1,NS
  171 S2(I)=SM
      SM=0.1
      IF(MODE.EQ.1) ASSIGN 900 TO NEXTOK
      IF(MODE.EQ.2) ASSIGN 300 TO NEXTOK
      IF(MODE.GE.3) ASSIGN 200 TO NEXTOK
      ASSIGN 600 TO NEXTNG
      GO TO 500
C
C**  **  **  **  **  **  **  **  **  **  **  **  CHAPTER 2  **  **  **
C
C     KINETIC
C.... Section for kinetic solution from adiabatic equilibrium estimates
C     (mode 3 and 4 only)
C
C.... Near-equilibrium solution  (kinetic with emv=1.0e-3)
C
  200 LEQUIL=.FALSE.
      IX=0
      EMV=1.0E-3
      XLO=EMV
C.... Increase minor species from equilibrium estimates
      DO 201 I=1,NS
        IF(S2(I).LT.SMALL) S2(I)=SMALL
  201 CONTINUE
      ASSIGN 230 TO NEXTOK
      ASSIGN 210 TO NEXTNG
      GO TO 500
C.... Failure on near-equil with emv=xlo, decrease emv by an order of
C     magnitude and attempt again, iterating this way up to 12 times.
  210 EMV=EMV*0.1
      XLO=EMV
      IX=IX+1
      IF(IX.EQ.12) GO TO 610
      TK=TSAVE
      DO 211 I=1,NS
  211 S2(I)=SSAVE(I)
      ASSIGN 230 TO NEXTOK
      ASSIGN 210 TO NEXTNG
      GO TO 500
C.... Have near-equil solution, so first try directly to obtain
C     required solution at given emv
  230 EMV=EMVSAV
      IF(MODE.EQ.3) ASSIGN 900 TO NEXTOK
      IF(MODE.EQ.4) ASSIGN 300 TO NEXTOK
      ASSIGN 250 TO NEXTNG
      GO TO 500
C
C     UPPER BRANCH MARCHING
C.... Have a kinetic solution but at emv .lt. emvsv.  start at known
C     solution and increase emv by factor to move towards a soln there,
C     if successful, repeat until emvsav is reached; if not successful
C     start half interval searching described below
C
  250 XLO=EMV
      EMV=XLO*FACTOR
      IF(EMV.GT.EMVSAV) EMV=EMVSAV
      XHI=EMV
      IX=0
      TK=TSAVE
      DO 251 I=1,NS
  251 S2(I)=SSAVE(I)
      ASSIGN 250 TO NEXTOK
      ASSIGN 270 TO NEXTNG
      GO TO 500
C
C     HALF-INTERVAL SEARCHING
C.... Have solution at xlo but not at xhi, hence start interval
C     searching by setting emv to the logarithmic average;
C     if iterating more than ten times, terminate.
C
  270 IX=IX+1
      TK=TSAVE
      DO 271 I=1,NS
  271 S2(I)=SSAVE(I)
      IF(IX.GT.10) GO TO 620
      EMV=SQRT(XLO*XHI)
      XHI=EMV
      ASSIGN 250 TO NEXTOK
      ASSIGN 270 TO NEXTNG
      GO TO 500
C
C***   ***   ***   ***   ***   ***   ***   ***   CHAPTER 3   ***   ***
C
C     NON-ADIABATIC
C.... Section for non-adiabatic solutions from adiabatic estimates
C     (mode 2 and 4 only)
C     Try directly to obtain non-adiabatic solution; if not successful,
C     start half-interval scaling from the adiabatic solution by
C     defining a scaling factor fq (0.0-1.0) to multiply the non-adiabatic
C     term (q) in the energy equation in spece
C
  300 LADIAB=.FALSE.
      XLO=0.0
      XHI=1.0
      FQ=1.0
      IX=0
  310 ASSIGN 320 TO NEXTOK
      ASSIGN 330 TO NEXTNG
      GO TO 500
C
  320 IF(QEQ(FQ,1.0)) GO TO 900
      FQ=1.0
      XLO=FQ
      IX=0
      GO TO 310
C
  330 IX=IX+1
      IF(IX.GT.10) GO TO 340
      TK=TSAVE
      DO 331 I=1,NS
  331 S2(I)=SSAVE(I)
      XHI=FQ
      FQ=0.5*(XLO+XHI)
      GO TO 310
C
  340 FQ=1.0
      GO TO 630
C
C****    ****    ****    ****    ****    ****    CHAPTER 4   ****    ***
C
C     FAILURE EXITS
C.... Failed equil or kinetic soln or equiv ratio outside (0.1,10)
C     return adiabatic, non-reacted mixture properties
C
  400 SM=0.0
      DO 401 I=1,NS
        S2(I)=S1(I)
  401 SM=SM+S2(I)
      TK=1000.0
      XLO=TK
      IHCPS=1
      TKINV=1.0E-3
      DO 403 I=1,10
        CALL HCPS
        HSUM=0.0
        DO 402 K=1,NS
          HSUM=HSUM+H0(K)*S2(K)
  402   CONTINUE
        TK=TK*(1.0+(HSUB0*RGASIN*TKINV-HSUM)/CPSUM)
        TKINV=1.0/TK
        XHI=ABS(TK-XLO)
        XLO=TK
        IF(XHI.LT.1.0) GO TO 404
  403 CONTINUE
      WRITE(14,140) I,XHI
  404 CONTINUE
      RHOP=PA/(RGAS*TK*SM)
      GO TO 900
C
C*****     *****     *****     *****     *****   CHAPTER 5   *****     *
C
C     PROBLEM CELL CALL TO SPECE
C.... Take the estimates generated in chapters 1,2 and attempt a solution
C     with full equations. If successful, update the save answers with the
C     solution and return to statement number nextok. If not, the action
C     depends on whether an equilbm or kinetic soln is sought.  Failed
C     equil soln, return to statement number nextng, while failure in a
C     kinetic soln will be followed by an attempt with lnrg=f ---> rt=0.0
C     and same estimates.  Setting rt=0.0 implies that a change in temp
C     field has no effect on species distribution for that particular
C     iteration, but does allow the species changes to influence the temp
C     change ---> partial decoupling of the energy equation. If successful
C     with lnrg=f, repeat with lnrg=t (full equations); if still no good,
C     return to statement number nextng.
C
  500 CALL SPECE
      IF(LCONVG) GO TO 540
C     SOLUTION FAILED TRY RT=0.0
      IF(LNRG) GO TO 520
      LNRG=.TRUE.
  510 GO TO NEXTNG, (100,130,170,210,250,270,330,600)
  520 IF(LEQUIL) GO TO 510
      LNRG=.FALSE.
C**MPD 530 LNRG=.FALSE.     >>>> UNREFERENCED STATEMENT LABEL <<<< ??
      GO TO 500
C     HAVE SOLN BUT WITHOUT RT TERMS, SEND THIS
C     SOLN AS ESTIMATES WITH RT CALCULATION
  540 IF(LNRG) GO TO 550
      LNRG=.TRUE.
      GO TO 500
C
C.... Solution is successful, update save answers and continue.
C
  550 TSAVE=TK
      DO 560 I=1,NS
  560 SSAVE(I)=S2(I)
C
      GO TO NEXTOK, (200,230,300,320,900)
C
C******      ******      ******      ******      CHAPTER 6    ******
C
C    ERROR MESSAGES
C
  600 CONTINUE
C     WRITE(14,601)
C 601 FORMAT(/10X,3(4H****),' FAILURE TO FIND EQUIL SOLN...',
C    1                      'AVG INLET PROPS RETURNED'/)
C     RETURN INITIAL GUESS....NOT THE UNREACTED MIXTURE PROPERTIES
      GO TO 670
  610 WRITE(14,611)
  611 FORMAT(/10X,3(4H****),' FAILURE TO FIND NEAR-EQUIL SOLN...',
     1                      'AVG INLET PROPS RETURNED'/)
      GO TO 650
  620 WRITE(14,621)
  621 FORMAT(/10X,3(4H****),' FAILURE TO OBTAIN KINETIC SOLN AFTER',
     1            ' TEN INTERVAL HALVING...AVG INLET PROPS RETURNED'/)
      GO TO 650
  630 WRITE(14,631)
  631 FORMAT(/10X,3(4H****),' NON-ADIABATIC SOLN FAILED...',
     1                      'ADIAB SOLN RETURNED'/)
      GO TO 670
C
C.... Restore failed problem mode prior to return
C
  650 IF(MODE.EQ.2) LADIAB=.FALSE.
      IF(MODE.EQ.3) LEQUIL=.FALSE.
      IF(MODE.EQ.4) LADIAB=.FALSE.
      GO TO 400
C
C.... Failed non-adiabatic solution...return adiabatic
c     equil or kinetic solution
C
  670 TK=TSAVE
      DO 671 I=1,NS
  671 S2(I)=SSAVE(I)
C
  900 NAMSUB='gxcrek'
      END
c