cgxchki.htm

C    element mass fraction and mole fraction
**********************************************************************
c
c
      SUBROUTINE GXCHKI(ICELLA,VALOUT)
C-----------------------------------------------------------------
C                    PHOENICS-CHEMKIN INTERFACE
C-----------------------------------------------------------------
C   Function:
C     The interface uses CHEMKIN and TRANLIB routines to calculate
C     thermodynamic and transport properties for gas mixtures, and
C     calculates the diffusion fluxes according to 2 models of
C     diffusion. The diffusion models are:
C       0) Standard PHOENICS model - ENULA = 0.0  : IMCOPT = 0 (default)
C       1) Mixture-averaged model  - ENULA = GRND9: IMCOPT = 1
C     Furthermore, the influence of chemical reactions is treated
C     according to 2 different algorithms:
C       1) Provide source terms calculated using CHEMKIN routines
C          and solve using the standard PHOENICS algorithms
C       2) Solve the equations using the PBP algorithm; but use
C          a Newton-Raphson algorithm operating simultaneously on
C          all species and on temperature (enthalpy) to generate
C          the corrections
C     The source-term options are indicated to the code by way of
C     the PATCH name:
C          Characters 1-5 must be CHEMK
C          Character 6 must be 1 to indicate which of the
C          options is to be used
C          Characters 7-8 are free to be used by the user
C
C   Assumptions:
C     Entered only for 1st phase
C     Only mixture-averaged model for transport (IMCOPT==0)
C     Only solution for TEM1 is implemented (for energy equ'n.)
C     Properties are calculated for every sweep
C     Non-orthogonality terms have not been included in the enhanced
C     diffusion model
C
C     Notes: The implementation of the thermal diffusion terms in the
C            species diffusion equations was modified on 21/10/03 so
C            ensure species mass conservation. The previous treatment
C            was discovered to be non-conservative when running library
C            cases C206 to C208 inclusive.
C
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C-----------------------------------------------------------------
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/pltcfile'
      INCLUDE '/phoenics/d_includ/parear'
      COMMON/INDAUX/L0ISL,L0IST,L0SL,L0ST,NSTO,NSOL,L0SLRS,L0TTRS,
     1                                                           IFL(12)
C
      COMMON/GENI/NXNY,IGFL1,NXNYST,IGFL4,KDUMM,IGFL5(4),NFM,NWHOLE,
     1            IGSH,IGFL13(30),NFTOT,IGFL44(6),ITEM,IGFL51(5),
     1            IPRPS,IGFL57(4)
      COMMON/INDUN/INDFL(12),INN,INS,INE,INW,INH,INL
      COMMON/IDA7/LOFLUX(150) /IDA8/LODFCO(150)
      COMMON/SLDPRP/SOLPRP,PORPRP,VACPRP
      COMMON/F01/M0F(600)
C
      CHARACTER*256 NKEEP, NTEMP
C
      PARAMETER (NCKMI0=84)
      PARAMETER (VISSCA=1.E-1,CONSCA=1.E-5,DIFSCA=1.E-4)
C
      INTEGER LBCKF, ICKMC0(NCKMI0), PVCKID
      LOGICAL LCKERR,GRN,EQZ,VARPRS,QEQ,NEZ,SOLTEM,HRM,
     1        DBGLOC,QGE,VARTEM,CONTEM,ERROR,TSTVPO,BLK,BYPASS,
     2        BLOCKD,SOLID,FIXCEL,MONPRT,DENSC,CFIRST,TPSEM
      CHARACTER*9  CODE
      CHARACTER*40 PRECIS
      CHARACTER*4  NAMD,CMFILE
C
      EQUIVALENCE (PREF,CHSOC),
     3            (ICKMC0(1),LBCKF),(PVCKID,LBPVSP(6)),
     4            (CMFILE,CSG4)
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*****double precision
      DOUBLE PRECISION
     1                 RU,RUC,PRES,CPMIX,RHOMIX,VISMIX,
     2                 CONMIX,VOLUME,ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,
     3                 UFAC,DFAC,U0,U,RTMP,QDOT,DTTP,DTMIN,DTMAX,
     4                 DTEMP0,GRHO
      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
     2       /CKMCR0/ATIM,RTIM,ATOL,RTOL,UFAC,DFAC,ABSOL,RELAT,
     2               RU,RUC,DTMIN,DTMAX,DTEMP0
     4       /CKMCL0/VARPRS,CONTEM,SIUNIT
      SAVE /CKMCI0/,/CKMCL0/,
     1     DODIFF,DOSORC,DOTDIF,SOLTEM,MODDIF,VARTEM,ENTFLX,RELQDT,
     2     ALQDOT,STOCMP,BYPASS,DENSC,CFIRST,TPSEM
C
      COMMON /CKMIND/ NDCMI,NDCMR,NDCMC,NDCML
      LOGICAL DODIFF,DOSORC,DOTDIF,DFN,WHF,MODDIF,ENTFLX,RELQDT,
     1        STOCMP,SIUNIT
      logical exists
      character*(1024) cval, delim*1
C*****double precision
      DATA PRECIS /' DOUBLE PRECISION Version is activated  '/
C*****END double precision
C
      dbgloc=debug.and.izstep>=izdb1 .and.izstep<=izdb2.and.
     1                   isweep>=iswdb1.and.isweep<=iswdb2
      IXL=IABS(IXL)
      NXNY=NX*NY
C
      IF(IGR==1) CFIRST=.TRUE.
      IF(IGR>1) THEN
        IF(VARPRS) L0P1=L0F(P1)
        IF(ICKOP1>1.AND.STORE(KE).AND.STORE(EP))
     1               CALL SUB2(L0KE,L0F(KE),L0EP,L0F(EP))
        IF(STORE(ITIND)) THEN
          L0TEM=L0F(ITIND)
        ELSE
          L0TEM=0
          VARTEM=.FALSE.
        ENDIF
        L0DEN=L0F(DEN1)
        IF(IQDOT>0) THEN
          L0QDOT=L0F(IQDOT)
          ALQDOT=DTFALS(IQDOT)
          RELQDT=QGE(ALQDOT,-1.0).AND.ALQDOT<0.0
          ALQDOT=-DTFALS(IQDOT)
        ELSE
          L0QDOT=0
        ENDIF
        IF(IENT>0) THEN
          L0ENT=L0F(IENT)
        ELSE
          L0ENT=0
        ENDIF
        IF(ICON>0) THEN
          L0CON=L0F(ICON)
        ELSE
          L0CON=0
        ENDIF
        IF(IMMWT>0) THEN
          L0MMWT=L0F(IMMWT)
        ELSE
          L0MMWT=0
        ENDIF
        IF(STORE(VPOR)) THEN
          L0VPOR=L0F(VPOR)
          TSTVPO=BLK(IZSTEP,7)
        ELSE
          L0VPOR=0
          TSTVPO=.FALSE.
        ENDIF
        IF(IPRPS==0) THEN
          L0PRPS=0
        ELSE
          L0PRPS=L0F(IPRPS)
        ENDIF
      ENDIF
C
      IF(IGR==8.OR.IGR==9.OR.IGR==13.OR.IGR==11
     1           .OR.IGR==19) THEN
        IF(WHF(LBCKF)) THEN
          IZOFF=(IZSTEP-1)*NXNY
        ELSE
          IZOFF=0
        ENDIF
        IF(STOCMP) THEN
          DO 10002 K=1,KK
            ICKMC(K0L0+K)=L0F(LBCKF+K-1)
10002     CONTINUE
        ENDIF
      ENDIF
      GO TO (1,25,25,25,25,25,25,8,9,25,11,25,13,25,25,25,25,25,19,
     1       25,25,25,23,25),IGR
   25 CONTINUE
      go to 999
C-----------------------------------------------------------------
C--- GROUP 1. Run title and other preliminaries
C
    1 CONTINUE
C
C     User may here change message transmitted to the VDU screen
      IF(ISC/=1) GO TO 999
        IF(.NOT.NULLPR) THEN
          CALL WRITBL
          CALL WRITST
          CALL WRYT40('GXCHKI file is GXCHKI.HTM of:    240613 ')
          CALL WRYT40(' PHOENICS-CHEMKIN Interface v1.0        ')
          CALL WRYT40(PRECIS)
          CALL WRIT40(' PHOENICS-CHEMKIN Interface v1.0        ')
          CALL WRIT40(PRECIS)
          CALL WRITST
        ENDIF
C
C--- Initialise /CKMCI0/
C
      DO 10000 I=1,NCKMI0
        ICKMC0(I)=0
10000 CONTINUE
C
      LBCKF=16
      LBCKF0=NINT(CHSOB)
      IF(LBCKF0/=0) LBCKF=LBCKF0
      STOCMP=LBCKF>0
      DTEMP0=TEMP0
C
      IF(STOCMP.AND.GRNDNO(9,CHSOA)) THEN
        ICKOPT=1
        IF(SELREF) WRITE(LUPR1,*)
     1    'SELREF has been set false in GXCHKI because CHSOA=GRND9'
        SELREF=.FALSE.
      ELSE
        ICKOPT=0
      ENDIF
C
C--- Use GETSPD to retrieve the Q1 settings which overwrite
C    the default parameters
C
      CALL GETSPD('CHEM','LEVEL1',2,RDUM,LEVEL1,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','LEVEL2',2,RDUM,LEVEL2,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','NJAC'  ,2,RDUM,NJAC  ,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','ITJAC' ,2,RDUM,ITJAC ,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','IRETIR',2,RDUM,IRETIR,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','NUMDT' ,2,RDUM,NUMDT ,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','NINIT' ,2,RDUM,NINIT ,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','CALL1' ,2,RDUM,ICALL1,.FALSE.,'   ',IERR)
      VARPRS=.FALSE.
      CALL GETSPD('CHEM','VARPRS',3,RDUM,IDUM,VARPRS,'   ',IERR)
      CONTEM=.FALSE.
      CALL GETSPD('CHEM','CONTEM',3,RDUM,IDUM,CONTEM,'   ',IERR)
      SIUNIT=.FALSE.
      CALL GETSPD('CHEM','SIUNIT',3,RDUM,IDUM,SIUNIT,'   ',IERR)
C   option to bypass twopnt-solver termination
      BYPASS=.FALSE.
      CALL GETSPD('CHEM','BYPASS',3,RDUM,IDUM,BYPASS,'   ',IERR)
C   Provision for density conversion from SI to CGS
      DENSC=.FALSE.
      CALL GETSPD('CHEM','DENSC',3,RDUM,IDUM,DENSC,'   ',IERR)
      TPSEM=.TRUE.
      CALL GETSPD('CHEM','TPSEM',3,RDUM,IDUM,TPSEM,'   ',IERR)
      CALL SUB4R(UFAC1,0.0,DFAC1,0.0,DTMIN1,0.0,DTMAX1,0.0)
      CALL GETSPD('CHEM','UFAC' ,1,UFAC1 ,IDUM,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','DFAC' ,1,DFAC1 ,IDUM,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','DTMIN',1,DTMIN1,IDUM,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','DTMAX',1,DTMAX1,IDUM,.FALSE.,'   ',IERR)
      RTOL=1.E-4
      ATOL=1.E-8
      ATIM=3.*ATOL
      RTIM=3.*RTOL
      CALL SUB4R(ATOL1,0.0,RTOL1,0.0,ATIM1,0.0,RTIM1,0.0)
      CALL GETSPD('CHEM','ATOL',1,ATOL1 ,IDUM,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','RTOL',1,RTOL1 ,IDUM,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','ATIN',1,ATIM1,IDUM,.FALSE.,'   ',IERR)
      CALL GETSPD('CHEM','RTIM',1,RTIM1,IDUM,.FALSE.,'   ',IERR)
      IF(NEZ(ATOL1)) ATOL=ATOL1
      IF(NEZ(RTOL1)) RTOL=RTOL1
      IF(NEZ(ATIM1)) ATIM=ATIM1
      IF(NEZ(RTIM1)) RTIM=RTIM1
C
C--- Make some debug settings
C
      IF(ICKOPT>=1) THEN
        CALL SUB2(LEVEL(1),LEVEL1,LEVEL(2),LEVEL2)
        LEVEL(1)=MAX(LEVEL(1),LEVEL(2))
      ENDIF
C
C--- Are there any source PATCHes, are there any dangerous CHSO
C    PATCHes?
C
      DOSORC=ICKOPT>0
      LCKERR=.FALSE.
      DO 10013 IP=1,NUMPAT
        IF(.NOT.DOSORC.AND.NAMPAT(IP)(1:5)=='CHEMK') THEN
          DOSORC=.TRUE.
          READ(NAMPAT(IP)(6:6),'(I1)') ICKOP1
        ELSE
          LCKERR=LCKERR.OR.NAMPAT(IP)(1:4)=='CHSO'
        ENDIF
10013 CONTINUE
C
C--- If CHSO has been set we have to exit
C
      IF(LCKERR) THEN
        CALL WRIT40(' GXCHKI - CHSO PATCHES CANNOT BE USED   ')
        CALL WRIT40('          WITH THE CHEMKIN INTERFACE    ')
        CALL SET_ERR(254, 'Error. See result file.',1)
C        CALL EAROUT(1)
      ENDIF
C
C--- Set a pointer to various variables
C
        IF(STORE(ITEM)) THEN
          ITIND=ITEM
        ELSEIF(STORE(TEMP1)) THEN
          ITIND=TEMP1
        ELSE
          ITIND=0
        ENDIF
        IMIXF=LBNAME('MIXF')
        IF(.NOT.STORE(IMIXF)) IMIXF=0
        SOLTEM=SOLVE(ITIND)
        VARPRS=VARPRS.AND.SOLVE(P1)
        IQDOT=LBNAME('QDOT')
        IF(.NOT.STORE(IQDOT)) IQDOT=0
        IENT =LBNAME('ENTH')
        IF(.NOT.STORE(IENT))  IENT=0
        ICON =LBNAME('COND')
        IF(.NOT.STORE(ICON))  ICON=0
        IMMWT =LBNAME('MMWT')
        IF(.NOT.STORE(IMMWT)) IMMWT=0
C
C--- Initialise CHEMKIN database (name may be constructed from CMFILE)
C
      NTEMP=' '
      DO 10201 ITRY=0,100,100
        JTRY=ITRY
        IF(CMFILE/=' ') THEN
          I=INDEX(CMFILE,' ')
          IF(I==0) THEN
            I=4
          ELSE
            I=I-1
          ENDIF
          IF(JTRY==0) THEN
            NKEEP=NMFIL(55)
            NMFIL(55)=CMFILE(1:I)//'ckln'
            CALL BLDNAM(NMFIL(55),NTEMP,LNMTMP)
          ELSE
            I55=IPEXT(55,JTRY/100,0)
            IC=INDEX(NMFLEX(I55),'cklink')
            IF(IC<=0.OR.IC>40) THEN
              IERR=1
              GO TO 10202
            ENDIF
            NKEEP=NMFLEX(I55)
            NMFLEX(I55)(IC:)=CMFILE(1:I)//'ckln'
            CALL BLDNAM(NMFLEX(I55),NTEMP,LNMTMP)
          ENDIF
        ENDIF
        IERR=-1
C#### JCL 24.06.13 copy file to slaves
        IF(MIMD.AND.NPROC>1.AND.JTRY==0)THEN
          CVAL=NMFIL(55)
          CALL PAROPENFL(90,CVAL,IERR,EXISTS)
          IF(EXISTS)THEN
            CALL GETSYSID(ISYSID)
            DELIM='/'
            IF(ISYSID<=2) DELIM=CHAR(92)
            CALL COPYBINARYFILE(DELIM,CVAL,IERR,IOS,EXISTS)
            IF(IERR/=0)GOTO 10201
          ELSE
            GOTO 10201
          ENDIF        
        ENDIF
        CALL OPENZZ(55+JTRY,IERR)
        IF(IERR==0) GO TO 10202
10201 CONTINUE
10202 CONTINUE
      IF(IERR/=0) THEN
        CALL WRIT40(' GXCHKI - CHEMKIN LINK FILE IS ABSENT   ')
        CALL SET_ERR(255,
     *        'Error. GXCHKI - CHEMKIN LINK FILE IS ABSENT.',1)
C        CALL EAROUT(1)
      ELSEIF(.NOT.NULLPR) THEN
        CALL WRITBL
        WRITE(LUPR1,'(''  GXCHKI - CHEMKIN LINK FILE: ''/10X,A70)')
     1                                                 NTEMP
        WRITE(LUPR3,'(''  GXCHKI - CHEMKIN LINK FILE: ''/10X,A70)')
     1                                                 NTEMP
      ENDIF
      LINK=LUNIT(55)
      LOUTCM=LUPR1
      CALL CKLEN(LINK,LOUTCM,LCKI,LCKR,LCKC)
      LCKERR=.FALSE.
C.... Allocate CHEMKIN arrays
      CALL CHEM_MEM(1,LCKR,LCKI,LCKC,1,NDCMR,NDCMI,NDCMC,NDCML)
      IF(LCKR>=NDCMR.OR.LCKI>=NDCMI.OR.LCKC>=NDCMC) THEN
        LCKERR=.TRUE.
        CALL WRIT40(' GXCHKI - Not enough room for CHEMKIN   ')
        CALL WRYT40(' GXCHKI - Not enough room for CHEMKIN   ')
        CALL WRIT40('   Currently:                           ')
        CALL WRIT3I('   NDCMR',NDCMR,'   NDCMI',NDCMI,
     1              '   NDCMC',NDCMC)
        CALL WRIT40('   Require:                             ')
        CALL WRIT3I('   NDCKR', LCKR,'   NDCKI', LCKI,
     1              '   NDCKC', LCKC)
        CALL SET_ERR(256,
     *        'Error. GXCHKI - Not enough room for CHEMKIN.',1)
C        CALL EAROUT(1)
      ENDIF
      CALL CKINIT(LCKI,LCKR,LCKC,LINK,LOUTCM,ICKMC,RCKMC,CCKMC1)
      CALL PHINDX(1,LINK,LOUTCM,PRECIS,0,0)
      CALL CKINDX(ICKMC,RCKMC,MM,KK,II,NFIT)
      CALL CLOSFL(55)
      IF(JTRY==0) THEN
        NMFIL(55)=NKEEP
      ELSE
        I55=IPEXT(55,JTRY/100,0)
        NMFLEX(I55)=NKEEP
      ENDIF
      IF(.NOT.NULLPR) THEN
        CALL WRITBL
        CALL WRIT40(' GXCHKI - CHEMKIN initialisation done   ')
        CALL WRITBL
        CALL WRIT40(PRECIS)
      ENDIF
C
C--- If MODified DIFfusion law has been selected, set UDIFF, UDIFNE
C    and USOURC. BUT the modifications, including the thermal diffusion
C    modifications are valid only for truly laminar flow, so we must
C    deactivate the modifications if the turbulent diffusivity for the
C    species are appreciable
C
      IF(GRNDNO(9,ENULA)) THEN
        IMCOPT=1
      ELSE
        IMCOPT=0
      ENDIF
      MODDIF=IMCOPT==1
      ENTFLX=MODDIF.AND.SOLTEM
      IF(NEZ(ENUT).AND.STOCMP) THEN
        DO 10009 K=1,KK
          LK=LBCKF+K-1
          MODDIF=MODDIF.AND.QGE(PRT(LK),1.E10)
10009   CONTINUE
      ENDIF
      IF(MODDIF.OR.ENTFLX) CALL SUB3L(UDIFF,.TRUE.,UDIFNE,.TRUE.,
     1                                USOURC,.TRUE.)
      IF(ICKOPT>0) USOLVE=.TRUE.
      ICALL1=MAX(0,ICALL1)
      IF(ICALL1==0) THEN
        ICALL1=2
      ELSE
        ICALL1=MIN(2,ICALL1)
      ENDIF
C
C--- Initialise TRANLIB if transport properties are required, ie.
C    look for activation of diffusion, solution and the specific
C    diffusivity option for CHEMKIN variables
C
      DODIFF=.FALSE.
      DO 10102 I=9,10
        DODIFF=DODIFF.OR.DFN(ITIND).AND.GRNDNO(I,-PRNDTL(ITIND))
        DO 10103 K=3,7,2
          DODIFF=DODIFF.OR.SOLVE(K).AND.GRNDNO(I,ENUL)
10103   CONTINUE
        IF(STOCMP) THEN
          DO 10101 K=1,KK
            LK=LBCKF+K-1
            DODIFF=DODIFF.OR.SOLVE(LK).AND.DFN(LK).AND.
     1                                 GRNDNO(I,-PRNDTL(LK))
10101     CONTINUE
        ENDIF
10102 CONTINUE
      MODDIF=MODDIF.AND.DODIFF
      DOTDIF=DOTDIF.AND.DODIFF
      IF(DODIFF) THEN
        CALL SUB2(KMCR,LCKR+1,KMCI,LCKI+1)
        NTEMP=' '
        DO 10211 ITRY=0,100,100
          IF(CMFILE/=' ') THEN
            I=INDEX(CMFILE,' ')
            IF(I==0) THEN
              I=4
            ELSE
              I=I-1
            ENDIF
            JTRY=ITRY
            IF(JTRY==0) THEN
              NKEEP=NMFIL(56)
              NMFIL(56)=CMFILE(1:I)//'mcln'
              CALL BLDNAM(NMFIL(56),NTEMP,LNMTMP)
            ELSE
              I56=IPEXT(56,JTRY/100,0)
              NKEEP=NMFLEX(I56)
              IC=INDEX(NMFLEX(I56),'tplink')
              IF(IC<=0.OR.IC>40) THEN
                IERR=1
                GO TO 10212
              ENDIF
              NKEEP=NMFLEX(I56)
              NMFLEX(I56)(IC:)=CMFILE(1:I)//'mcln'
              CALL BLDNAM(NMFLEX(I56),NTEMP,LNMTMP)
            ENDIF
          ENDIF
          IERR=-1
C#### JCL 24.06.13 copy file to slaves
          IF(MIMD.AND.NPROC>1.AND.JTRY==0)THEN
            CVAL=NMFIL(56)
            CALL PAROPENFL(90,CVAL,IERR,EXISTS)
            IF(EXISTS)THEN
              CALL GETSYSID(ISYSID)
              DELIM='/'
              IF(ISYSID<=2) DELIM=CHAR(92)
              CALL COPYBINARYFILE(DELIM,CVAL,IERR,IOS,EXISTS)
              IF(IERR/=0)GOTO 10211
            ELSE
              GOTO 10211
            ENDIF        
          ENDIF
          CALL OPENZZ(56+ITRY,IERR)
          IF(IERR==0) GO TO 10212
10211   CONTINUE
10212   CONTINUE
        IF(IERR/=0) THEN
          CALL WRIT40(' GXCHKI - TRANLIB LINK FILE IS ABSENT   ')
        CALL SET_ERR(257,
     *        'Error. GXCHKI - TRANLIB LINK FILE IS ABSENT.',1)
C          CALL EAROUT(1)
        ELSEIF(.NOT.NULLPR) THEN
          CALL WRITBL
          WRITE(LUPR1,'(''  GXCHKI - TRANLIB LINK FILE: ''/10X,A70)')
     1                                                 NTEMP
          WRITE(LUPR3,'(''  GXCHKI - TRANLIB LINK FILE: ''/10X,A70)')
     1                                                 NTEMP
        ENDIF
        LINK=LUNIT(56)
        CALL MCLEN(LINK,LOUTCM,LMCI,LMCR)
        IF(LMCR>=NDCMR-LCKR.OR.LMCI>=NDCMI-LCKI) THEN
          NEEDR=3*(LCKR+LMCR)/2
          NEEDI=3*(LCKI+LMCI)/2
          CALL CHEM_MEM(2,NEEDR,NEEDI,NDCMC,NDCML,NDCMR,NDCMI,NDCMC,
     1                                                        NDCML)
        ENDIF
        IF(LMCR>=NDCMR-LCKR.OR.LMCI>=NDCMI-LCKI) THEN
          LCKERR=.TRUE.
          CALL WRIT40(' GXCHKI - Not enough room for TRANLIB   ')
          CALL WRYT40(' GXCHKI - Not enough room for TRANLIB   ')
          CALL WRIT40('   Currently:                           ')
          CALL WRIT3I('   NDCMR',NDCMR,'   NDCMI',NDCMI,
     1              '   NDCMC',NDCMC)
          CALL WRIT40('   Require:                             ')
          CALL WRIT2I('   NDCKR', LCKR+LMCR,'   NDCKI', LCKI+LMCI)
        CALL SET_ERR(258,
     *        'Error. GXCHKI -  Not enough room for TRANLIB.',1)
C          CALL EAROUT(1)
        ENDIF
        CALL MCINIT(LINK,LOUTCM,LMCI,LMCR,ICKMC(KMCI),RCKMC(KMCR))
        CALL PHINDX(2,LINK,LOUTCM,PRECIS,NTDIF,K0TD)
        K0TD=K0TD+LCKI
        IF(.NOT.NULLPR) THEN
          CALL WRITBL
          CALL WRIT40(' GXCHKI - TRANLIB initialisation done   ')
        ENDIF
        CALL CLOSFL(56)
        IF(JTRY==0) THEN
          NMFIL(56)=NKEEP
        ELSE
          I56=IPEXT(56,JTRY/100,0)
          NMFLEX(I56)=NKEEP
        ENDIF
      ELSE
        CALL SUB3(IMCOPT,0,LMCR,0,LMCI,0)
      ENDIF
C
C--- Need an interface buffer and atomic weights
C
      IF(IMCOPT<2) THEN
C#### JCL/MRM 25.10.10 This should be number of reactions, II (US/022/10/10)      
C#### JCL 03.11.10 further correction - may need to be 2*no of species
C####        NBUF=KK*2
        NBUF=MAX(KK*2,II)
      ELSE
        NBUF=KK*(1+KK)
      ENDIF
      KK1=KK+1
      LENTWP=7*KK1+7*KK1+2
      K0MWT =LCKR+LMCR
C Introduce atomic-weight storage
C####      K0AWT= K0MWT+MM
C#### MRM 09.05.08 Correct error in storage allocation for K0MWT
      K0AWT= K0MWT+KK1
      K0VAR =K0AWT+KK1
C
      K0BUF =K0VAR+KK1
C#### JCL/MRM 25.10.10 seem to need extra KK1 here      
      K0H0K =K0BUF+NBUF+KK1
      NRTOT =K0H0K+KK1
      IF(ICKOPT==1) THEN
        K0ABV =NRTOT
        K0BLO =K0ABV +KK1
        K0TPB =K0BLO +KK1
        K0TWP =K0TPB +KK1
        K0SCR =K0TWP +LENTWP
        K0VARN=K0SCR +8*MAX(II,KK1)
        K0RES =K0VARN+KK1
        K0RESN=K0RES +KK1
        K0A   =K0RESN+KK1
        K0SU  =K0A   +KK1*KK1
        K0AP  =K0SU  +KK1
        NRTOT =K0AP  +KK1
      ENDIF
      IF(NRTOT>NDCMR) THEN
        NEEDR=NRTOT+1
        CALL CHEM_MEM(2,NEEDR,NDCMI,NDCMC,NDCML,NDCMR,NDCMI,NDCMC,NDCML)
      ENDIF
      IF(NRTOT>NDCMR) THEN
        CALL WRYT40(' GXCHKI - RCKMC too small for buffers   ')
        CALL WRIT40(' GXCHKI - RCKMC too small for buffers   ')
        CALL WRIT1I(' Add no.',NRTOT-NDCMR)
        CALL SET_ERR(259,
     *        'Error. GXCHKI -  RCKMC too small for buffers.',1)
C        CALL EAROUT(1)
      ENDIF
C
C--- Storage for indices to species, and elemental composition
C
      K0L0  =LCKI +LMCI
      K0ELT =K0L0 +KK
      K0IPV =K0ELT+MM*KK
      K0WDT =K0IPV+KK
      K0ROR =K0WDT+2*KK
      K0EBUS=K0ROR+2*II
      NITOT =K0EBUS+KK*II
      IF(NITOT>NDCMI) THEN
        NEEDI=NITOT+1
        CALL CHEM_MEM(2,NDCMR,NEEDI,NDCMC,NDCML,NDCMR,NDCMI,NDCMC,NDCML)
      ENDIF
      IF(NITOT>NDCMI) THEN
        CALL WRYT40(' GXCHKI - ICKMC too small for buffers   ')
        CALL WRIT40(' GXCHKI - ICKMC too small for buffers   ')
        CALL WRIT1I(' Add no.',NITOT-NDCMI)
        CALL SET_ERR(260,
     *        'Error. GXCHKI -  ICKMC too small for buffers.',1)
C        CALL EAROUT(1)
      ENDIF
C
C--- Storage for the chemical species and element names
C
      K0SYM =LCKC
      K0SYME=K0SYM+KK
c      K0AWT=K0SYME+MM
c      NCTOT=K0AWT+MM
      NCTOT=K0SYME+MM
      IF(NCTOT>NDCMC) THEN
        NEEDC=NCTOT+1
        CALL CHEM_MEM(2,NDCMR,NDCMI,NEEDC,NDCML,NDCMR,NDCMI,NDCMC,NDCML)
      ENDIF
      IF(NCTOT>NDCMC) THEN
        CALL WRYT40(' GXCHKI - CCKMC too small for buffers   ')
        CALL WRIT40(' GXCHKI - CCKMC too small for buffers   ')
        CALL WRIT1I(' Add no.',NCTOT-NDCMC)
        CALL SET_ERR(261,
     *        'Error. GXCHKI -  CCKMC too small for buffers.',1)
C        CALL EAROUT(1)
      ENDIF
C
C--- Storage for the ACTIVE flags for TWOPNT
C
      IF(ICKOPT==1) THEN
        NLTOT=KK
        IF(NLTOT>NDCML) THEN
          NEEDL=NLTOT+1
          CALL CHEM_MEM(2,NDCMR,NDCMI,NDCMC,NEEDL,NDCMR,NDCMI,NDCMC,
     1                                                        NDCML)
        ENDIF
        IF(NLTOT>NDCML) THEN
          CALL WRYT40(' GXCHKI - LCKMC too small for buffers   ')
          CALL WRIT40(' GXCHKI - LCKMC too small for buffers   ')
          CALL WRIT1I(' Add no.',NLTOT-NDCML)
        CALL SET_ERR(262,
     *        'Error. GXCHKI -  LCKMC too small for buffers.',1)
C          CALL EAROUT(1)
        ELSE
          DO 10005 K=1,KK
            LCKMC(K)=.FALSE.
10005     CONTINUE
        ENDIF
      ENDIF
C
C--- Set value of PRESS0 from CHEMKIN data; get no. of elements etc.;
C    get molecular weights; species names. If SIUNIT, then convert
C    pressure and RU to SI units.
C
      IF(EQZ(PREF)) PREF=1.0
      CALL CKRP(ICKMC,RCKMC,RU,RUC,PRES)
      IF(SIUNIT) THEN
        RU=RU*1.E-4
        PRES=PRES*1.E-1
      ENDIF
      PRESS0=PRES*PREF
      CALL CKWT(ICKMC,RCKMC,RCKMC(K0MWT+1))
C atomic-weight storage
      CALL CKAWT(ICKMC,RCKMC,RCKMC(K0AWT+1))
      CALL CKSYMS(CCKMC1,LCMOUT,CCKMC2(K0SYM+1),LCKERR)
      CALL CKSYME(CCKMC1,LCMOUT,CCKMC2(K0SYME+1),LCKERR)
      VARTEM=.NOT.CONTEM
      RTMP=273.0
      CALL CKHMS(RTMP,ICKMC,RCKMC,RCKMC(K0H0K+1))
      CALL CKNCF(MM,ICKMC,RCKMC,ICKMC(K0ELT+1))
      IF(LCKERR) THEN
        CALL WRYT40(' GXCHKI - Failed to load species names  ')
        CALL WRIT40(' GXCHKI - Failed to load species names  ')
        CALL SET_ERR(263,
     *        'Error. GXCHKI -  Failed to load species names.',1)
C        CALL EAROUT(1)
      ENDIF
C
C--- Check that ITEM comes after species
C
      IF(STORE(ITIND).AND.ITIND species  ')
        CALL WRIT40('   indices; please reorder the variables')
        CALL SET_ERR(264,
     *        'Error. GXCHKI - TEM1 index must be > species indices.',1)
C        CALL EAROUT(1)
      ELSE
        LCKERR=.FALSE.
        IF(STOCMP) THEN
          K1=LBCKF
          DO 10014 K=1,KK-1
            LCKERR=LCKERR.OR..NOT.STORE(K1)
            K1=K1+1
10014     CONTINUE
          LCKERR=LCKERR.OR..NOT.STORE(LBCKF+KK-1)
          IF(LCKERR) THEN
            CALL WRIT40(' GXCHKI - too few variables for species ')
            CALL WRIT40('   increase the no. of variables        ')
            CALL SET_ERR(265,
     *        'Error. GXCHKI - too few variables for species.',1)
C            CALL EAROUT(1)
          ENDIF
          IF(MOD(ISLN(LBCKF+KK-1),3)==0) THEN
            ISLN(LBCKF+KK-1)=ISLN(LBCKF+KK-1)/3
            F(L0ISL+LBCKF+KK-1)=-F(L0ISL+LBCKF+KK-1)
          ENDIF
        ENDIF
      ENDIF
C
C--- Declare storage for properties
C
C--- Settings to ensure PBP flags are set so that Ap is set,
C    and to ensure whole-field solver storage is deactivated
C
      IF(ICKOPT==1) THEN
C#### SCM 01.10.10 Update wholefield count if we reset ISLN
        IF(STOCMP) THEN
          K1=LBCKF
          DO 10021 K=1,KK
            IF(MOD(ISLN(K1),7)/=0) ISLN(K1)=ISLN(K1)*7
            IF(MOD(ISLN(K1),5)==0) THEN 
              ISLN(K1)=ISLN(K1)/5
              NWHOLE=NWHOLE-1
              NWHOLE=MAX0(1,NWHOLE)
            ENDIF
            K1=K1+1
10021     CONTINUE
        ENDIF
        IF(SOLTEM) THEN
          IF(MOD(ISLN(ITIND),7)/=0) ISLN(ITIND)=ISLN(ITIND)*7
          IF(MOD(ISLN(ITIND),5)==0) THEN
            ISLN(ITIND)=ISLN(ITIND)/5
            NWHOLE=NWHOLE-1
            NWHOLE=MAX0(1,NWHOLE)
          ENDIF
        ENDIF
      ENDIF
      NFTOT0=NFTOT
      CALL GXMAKE(L0DTF,KK1,'DTF ')
      NXYZ=NXNY
      NXYZ2=NXNY
      IF(WHF(LBCKF)) NXYZ=NXYZ*NZ
      IF(NZ>1) NXYZ2=2*NXNY
      CALL GXMAKE(L0HK  ,NXYZ2*KK,'HK  ')
      CALL GXMAKE(L0H   ,NXYZ,'H   ')
      CALL GXMAKE(L0CP  ,NXYZ,'CP  ')
      IF(DODIFF) THEN
        IF(IMCOPT<2) CALL GXMAKE(L0DK ,MAX(NXYZ,NXYZ2)*KK,'DK  ')
        IF(IMCOPT>=1) THEN
          IF(IMCOPT>=2) THEN
            CALL GXMAKE(L0DJK ,NXNY*KK*KK,'DJK ')
            IF(NZ>1) CALL GXMAKE(L0DJKH,NXNY*KK*KK,'DJKH')
          ENDIF
          IF(NX>1) CALL GXMAKE(L0VDE,NXYZ*KK,'VDFE')
          IF(NY>1) CALL GXMAKE(L0VDN,NXYZ*KK,'VDFN')
          IF(NZ>1) CALL GXMAKE(L0VDH,NXYZ*KK,'VDFH')
        ENDIF
C
C--- Determine whether or not we require the thermal diffusivities
C    on the basis of are MODifying the DIFfusion law, and whether
C    any of the species undergo thermal diffusion
C
        DOTDIF=GRNDNO(9,ENULB).AND.NTDIF>0.AND.MODDIF
        IF(DOTDIF) CALL GXMAKE(L0DTK,MAX(NXYZ,NXYZ2)*KK,'DTK ')
        IF(MODDIF.OR.ENTFLX) THEN
          CALL GXMAKE(L0HSU,MAX(NXYZ,NXYZ2),'HSU ')
          CALL GXMAKE(L0HAP,MAX(NXYZ,NXYZ2),'HAP ')
          IF(NZ>1) THEN
            CALL GXMAKE(L0SUHI,NXNY,'SUHI')
            CALL GXMAKE(L0APHI,NXNY,'APHI')
          ENDIF
        ENDIF
      ENDIF
      IF(DOSORC) THEN
        CALL GXMAKE(L0CDK,NXYZ*KK,'CDOT')
        CALL GXMAKE(L0TDK,NXYZ*KK,'TAUD')
        IF(ICKOPT==0) THEN
          CALL GXMAKE(L0HCO ,NXYZ,'HCO ')
          CALL GXMAKE(L0HVAL,NXYZ,'HVAL')
        ELSEIF(ICKOPT==1) THEN
          L0SUK=L0CDK
          L0APK=L0TDK
cccc          USOLVE=.TRUE.
          USOLVE=.TRUE.
        ENDIF
      ENDIF
C      call writ2i('l0suk   ',l0suk,'l0apk   ',l0apk)
C
C--- Look for Inlet Diffusion PATCHes and create PATCH-wise
C    variables
C
      IF(SOLTEM) THEN
        DO 10022 IP=1,NUMPAT
          IF(NAMPAT(IP)(1:6)=='CHEMID') THEN
            CALL MAKEPV(PVCKID,IP)
          ENDIF
10022   CONTINUE
      ENDIF
C
C--- Print memory requirements into the RESULT file
C
      IF(.NOT.NULLPR) THEN
        CALL WRITBL
        CALL WRIT40(' GXCHKI - The following memory has been ')
        CALL WRIT40('          taken for CHEMKIN & TRANLIB   ')
        CALL WRITBL
        CALL WRIT2I('REAL:NO.',NRTOT,'INT.:NO.',NITOT)
        CALL WRIT2I('CHR.:NO.',NCTOT,'In F:NO.',NFTOT-NFTOT0)
        IF(ICKOPT==1) CALL WRIT1I('LOG.:NO.',NLTOT)
      ENDIF
C
C--- Set up indices for storing WDoTs and Rates-Of-Reaction
C
      IF(STOCMP) THEN
        IGRND=NINT(GRND)
        ICKMC(K0WDT+1)=IGRND
        DO 10031 K=2,2*KK
          ICKMC(K0WDT+K)=0
10031   CONTINUE
        ICKMC(K0ROR+1)=IGRND
        DO 10032 K=2,2*II
          ICKMC(K0ROR+K)=0
10032   CONTINUE
        DO 10025 K=16,NPHI
          NAMD=' '
          K0=0
          IF(NAME(K)(4:4)=='+') THEN
            NAMD=NAME(K)(1:3)//' '
            I=LBNAME(NAMD)
            K0=K0WDT
          ELSEIF(NAME(K)(3:3)=='+') THEN
            NAMD=NAME(K)(1:2)//'  '
            I=LBNAME(NAMD)
            K0=K0WDT
          ELSEIF(NAME(K)(2:2)=='+') THEN
            NAMD=NAME(K)(1:1)//'   '
            I=LBNAME(NAMD)
            K0=K0WDT
          ELSEIF(NAME(K)(4:4)=='&') THEN
            READ(NAME(K),'(I3)') I
            K0=K0ROR
          ELSEIF(NAME(K)(3:3)=='&') THEN
            READ(NAME(K),'(I2)') I
            K0=K0ROR
          ELSEIF(NAME(K)(2:2)=='&') THEN
            READ(NAME(K),'(I1)') I
            K0=K0ROR
          ENDIF
          IF(K0>0) THEN
            IF(K0==K0WDT) THEN
              IF(ILBCKF+KK-1) THEN
                CALL WRIT40(' GXCHKI - Failed to get WDT var. name   ')
                CALL SET_ERR(266,
     *            'Error. GXCHKI - Failed to get WDT var. name.',1)
C                CALL EAROUT(1)
              ELSE
                I=I-LBCKF+1
              ENDIF
            ELSEIF(K0==K0ROR) THEN
              IF(I<1.OR.I>II) THEN
                CALL WRIT40(' GXCHKI - Failed to get ROR var. name   ')
                CALL SET_ERR(266,
     *            'Error. GXCHKI - Failed to get ROR var. name.',1)
C                CALL EAROUT(1)
              ENDIF
            ENDIF
            IF(ICKMC(K0+1)==IGRND) ICKMC(K0+1)=0
            ICKMC(K0+I)=K
          ENDIF
10025   CONTINUE
      ENDIF
C
      IF(.NOT.NULLPR) THEN
        CALL WRITBL
        CALL WRIT40(' The chemical species treated are:      ')
        WRITE(LUPR1, '(''  No.   Name  '',
     1               ''                   Mol.Wt.  (Status)'')')
        IF(STOCMP) THEN
          CODE='(Solved) '
        ELSE
          CODE='(Virtual)'
        ENDIF
      ENDIF
      DO 10012 K=1,KK
        IF(STOCMP) THEN
          IF(.NOT.SOLVE(LBCKF+K-1)) CODE='(Derived)'
        ENDIF
        IL=INDEX(CCKMC2(K0SYM+K),' ')
        IF(IL>0.AND.IL<=5) THEN
          IF(STOCMP) NAME(LBCKF+K-1)=CCKMC2(K0SYM+K)
          IF(.NOT.NULLPR) WRITE(LUPR1,'(2X,I3,4X,A16,8X,F6.2,3X,A9)')
     1         K,CCKMC2(K0SYM+K),RCKMC(K0MWT+K),CODE
        ELSE
          IF(.NOT.NULLPR.AND.STOCMP) WRITE(LUPR1,'(2X,I3,4X,A16,
     1     '' ('',A4,'') '',F6.2,3X,A9)') K,CCKMC2(K0SYM+K),
     2               NAME(LBCKF+K-1),RCKMC(K0MWT+K),CODE
        ENDIF
10012 CONTINUE
      IF(.NOT.NULLPR) THEN
        CALL WRITBL
        CALL WRIT1R(' PRESS0 ',PRESS0)
        CALL WRITBL
        IF(CONTEM) THEN
          CALL WRIT40(' CONTEM=T  constant temperature option  ')
          CALL WRITBL
          CALL WRIT1R(' TEMP0  ',TEMP0)
          CALL WRITBL
        ENDIF
        IF(.NOT.SIUNIT) THEN
          CALL WRIT40(' CGS units: cm, g, erg, mol, dynes, s, K')
          CALL WRITST
          CALL WRITBL
        ENDIF
      ENDIF
C
C--- Copy the DTFALS values for the CHEMKIN species
C
      IF(STOCMP) THEN
        K1=LBCKF
        DO 10016 K=1,KK
          F(L0DTF+K)=DTFALS(K1)
          K1=K1+1
10016   CONTINUE
      ENDIF
      IF(SOLTEM) F(L0DTF+KK+1)=DTFALS(ITIND)
C
C--- Special initialisations for the P-B-P Newton-Raphson solution
C
      IF(ICKOPT==1) THEN
C
C--- Settings for the TWOPNT driver
C
        IF(NJAC==0)   NJAC  =20
        IF(ITJAC==0)  ITJAC =0
        IF(IRETIR==0) IRETIR=50
        IF(NUMDT==0)  NUMDT =50
        NUMDT=MAX(0,NUMDT)
        IF(NINIT<0)  NINIT =50
        NINIT=MAX(0,NINIT)
        UFAC=UFAC1
        DFAC=DFAC1
        IF(EQZ(UFAC1)) UFAC=2.
        IF(EQZ(DFAC1)) DFAC=4.1
        UFAC=MAX(REAL(UFAC),1.01)
        DFAC=MAX(REAL(DFAC),1.01)
        DTMIN=DTMIN1
        DTMAX=DTMAX1
        IF(EQZ(DTMIN1)) DTMIN=1.E-10
        IF(EQZ(DTMAX1)) DTMAX=0.1
C
C--- Determine U, the smallest increment to 1.0 which the computer
C    can detect, ie. the computer's resolution
C
        U = 1.0
10017   CONTINUE
        U  = U*0.5
        U0 = 1.0 + U
c!!!! do NOT replace the next line by IF(QNE ....)) !!!!
        IF(U0 /= 1.0) GO TO 10017
        ABSOL = SQRT(2.0*U)
        RELAT = ABSOL
        IF(SOLTEM) THEN
          NY0=1
          RCKMC(K0ABV+1)=VARMAX(ITIND)
          RCKMC(K0BLO+1)=VARMIN(ITIND)
        ELSE
          NY0=0
        ENDIF
        K1=LBCKF
        ATOL=ENDIT(LBCKF)
        IF(STOCMP) THEN
          DO 10018 K=1,KK-1
            ATOL=MIN(REAL(ATOL),ENDIT(K1))
            RCKMC(K0ABV+NY0+K)=VARMAX(K1)
            RCKMC(K0BLO+NY0+K)=VARMIN(K1)
            K1=K1+1
10018     CONTINUE
        ENDIF
C
        IF(.NOT.STOCMP.AND.IMIXF>0)
     1        CALL CKNU(KK,ICKMC,RCKMC,ICKMC(K0EBUS+1))
C
      ENDIF
C
      GO TO 999
C-----------------------------------------------------------------
C   * Make changes for this group only in group 19.
C--- GROUP 7. Variables stored, solved & named
C-----------------------------------------------------------------
C
C--- GROUP 8. Terms (in differential equations) & devices
C
    8 GO TO (81,82,83,84,85,86,87,88,89,810,811,812,813,814,815)
     1,ISC
   81 CONTINUE
C   * ------------------- SECTION  1 ---------------------------
C    For U1AD<=GRND--- phase 1 additional velocity. Index VELAD
      GO TO 999
   82 CONTINUE
C   * ------------------- SECTION  2 ---------------------------
C    For U2AD<=GRND--- phase 2 additional velocity. Index VELAD
      GO TO 999
   83 CONTINUE
C   * ------------------- SECTION  3 ---------------------------
C    For V1AD<=GRND--- phase 1 additional velocity. Index VELAD
      GO TO 999
   84 CONTINUE
C   * ------------------- SECTION  4 ---------------------------
C    For V2AD<=GRND--- phase 2 additional velocity. Index VELAD
      GO TO 999
   85 CONTINUE
C   * ------------------- SECTION  5 ---------------------------
C    For W1AD<=GRND--- phase 1 additional velocity. Index VELAD
      GO TO 999
   86 CONTINUE
C   * ------------------- SECTION  6 ---------------------------
C    For W2AD<=GRND--- phase 2 additional velocity. Index VELAD
      GO TO 999
   87 CONTINUE
C   * ------------------- SECTION 7 ---- Volumetric source for gala
      GO TO 999
   88 CONTINUE
C   * ------------------- SECTION 8 ---- Convection fluxes
      GO TO 999
   89 CONTINUE
C   * ------------------- SECTION 9 ---- Diffusion coefficients
C--- Entered when UDIFF =.TRUE.; block-location indices are LAE
C    for east, LAW for west, LAN for north, LAS for
C    south, LD11 for high, and LD11 for low.
C    User should provide INDVAR and NDIREC IF's as above.
C    EARTH will apply the DIFCUT and GP12 modifications after the user
C    has made his settings.
C
C--- Modify the diffusion coefficients to give derivatives of mole-
C    fractions. See the section for diffusion neighbours.
C
      IF(MODDIF.AND.INDVAR>=LBCKF.AND.
     1   INDVAR<=LBCKF+KK-2) THEN
        IF(NDIREC==1) THEN
          CALL SUB3(KDP,L0F(LAN),KDM,L0F(LAS),L0VDI,L0VDN)
          CALL SUB2(INI,INN,L0DX,L0F(LXYDY))
        ELSEIF(NDIREC==3) THEN
          CALL SUB3(KDP,L0F(LAE),KDM,L0F(LAW),L0VDI,L0VDE)
          CALL SUB2(INI,INE,L0DX,L0F(LXYDX))
        ELSEIF(NDIREC==5) THEN
          CALL SUB3(KDP,L0F(LD11),KDM,0,L0VDI,L0VDH)
          CALL SUB3(INI,INH,L0DX,L0F(LXYDZ),L0DXH,L0F(LXYDZG))
        ENDIF
        PRES=PRESS0
        RTMP=TEMP0
        TMPNEI=TEMP0
        PRSNEI=PRESS0
        K=INDVAR-LBCKF+1
        K0IV=ICKMC(K0L0+K)
        IVOFF =IZOFF*KK+(INDVAR-LBCKF)*NXNY
        IVTOFF=IZOFF*kk+(INDVAR-LBCKF)*NXNY
        K0IV=ICKMC(K0L0+K)
        DO 80911 ICELL=1,NXNY
          IF(VARTEM) RTMP=TEMP0+F(L0TEM+ICELL)
          IF(VARPRS) PRES=PRESS0+F(L0P1+ICELL)
          WTMEAN=0.0
          DO 80905 KA=1,KK
            WTMEAN=WTMEAN+F(ICKMC(K0L0+KA)+ICELL)/RCKMC(K0MWT+KA)
80905     CONTINUE
          WTMEAN=1./WTMEAN
c          write(lupr1,'('' species mass fractions:''/1p5e12.5)')
c     1           (f(ickmc(k0l0+ka)+icell),ka=1,kk)
c          call writ1r(' wtmean ',wtmean)
          FLUXTD=F(KDP+ICELL)
          IF(KDP>0) F(KDP+ICELL)=FLUXTD*WTMEAN
          IF(KDM>0) F(KDM+ICELL)=F(KDM+ICELL)*WTMEAN
          IF(DOTDIF) THEN
C
C--- For thermal diffusion, we must first calculate an equivalent
C    mean diffusion coefficient
C
            INEI=NINT(F(INI+ICELL))
            IF(NDIREC==5) THEN
              INEI1=ICELL+KK*NXNY
              INEI2=ICELL+L0DXH-L0DX
              DXNEI=2.*F(L0DXH+ICELL)-F(L0DX+ICELL)
              INEI3=ITWO(ICELL,INEI,.NOT.STORE(DEN1).AND.RHO1>=0.)
            ELSE
              INEI1=INEI
              DXNEI =F(L0DX+INEI)
              INEI3=INEI
            ENDIF
            RHO   =F(L0DEN+ICELL)
            RHONEI=F(L0DEN+INEI3)
            IF(VARPRS) PRSNEI=PRESS0+F(L0P1+INEI)
            IF(VARTEM) TMPNEI=TEMP0+F(L0TEM+INEI)
            DTK   =F(L0DTK+IVTOFF+ICELL)+TINY
            DTKNEI=F(L0DTK+IVTOFF+INEI1)+TINY
            DCK   =F(L0DK+IVOFF+ICELL)+TINY
            DCKNEI=F(L0DK+IVOFF+INEI1)+TINY
            DX    =F(L0DX+ICELL)
            IF(HRM(INDVAR)) THEN
              AVEDTK=(DX/(RHO*DCK)+DXNEI/(RHONEI*DCKNEI))/
     1               (DX/(RHO*DTK*DCK)+DXNEI/(RHONEI*DTKNEI*DCKNEI))
            ELSE
              AVEDTK=(RHO*DTK*DCK+RHONEI*DTKNEI*DCKNEI)/
     1               (RHO*DCK+RHONEI*DCKNEI)
            ENDIF
            GRAT=TMPNEI/RTMP
            FLUXTD=-AVEDTK*FLUXTD*WTMEAN*LOG(TMPNEI/RTMP)
            F(L0VDI+IVTOFF+ICELL)=FLUXTD
          ENDIF
80911   CONTINUE
C
      ENDIF
      GO TO 999
  810 CONTINUE
C   * ------------------- SECTION 10 --- Convection neighbours
      GO TO 999
  811 CONTINUE
C   * ------------------- SECTION 11 --- Diffusion neighbours
C
C--- Modify the diffusion neighbours to take account of mole-
C    fraction derivatives and thermal diffusion
C
      IF((ENTFLX.OR.MODDIF).AND.INDVAR>=LBCKF.AND.
     1   INDVAR<=LBCKF+KK-2) THEN
        IF(WHF(LBCKF)) THEN
          CALL SUB2(KSU,-LBZ(L3SU,IZ),KAP,-LBZ(L3AP,IZ))
        ELSE
          CALL SUB2(KSU,L0F(LSU),KAP,L0F(LAP))
        ENDIF
        K=INDVAR-LBCKF+1
        IF(INDVAR==LBCKF.AND.SOLTEM) THEN
          CALL FN1(-L0HAP,0.0)
          CALL FN1(-L0HSU,0.0)
          IF(NZSTEP>1) THEN
            CALL FN1(-L0HAP-NXNY,0.0)
            CALL FN1(-L0HSU-NXNY,0.0)
          ENDIF
        ENDIF
        IF(NDIREC==1) THEN
          CALL SUB4(KDP,L0F(LAN),KDM,L0F(LAS),L0VDI,L0VDN,
     1              KDNEI,L0F(LD8))
          CALL SUB2(INI,INN,L0DX,L0F(LXYDY))
        ELSEIF(NDIREC==3) THEN
          CALL SUB4(KDP,L0F(LAE),KDM,L0F(LAW),L0VDI,L0VDE,
     1              KDNEI,L0F(LD8))
          CALL SUB2(INI,INE,L0DX,L0F(LXYDX))
        ELSEIF(NDIREC==5) THEN
C Must use LAH here not LD11
          CALL SUB4(KDP,L0F(LAH),KDM,0,L0VDI,L0VDH,KDNEI,L0F(LD7))
          CALL SUB4(INI,INH,L0DX,L0F(LXYDZ),K0IHH,
     1      L0HK+(K-1)*NXNY+KK*NXNY,K0HHKK,L0HK+(KK-1)*NXNY+KK*NXNY)
          CALL FN1(-L0APHI,0.0)
          CALL FN1(-L0SUHI,0.0)
        ENDIF
        PRES=PRESS0
        PRESP=PRESS0
        RTMP=TEMP0
        RTMP=TEMP0
        CALL SUB4(K0IV,ICKMC(K0L0+K),K0IH,L0HK+(K-1)*NXNY,
     1            K0IHKK,L0HK+(KK-1)*NXNY,IVTOFF,IZOFF*KK+(K-1)*NXNY)
        DO 81111 ICELL=1,NXNY
          INEI=NINT(F(INI+ICELL))
          IF(MODDIF.OR.ENTFLX) THEN
            IF(VARPRS) THEN
              PRES =PRESS0+F(L0P1+ICELL)
              PRESP=PRESS0+F(L0P1+INEI)
            ENDIF
            IF(VARTEM) THEN
              RTMP  =TEMP0+F(L0TEM+ICELL)
              TMPNEI=TEMP0+F(L0TEM+INEI)
            ENDIF
          ENDIF
          IF(MODDIF) THEN
            WTMEAN=0.0
            WTMNEI=0.0
            DO 81105 KA=1,KK
              WTMEAN=WTMEAN+F(ICKMC(K0L0+KA)+ICELL)/RCKMC(K0MWT+KA)
              WTMNEI=WTMNEI+F(ICKMC(K0L0+KA)+INEI )/RCKMC(K0MWT+KA)
81105       CONTINUE
            WTMEAN=1./WTMEAN
            WTMNEI=1./WTMNEI
            WTMRAT=WTMEAN/WTMNEI
          ELSE
            WTMRAT=1.0
          ENDIF
          FLUX=0.0
          IF(DOTDIF) THEN
            FLUXTD=F(L0VDI+IVTOFF+ICELL)
            IF(NDIREC==5) THEN
              KAPNEI=L0APHI+ICELL
              KSUNEI=L0SUHI+ICELL
            ELSE
              KAPNEI=KAP+INEI
              KSUNEI=KSU+INEI
            ENDIF
            IF(FLUXTD<0.0) THEN
              F(KAP+ICELL)=F(KAP+ICELL)-FLUXTD/F(K0IV+ICELL)
            ELSEIF(FLUXTD>0.0) THEN
              F(KAPNEI) =F(KAPNEI)+FLUXTD/F(KDNEI+ICELL)
            ENDIF
            F(KSU+ICELL)=F(KSU+ICELL) - FLUXTD
            F(KSUNEI)   =F(KSUNEI)    + FLUXTD
            FLUX=FLUX+FLUXTD
          ENDIF
          IF(MODDIF) F(KDNEI+ICELL)=F(KDNEI+ICELL)/WTMRAT
          FLUX=FLUX+F(KDP+ICELL)*(F(K0IV+ICELL)-F(KDNEI+ICELL))
C--- Accumulate the enthalpy flux that results from mass-diffusion
C    (note that the arrays for SU and AP terms must be zeroed on
C    entry for the first CHEMKIN species)
C
          IF(ENTFLX) THEN
C
C--- The enthalpy fluxes must be upwinded
C
            ENTI  =F(K0IH+ICELL)  -F(K0IHKK+ICELL)
            IF(NDIREC==5) THEN
              ENTNEI=F(K0IHH+ICELL) -F(K0HHKK+ICELL)
              INEIH=NXNY+ICELL
            ELSE
              ENTNEI=F(K0IH+INEI)   -F(K0IHKK+INEI)
              INEIH=INEI
            ENDIF
            IF(FLUX<0.0) THEN
              F(L0HAP+ICELL)=F(L0HAP+ICELL)-FLUX*ENTI/RTMP
              F(L0HSU+ICELL)=F(L0HSU+ICELL)-FLUX*(ENTNEI-ENTI)
            ELSEIF(FLUX>0.0) THEN
              F(L0HAP+INEIH)=F(L0HAP+INEIH)+FLUX*ENTNEI/TMPNEI
              F(L0HSU+INEIH)=F(L0HSU+INEIH)+FLUX*(ENTNEI-ENTI)
            ENDIF
          ENDIF
81111   CONTINUE
      ENDIF
C
      GO TO 999
  812 CONTINUE
C   * ------------------- SECTION 12 --- Linearised sources
C
C--- Add the enthalpy flux resulting from mass-diffusion to the
C    SU for the energy conservation equation, and make the
C    corresponding additions to AP
C
      IF((ENTFLX.OR.MODDIF).AND.SOLTEM.AND.INDVAR==ITIND) THEN
        IF(WHF(LBCKF)) THEN
          CALL SUB2(KAP,LBZ(L3AP,IZ),KSU,LBZ(L3SU,IZ))
        ELSE
          CALL SUB2(KAP,-L0F(LAP),KSU,-L0F(LSU))
        ENDIF
        CALL FN60(KAP,-L0HAP)
        CALL FN60(KSU,-L0HSU)
      ENDIF
C
      GO TO 999
  813 CONTINUE
C   * ------------------- SECTION 13 --- Correction coefficients
      GO TO 999
  814 CONTINUE
C   * ------------------- SECTION 14 --- User's own solver
C
C--- The CHEMKIN variables are solved on a point-by-point (PBP)
C    basis using the Sandia TWOPNT solver (Newton-Raphson iteration)
C    If this is not a CHEMKIN variable (assumed to be the default)
C    then we want to go into the PHOENICS solver so we set USER=.F.
C    For CHEMKIN variables we must initialise LD7, into which the
C    PHOENICS solver would put the corrections, to zero.
C
        USER=.FALSE.
C
C--- Set up indices for storage of Rates-Of-Reaction and WDoTs
C
      IF(STOCMP) THEN
        IF(ICKMC(K0WDT+1)/=IGRND) THEN
          DO 81401 K=1,KK
            I=ICKMC(K0WDT+K)
            IF(I>0) ICKMC(K0WDT+KK+K)=L0F(I)
81401     CONTINUE
        ENDIF
        IF(ICKMC(K0ROR+1)/=IGRND) THEN
          DO 81402 K=1,II
            I=ICKMC(K0ROR+K)
            IF(I>0) ICKMC(K0ROR+II+K)=L0F(I)
81402     CONTINUE
        ENDIF
      ENDIF
C
      IF((SOLTEM.AND.INDVAR==ITIND.OR.
     1   .NOT.SOLTEM.AND.INDVAR==LBCKF+KK-2)
     2   .AND.ICKOPT==1) THEN
C
C--- Initialise ENUFSW for later setting on the basis of residual
C    on entry to PHOTWO, set USER to cut out the PHOENICS solver,
C    and ensure that the dummy array (LD7) that usually contains
C    the corrections is set to zero
C
        ENUFSW=.TRUE.
        DO 81411 K=1,KK
          LK=LBCKF+K-1
          IF(ISL(LK)/=0) F(L0SLRS+IABS(ISL(LK)))=0.0
81411   CONTINUE
        IF(SOLTEM) F(L0SLRS+ISL(ITIND))=0.0
        IF(IZSTEP==1) THEN
          DO 81412 K=1,KK
            LK=LBCKF+K-1
            IF(ISL(LK)/=0) F(L0TTRS+IABS(ISL(LK)))=0.0
81412     CONTINUE
          IF(SOLTEM) F(L0TTRS+ISL(ITIND))=0.0
        ENDIF
        USER=.TRUE.
        CALL FN1(LD7,0.0)
C
C--- Final CHEMKIN variable so now loop over cells solving on a point-
C    by-point basis: first add in the extra diffusion sources for
C    enthalpy and set up indices dependent on whether or not the
C    temperature variable is being solved
C
        CALL SUB3(L0SU,L0F(LSU),L0AP,L0F(LAP),L0VOL,L0F(VOL))
        IF(INDVAR==ITIND) THEN
          NATJ=KK
          NY0=1
        ELSE
          NATJ=KK-1
          NY0=0
        ENDIF
        I0=L0SUK
        ISA=L0APK-L0SUK
        PRES=PRESS0
        RTMP=TEMP0
        DTTP=DTFALS(LBCKF)/(UFAC-1.)
C
C--- Now loop over cells calling the solver-driver PHOTWO to generate
C    the solutions
C
        DO 81441 ICELL=1,NXNY
C
C--- Before doing anything else look to see whether the cell is
C    blocked. If blocked we do nothing for the chemical species
C    but still do a PBP solution for the TEM1 variable
C
          CALL SUB2L(BLOCKD,.FALSE.,SOLID,.FALSE.)
          IF(IPRPS>0) THEN
            SOLID =F(L0PRPS+ICELL)>=SOLPRP
            BLOCKD=F(L0PRPS+ICELL)>=PORPRP
          ENDIF
C         IF(TSTVPO) BLOCKD=BLOCKD.OR.F(L0VPOR+ICELL)<=1.E-10
          IF(SOLID.AND.INDVAR==ITIND) THEN
C
C--- If the cell is a solid then we solve only for temperature, and
C    if temperature is solved we skip over all solution
C
            RESID=ABS(F(L0SU+ICELL))
            DELVAR=F(L0SU+ICELL)/F(L0AP+ICELL)
            F(L0TEM+ICELL)=F(L0TEM+ICELL)+DELVAR
            GO TO 81427
          ELSEIF(BLOCKD.OR.SOLID) THEN
C
C--- If the cell is completely blocked, do nothing
C
            GO TO 81439
          ENDIF
C--- If FIXVAL cell then (a) fix temperature and species; and
C    (b) skip to end of loop bypassing call to photwo
            FIXCEL=.FALSE.
            IF(NY0>0) THEN
              IF(F(L0AP+ICELL)>=FIXVAL) THEN
                FIXCEL=.TRUE.
                DELVAR=F(L0SU+ICELL)/F(L0AP+ICELL)
                F(L0TEM+ICELL)=F(L0TEM+ICELL)+DELVAR
              ENDIF
            ENDIF
            DO 81499 K=1,NATJ-1
              IF(F(I0+ISA+K)>=FIXVAL) THEN
                FIXCEL=.TRUE.
                DELVAR=F(I0+K)/F(I0+ISA+K)
                F(ICKMC(K0L0+K)+ICELL)=F(ICKMC(K0L0+K)+ICELL)+DELVAR
              ENDIF
81499       CONTINUE
            IF(FIXCEL) GO TO 81439
C---
C
C--- First copy the sources, Su, the coefficients, Ap, and the CHEMKIN
C    variables solved into buffer stores, set the pressure, if local
C    pressure is to be used in CHEMKIN, and get the cell volume for
C    calculating the reaction source and internal false time-step
C    terms
C
            IF(VARPRS) PRES=PRESS0+F(L0P1+ICELL)
            IF(VARTEM) RTMP=F(L0TEM+ICELL)
            DO 81405 K=1,KK
              RCKMC(K0SU+NY0+K)  =F(I0+K)
              RCKMC(K0AP+NY0+K)  =F(I0+ISA+K)
              RCKMC(K0VAR+NY0+K) =F(ICKMC(K0L0+K)+ICELL)
              RCKMC(K0VARN+NY0+K)=F(ICKMC(K0L0+K)+ICELL)
81405       CONTINUE
            IF(INDVAR==ITIND) THEN
              RCKMC(K0VAR+1)=F(L0TEM+ICELL)
              RCKMC(K0VARN+1)=F(L0TEM+ICELL)
              RCKMC(K0SU+1)=F(L0SU+ICELL)
              RCKMC(K0AP+1)=F(L0AP+ICELL)
            ELSE
              RCKMC(K0SU+NATJ)=F(L0SU+ICELL)
              RCKMC(K0AP+NATJ)=F(L0AP+ICELL)
            ENDIF
            VOLUME=F(L0VOL+ICELL)
            if(dbgloc) call writ1i('   icell',icell)
C
C--- Call the PHOenics TWOpnt driver routine PHOTWO to generate the
C    solution for the current cell
C
            CALL PHOTWO(LOUTCM,KK,NATJ,1,1,SOLTEM,LENTWP,LEVEL,
     1                  LCKMC,ICKMC,RCKMC,VOLUME,RCKMC(K0MWT+1),
     2                  ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,RCKMC(K0ABV+1),
     3                  RCKMC(K0BLO+1),RCKMC(K0TPB+1),RCKMC(K0TWP+1),
     4                  ICKMC(K0IPV+1),RCKMC(K0SCR+1),RCKMC(K0VAR+1),
     5                  RCKMC(K0VARN+1),RCKMC(K0RES+1),RCKMC(K0RESN+1),
     6                  RCKMC(K0A+1),RCKMC(K0SU+1),RCKMC(K0AP+1),
     7                  RCKMC(K0H0K+1),UFAC,DFAC,PRES,RTMP,NINIT,NJAC,
     8                  ITJAC,IRETIR,NUMDT,DTTP,DTMIN,DTMAX,
     9                  QDOT,ICALL1,DTEMP0,DBGLOC.OR.FLAG,ERROR,TPSEM)
C
C--- First check for errors reported by PHOTWO
C
            IF(ERROR) THEN
              IF(BYPASS) THEN
                IF(TPSEM) THEN
                 CALL WRITBL
                 CALL WRITST
                 CALL WRIT40(' GXCHKI - PHOTWO failed to converge at  ')
                 CALL WRIT4I('  ISWEEP',ISWEEP,'   ITHYD',ITHYD,
     1                       '  IZSTEP',IZSTEP,'   ICELL',ICELL)
                ENDIF
              ELSE
                CALL WRITBL
                CALL WRITST
                CALL WRIT40(' GXCHKI - PHOTWO reports a fatal error  ')
                CALL WRIT1I('INDVAR  ',INDVAR)
                CALL WRIT4I('  ISWEEP',ISWEEP,'   ITHYD',ITHYD,
     1                      '  IZSTEP',IZSTEP,'   ICELL',ICELL)
                CALL WRIT40('     Solution for this cell is:         ')
                WRITE(LUPR1,'('' Spec.       Value       VARMIN'',
     1                        ''      VARMAX         Ap       Res'')')
                IF(NY0>0) WRITE(LUPR1,'(1X,A4,3X,1P5E12.3)') 'TEM1',
     1                RCKMC(K0VAR+1), VARMIN(ITIND), VARMAX(ITIND),
     2                RCKMC(K0AP+1), RCKMC(K0SCR+2*KK+1)
                DO 81491 K=NY0+1,NATJ+1
                   KIV=LBCKF+K-NY0-1
                   IF(SOLVE(KIV)) THEN
                     WRITE(LUPR1,'(1X,A4,3X,1P5E12.3)') NAME(KIV),
     1                      RCKMC(K0VAR+K), VARMIN(KIV), VARMAX(KIV),
     2                      RCKMC(K0AP+K), RCKMC(K0SCR+2*KK+K)
                   ELSE
                     WRITE(LUPR1,'(1X,A4,3X,1P5E12.3)') NAME(KIV),
     1                      RCKMC(K0VAR+K), VARMIN(KIV), VARMAX(KIV)
                   ENDIF
81491           CONTINUE
                CALL WRITBL
                CALL WRITST
                CALL WRITBL
                CALL SET_ERR(267, 'Error. See result file.',1)
C                CALL EAROUT(1)
              ENDIF
            ENDIF
C
C--- Copy the new solution back into the PHOENICS variable stores,
C    calculating the mass-fraction of the last species form the sum
C    of the other species mass-fractions
C
            SUMY=0.0
            DO 81431 K=1,KK-1
              SUMY=SUMY+RCKMC(K0VAR+NY0+K)
              RCKMC(K0BUF+K)=RCKMC(K0VAR+NY0+K)
              F(ICKMC(K0L0+K)+ICELL)=RCKMC(K0VAR+NY0+K)
81431       CONTINUE
            RCKMC(K0BUF+KK)=1.-SUMY
            F(ICKMC(K0L0+KK)+ICELL)=RCKMC(K0BUF+KK)
            IF(SOLTEM) F(L0TEM+ICELL)=RCKMC(K0VAR+1)
C
C--- Load the print-out arrays for Rates-Of-Reaction, WDoTs &
C    Note that the first KK elements of the scratch array (K0SCR+...)
C    contain the rates-of-production (WDoTs)
C
            IF(L0QDOT>0) F(L0QDOT+ICELL)=QDOT
            IF(ICKMC(K0WDT+1)/=IGRND) THEN
              DO 81421 K=1,KK
                I=ICKMC(K0WDT+KK+K)
                IF(I>0) F(I+ICELL)=RCKMC(K0SCR+K)*RCKMC(K0MWT+K)
81421         CONTINUE
            ENDIF
            IF(ICKMC(K0ROR+1)/=IGRND) THEN
              CALL CKQYP(PRES,RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,
     1                   RCKMC(K0BUF+KK+1))
              DO 81422 K=1,II
                I=ICKMC(K0ROR+II+K)
                IF(I>0) F(I+ICELL)=RCKMC(K0BUF+KK+K)
81422         CONTINUE
            ENDIF
C
C--- Set ENUFSW on the basis of residual for the current cell on entry
C    to the solution procedure, and on the basis of the change in the
C    variables solved using PHOTWO for this iteration. If both are
C    small enough when compared with the RESREF and ENDIT variables
C    signal that ENUFSW need not be .FALSE.
C
            DO 81425 K=NY0+1,NATJ
              KIV=LBCKF+K-NY0-1
              IF(RESREF(KIV)<0) ENUFSW=ENUFSW.AND.
     1       ABS(RCKMC(K0SCR+2*KK+K))/(TINY+ABS(RESREF(KIV)))<1.0
              DELVAR=RCKMC(K0VAR+K)-RCKMC(K0VARN+K)
              ENUFSW=ENUFSW.AND.ABS(DELVAR)=LBCKF.AND.INDVAR<=LBCKF+KK-2.AND.
     1          ICKOPT==1) THEN
C
C--- If this is a CHEMKIN variable, but not the last, then just
C    save the SU and AP, zero the correction store, and set USER
C    to indicate that the PHOENICS solver should do nothing
C
          USER=.TRUE.
          CALL FN1(LD7,0.0)
          I=L0SUK+INDVAR-LBCKF+1
          ISA=L0APK-L0SUK
          L0SU=L0F(LSU)
          L0AP=L0F(LAP)
          DO 81451 ICELL=1,NXNY
C
            CALL SUB2L(BLOCKD,.FALSE.,SOLID,.FALSE.)
            IF(IPRPS>0) THEN
              SOLID =F(L0PRPS+ICELL)>=SOLPRP
              BLOCKD=F(L0PRPS+ICELL)>=PORPRP
            ENDIF
C           IF(TSTVPO) BLOCKD=BLOCKD.OR.F(L0VPOR+ICELL)<=1.E-10
            IF(.NOT.(BLOCKD.OR.SOLID)) THEN
C
C--- If the cell is completely blocked, do nothing
C
              F(I)=F(L0SU+ICELL)
              F(I+ISA)=F(L0AP+ICELL)
            ELSE
              F(I)=0.0
              F(I+ISA)=1.E10
            ENDIF
            I=I+KK
81451     CONTINUE
        ENDIF
      GO TO 999
  815 CONTINUE
C   * ------------------- SECTION 15 --- Change solution
      GO TO 999
C
C   * See the equivalent section in GREX for the indices to be
C     used in sections 7 - 15
C
C   * Make all other group-8 changes in GROUP 19.
C-----------------------------------------------------------------
C
C--- GROUP 9. Properties of the medium (or media)
C
C   The sections in this group are arranged sequentially in their
C   order of calling from EARTH. Thus, as can be seen from below,
C   the temperature sections (10 and 11) precede the density
C   sections (1 and 3); so, density formulae can refer to
C   temperature stores already set.
    9 CONTINUE
C-----------------------------------------------------------------
C
C--- The next section is done for all sections
C
      ICELL=ICELLA
      PRES=PRESS0
      RTMP=TEMP0
      IF(VARPRS) PRES=PRES+F(L0P1+ICELL)
      IF(VARTEM) RTMP=RTMP+F(L0TEM+ICELL)
      SUMY=0.0
      SUMW=0.0
C     IF(ICKOP1/=2.OR.IMIXF==0) THEN
      IF(STOCMP) THEN
        DO 9001 K=1,KK-1
          RCKMC(K0BUF+K)=F(ICKMC(K0L0+K)+ICELL)
          SUMW=SUMW+RCKMC(K0BUF+K)/RCKMC(K0MWT+K)
          SUMY=SUMY+RCKMC(K0BUF+K)
 9001   CONTINUE
      ENDIF
C
      RCKMC(K0BUF+KK)=1.-SUMY
      SUMW=SUMW+RCKMC(K0BUF+KK)/RCKMC(K0MWT+KK)
      IF(IMMWT>0) F(L0MMWT+ICELL)=1./SUMW
      IF(IGSH>0) IZOFF=IZOFF+NXNY
C--- Check whether we are in a solid (special actions taken)
C    probably done in GXPRPS
C
      IF(ISC==14) THEN
C   * ------------------- SECTION 14 ---------------------------
C    For SOLVE(TEMP1)-------------------------- phase-1 specific heat
C
C....
        CALL SUB2(L0HZ,L0H+IZOFF,L0CPZ,L0CP+IZOFF)
        CALL CKHMS(RTMP,ICKMC,RCKMC,RCKMC(K0BUF+KK+1))
        HMIX=0.0
        IF(IGSH>0) THEN
          CALL SUB2(K0L,L0HK,K0,L0HK+KK*NXNY)
          KADD=0
          DO 90001 K=1,KK
C            CALL FN0(-K0L-KADD,-K0-KADD)
            F(K0L+KADD+ICELL)=F(K0+KADD+ICELL)
            KADD=KADD+NXNY
90001     CONTINUE
        ELSE
          CALL SUB2(K0,L0HK,K0H,L0HK+KK*NXNY)
        ENDIF
C
C--- Calculate the sensible enthalpies, store them and calculate
C    the mixture enthalpy
C
 
        DO 9002 K=1,KK
          HSENS=RCKMC(K0BUF+KK+K)-RCKMC(K0H0K+K)
          F(K0+ICELL)=HSENS
          HMIX=HMIX+RCKMC(K0BUF+K)*HSENS
          K0=K0+NXNY
 9002   CONTINUE
c        VALOUT=HMIX/RTMP
        VALOUT=HMIX/f(l0tem+icell)
        F(L0HZ+ICELL)=HMIX
        IF(L0ENT>0) THEN
          HMIXT=0.0
          DO 9021 K=1,KK
            HMIXT=HMIXT+RCKMC(K0BUF+K)*RCKMC(K0BUF+KK+K)
 9021     CONTINUE
          F(L0ENT+ICELL)=F(L0DEN+ICELL)*HMIXT
        ENDIF
        CALL CKCPBS(RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,CPMIX)
        F(L0CPZ+ICELL)=CPMIX
        IF(IGSH==0.AND.NZ>1) THEN
          CALL SUB3(K0,L0HK,K0H,L0HK+KK*NXNY,KADD,0)
          DO 90011 K=1,KK
C            CALL FN0(-K0H-KADD,-K0-KADD)
            F(K0H+KADD+ICELL)=F(K0+KADD+ICELL)
            KADD=KADD+NXNY
90011     CONTINUE
        ENDIF
C
        if(dbgloc.and.dbt) then
          call writ3i('   icell',icell,'  isweep',isweep,
     1                '  izstep',izstep)
          call writ2r('    hmix',real(hmix),'   cpmix',real(cpmix))
          write(lupr1,'('' individual enthalpies are:''/1p5e12.5)')
     1        (rckmc(k0buf+kk+k),k=1,kk)
        endif
C
      ELSEIF(ISC==1) THEN
C   * ------------------- SECTION  1 ---------------------------
C    For RHO1<=GRND--- density for phase 1       Index DEN1
C
        CALL CKRHOY(PRES,RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,RHOMIX)
        VALOUT=RHOMIX
C
        if(dbgloc.and.dbrho) then
          call writ3i('   icell',icell,'  isweep',isweep,
     1                '  izstep',izstep)
          call writ1r('     rho',valout)
        endif
C
      ELSEIF(ISC==6) THEN
C   * ------------------- SECTION  6 ---------------------------
C    For ENUL<=GRND--- reference laminar kinematic viscosity
C                                                  Index VISL
C
        VALOUT=0.0
        IF(DODIFF) THEN
          VISMIX=0.0
          CALL CKYTX(RCKMC(K0BUF+1),ICKMC,RCKMC,RCKMC(K0BUF+1))
          CALL MCAVIS(RTMP,RCKMC(K0BUF+1),RCKMC(KMCR),VISMIX)
          VALOUT=VISMIX/F(L0DEN+ICELL)
          IF(SIUNIT) VALOUT=VALOUT*VISSCA
C
          if(dbgloc.and.dbemu) then
            call writ3i('   icell',icell,'  isweep',isweep,
     1                  '  izstep',izstep)
            call writ1r('     enu',valout)
          endif
        ENDIF
C
      ELSEIF(ISC==7) THEN
C   * ------------------- SECTION  7 ---------------------------
C    For PRNDTL( )<=GRND--- laminar PRANDTL nos., or diffusivity
C                                                  Index LAMPR
        IF(DODIFF) THEN
          CALL CKYTX(RCKMC(K0BUF+1),ICKMC,RCKMC,RCKMC(K0BUF+1))
          IF(INDVAR==ITIND) THEN
            CALL MCACON(RTMP,RCKMC(K0BUF+1),RCKMC(KMCR),CONMIX)
            VALOUT=CONMIX
            IF(SIUNIT) VALOUT=VALOUT*CONSCA
            IF(L0CON>0) F(L0CON+ICELL)=VALOUT
C
            if(dbgloc.and.dbgam) then
              call writ3i('   icell',icell,'  isweep',isweep,
     1                  '  izstep',izstep)
              call writ1r('    cond',valout)
            endif
C
          ELSEIF(IMCOPT<2) THEN
            IF(INDVAR==LBCKF) THEN
        if(dbgloc.and.dbgam) then
          call writ3i('   icell',icell,'  isweep',isweep,
     1                '  izstep',izstep)
          write(lupr1,'('' mole fractions are:''/1p5e12.5)')
     1        (rckmc(k0buf+k),k=1,kk)
        endif
              CALL MCADIF(PRES,RTMP,RCKMC(K0BUF+1),RCKMC(KMCR),
     1                    RCKMC(K0BUF+KK+1))
              K0=L0DK+IZOFF*KK+ICELL
              DO 9003 K=1,KK
                F(K0)=RCKMC(K0BUF+KK+K)
                K0=K0+NXNY
 9003         CONTINUE
              IF(SIUNIT) THEN
                K0=L0DK+IZOFF*KK+ICELL
                DO 90035 K=1,KK
                  F(K0)=F(K0)*DIFSCA
                  K0=K0+NXNY
90035           CONTINUE
              ENDIF
            ENDIF
            K0=L0DK+IZOFF*KK
            KIV=INDVAR-LBCKF
            VALOUT=F(K0+KIV*NXNY+ICELL)
            IF(MODDIF) THEN
              VALOUT=VALOUT*SUMW
            ENDIF
C
            if(dbgloc.and.dbgam) then
              call writ3i('   icell',icell,'  isweep',isweep,
     1                    '  izstep',izstep)
              write(lupr1,'('' diffusion coefficients:''/1p5e12.5)')
     1           (rckmc(k0buf+kk+k),k=1,kk)
              call writ1l(' moddif ',moddif)
cccc              call writ2r(' sumw   ',sumw,' pres   ',pres)
              call writ1r(' valout ',valout)
            endif
C
            IF(DOTDIF) THEN
C
C--- Calculate the thermal diffusion ratios and store in L0DTK
C
              IF(INDVAR==LBCKF) THEN
                CALL MCATDR(RTMP,RCKMC(K0BUF+1),ICKMC(KMCI),
     1                      RCKMC(KMCR),RCKMC(K0BUF+KK+1))
                K0  =L0DTK+IZOFF*KK+ICELL
                DO 9004 K=1,KK
                  F(K0)=RCKMC(K0BUF+KK+K)
                  K0=K0+NXNY
 9004           CONTINUE
C
                if(dbgloc.and.dbgam) then
                  call writ3i('   icell',icell,'  isweep',isweep,
     1                    '  izstep',izstep)
                  write(lupr1,'('' thermal diffusion ratios:''
     1               /1p5e12.5)') (rckmc(k0buf+kk+k),k=1,kk)
                endif
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      GO TO 999
C-----------------------------------------------------------------
C
C--- GROUP 11. Initialization of variable or porosity fields
C                                                  Index VAL
   11 CONTINUE
C    For RHO1<=GRND--- density for phase 1
C
      IF(NPATCH(1:6)=='ICHEMK') THEN
        L0VAL=L0F(VAL)
        IF(INDVAR==DEN1) THEN
          PRES=PRESS0
          RTMP=TEMP0
          DO 11011 ICELL=1,NXNY
            IF(VARPRS) PRES=PRESS0+F(L0P1+ICELL)
            IF(VARTEM) RTMP=TEMP0+F(L0TEM+ICELL)
            DO 11001 K=1,KK
              RCKMC(K0BUF+K)=F(ICKMC(K0L0+K)+ICELL)
11001       CONTINUE
            CALL CKRHOY(PRES,RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,RHOMIX)
            F(L0VAL+ICELL)=RHOMIX
11011     CONTINUE
        ELSEIF(INDVAR==W1) THEN
          CALL POISSL(KVEL,IMUL)
        ENDIF
      ENDIF
C
      GO TO 999
C-----------------------------------------------------------------
C
C--- GROUP 12. Convection and diffusion adjustments
C
c   12 CONTINUE
c      GO TO 999
C-----------------------------------------------------------------
C
C--- GROUP 13. Boundary conditions and special sources
C                                       Index for Coefficient - CO
C                                       Index for Value       - VAL
   13 CONTINUE
C
      IF(NPATCH(1:5)/='CHEMK'.AND.NPATCH(1:4)/='NOCP'.AND.
     2   NPATCH(1:6)/='CHEMID'.AND.NPATCH(1:5)/='CKWAL') THEN
        CALL WRIT40(' GXCHKI - Should not be here for PATCH  ')
        WRITE(LUPR1,'(''  NAME = '',A8)') NPATCH
        CALL SET_ERR(268,
     *        'Error. GXCHKI - Should not be here for PATCH.',1)
C        CALL EAROUT(1)
      ENDIF
C
C!!!! Note to users: The next statement indicates how the Eddy-Break-Up
C     model MIGHT be implemented; but it has not been tested. Users
C     wishing to work with it should "uncomment".
c      IF(NPATCH(1:6)=='CHEMK1'.OR.
c     1   NPATCH(1:6)=='CHEMK2'.OR.
c     1   NPATCH(1:6)=='CHEMK3') THEN
C
      IF(NPATCH(1:6)=='CHEMK1') THEN
C
C--- Non-stiff / slightly-stiff sources
C     CKCTYP - calculates rates of creation and characteristic
C              destruction time-scales
C     CKWYP  - calculates rates of formation of species
C     CKCTYR - calculates rates of creation and characteristic
C              timescales for destruction of species
C
        CALL SUB2(L0HVZ,L0HVAL+IZOFF*KK,L0HCZ,L0HCO+IZOFF*KK)
        IF(ISC<11) THEN
          IF(INDVAR==LBCKF) THEN
            DO 1301 K=1,KK
              F(L0DTF+K)=DTFALS(LBCKF+K-1)
 1301       CONTINUE
            IF(SOLTEM) F(L0DTF+KK+1)=DTFALS(ITIND)
            CALL L0F1(CO,ICELL,IADD,'GXCHKI')
            ICELL=ICELL+CO
            PRES=PRESS0
            RTMP=TEMP0
            DO 1305 IX=IXF,IXL
              ICELL=ICELL+IADD
              DO 1304 IY=IYF,IYL
                ICELL=ICELL+1
                BLOCKD=.FALSE.
                IF(IPRPS>0) BLOCKD=F(L0PRPS+ICELL)>=SOLPRP
C               IF(TSTVPO) BLOCKD=BLOCKD.OR.F(L0VPOR+ICELL)<=1.E-10
                IF(BLOCKD) THEN
                  K0C=L0CDK+IZOFF*KK
                  K0T=L0TDK+IZOFF*KK
                  DO 1391 K=1,KK
                    F(K0C+ICELL)=0.0
                    F(K0T+ICELL)=0.0
                    K0C=K0C+NXNY
                    K0T=K0T+NXNY
 1391             CONTINUE
                  F(L0HVZ+ICELL)=0.0
                  F(L0HCZ+ICELL)=0.0
                  GO TO 1304
                ENDIF
                IF(VARPRS) PRES=PRES+F(L0P1+ICELL)
                IF(VARTEM) RTMP=TEMP0+F(L0TEM+ICELL)
                SUMY=0.0
                SUMW=0.0
                DO 1302 K=1,KK-1
                  RCKMC(K0VAR+K)=F(ICKMC(K0L0+K)+ICELL)
                  RCKMC(K0BUF+K)=0.0
                  RCKMC(K0BUF+KK+K)=0.0
                  SUMY=SUMY+RCKMC(K0VAR+K)
                  SUMW=SUMW+RCKMC(K0VAR+K)/RCKMC(K0MWT+K)
 1302           CONTINUE
                RCKMC(K0VAR+KK)=1.-SUMY
                SUMW=SUMW+RCKMC(K0VAR+KK)/RCKMC(K0MWT+KK)
                IF(SOLTEM) THEN
                  IF(ICKOP1==1.OR.ICKOP1==3) THEN
c---- calculate wdotk, i.e. molar production rates of species k
                    CALL CKWYP(PRES,1.001*RTMP,RCKMC(K0VAR+1),
     1                         ICKMC,RCKMC,RCKMC(K0BUF+KK+1))
                  ENDIF
C
                  HSORP=0.0
C---- calculate qdot = sum(wdotk*wtk*h0k) over all k
                  DO 1306 K=1,KK
                    HSORP=HSORP+RCKMC(K0BUF+KK+K)*RCKMC(K0MWT+K)*
     1                          RCKMC(K0H0K+K)
 1306             CONTINUE
                  IF(RELQDT) HSORP=ALQDOT*HSORP-
     1                                     (1.-ALQDOT)*F(L0QDOT+ICELL)
                ENDIF
                IF(ICKOP1==1.OR.ICKOP1==3) THEN
                  GRHO=F(L0DEN+ICELL)
                  CALL CKCTYR(GRHO,RTMP,RCKMC(K0VAR+1),ICKMC,RCKMC,
     1                        RCKMC(K0BUF+1),RCKMC(K0BUF+KK+1))
C                  CALL CKCTYP(PRES,RTMP,RCKMC(K0VAR+1),ICKMC,RCKMC,
C     1                        RCKMC(K0BUF+1),RCKMC(K0BUF+KK+1))
                ENDIF
                RHO=F(L0DEN+ICELL)
                HSOR=0.0
                K0C=L0CDK+IZOFF*KK
                K0T=L0TDK+IZOFF*KK
                DTFT=F(L0DTF+KK+1)
                DO 1303 K=1,KK
                  CDOTK=RCKMC(K0BUF+K)
                  DTF=F(L0DTF+K)
                  IF(RCKMC(K0BUF+KK+K)>GREAT.OR.
     1               RCKMC(K0BUF+KK+K)<=0.0) THEN
                    TAUK = GREAT
                  ELSE
                    TAUK =RCKMC(K0BUF+KK+K)
                  ENDIF
                  IF(TAUK>1.E-20.OR.CDOTK>1.E-20) THEN
                    RHODT=RHO/TAUK
                    YCONR=RHO/(SUMW*RCKMC(K0MWT+K))
                    TAUCK=YCONR*(RCKMC(K0VAR+K)+1.E-4)/(TINY+CDOTK)
                    DTF=MIN(DTF,TAUCK)
                    DTF=MIN(DTF,TAUK)
                    F(L0DTF+K)=DTF
                    F(K0C+ICELL)=RHODT
                    F(K0T+ICELL)=CDOTK*RCKMC(K0MWT+K)/RHODT
                    IF(SOLTEM) HSOR=HSOR+RHODT*RCKMC(K0H0K+K)*
     1                           (F(K0T+ICELL)-F(ICKMC(K0L0+K)+ICELL))
                  ELSE
                    F(K0C+ICELL)=0.0
                    F(K0T+ICELL)=0.0
                  ENDIF
                  K0C=K0C+NXNY
                  K0T=K0T+NXNY
 1303           CONTINUE
                IF(SOLTEM) THEN
                  IF(RELQDT) HSOR=ALQDOT*HSOR-
     1                                     (1.-ALQDOT)*F(L0QDOT+ICELL)
                  DHSBDT=(HSORP-HSOR)/(0.001*RTMP)
                  if(hsor>0.0) then
                    IF(DHSBDT<0.0) THEN
                      F(L0HVZ+ICELL)=RTMP-HSOR/DHSBDT
                      F(L0HCZ+ICELL)=-DHSBDT
                    ELSE
                      F(L0HVZ+ICELL)=-HSOR/FIXFLU
                      F(L0HCZ+ICELL)=FIXFLU
                    ENDIF
                  else
                    IF(DHSBDT>0.0) THEN
                      F(L0HVZ+ICELL)=RTMP-HSOR/DHSBDT
                      F(L0HCZ+ICELL)=DHSBDT
                    ELSE
                      F(L0HVZ+ICELL)=-HSOR/FIXFLU
                      F(L0HCZ+ICELL)=FIXFLU
                    ENDIF
                  endif
                ENDIF
                IF(L0QDOT>0) F(L0QDOT+ICELL)=-HSOR
 1304         CONTINUE
 1305       CONTINUE
          ENDIF
          IF(INDVAR==ITEM) THEN
            CALL FN0(CO,-L0HCZ)
          ELSEIF(INDVAR>=LBCKF.AND.INDVAR0.0) CALL FN34(LAP,LM1,RECT)
          ENDIF
        ELSE
          IF(INDVAR==ITEM) THEN
            CALL FN0(VAL,-L0HVZ)
          ELSEIF(INDVAR>=LBCKF.AND.INDVAR=U1.AND.(EQZ(VCO)
     1     .OR.VCO>=0.8*FIXVAL).OR.
     1     INDVAR==P1.AND.QEQ(PCO,FIXFLU)) THEN
C
C--- Setting of inlet velocity: get incoming mass-flux
C
          IF(GRN(PCO)) THEN
            CALL WRIT40(' GXCHKI - Cannot do INLET with GRND P-CO')
            CALL SET_ERR(269,
     *        'Error. GXCHKI - Cannot do INLET with GRND P-CO.',1)
C            CALL EAROUT(1)
          ELSEIF(GRN(RHO1)) THEN
C
C--- Get the incoming gas composition and temperature to calculate the
C    density of the incoming gas: note that the incoming temperature
C    for GRND set enthalpy fluxes comes form TINLET
C
            CALL GETCOV(NPATCH,ITIND,TCO,TVAL)
            IF(GRN(TVAL)) THEN
              CALL GETSPD(NPATCH,'TINLET',1,TINLET,IDUM,.FALSE.,
     1                    '   ',IERR)
              IF(IERR==1) TINLET=TMP1A
                RTMP=TEMP0+TINLET
            ELSE
              RTMP=TEMP0+TVAL
            ENDIF
            LCKERR=.FALSE.
            SUMY=0.0
            DO 13211 K=1,KK-1
              CALL GETCOV(NPATCH,LBCKF+K-1,YCO,YVAL)
              LCKERR=LCKERR.OR.GRN(YVAL)
              RCKMC(K0BUF+K)=YVAL
              SUMY=SUMY+YVAL
13211       CONTINUE
            RCKMC(K0BUF+KK)=1.-SUMY
            IF(LCKERR) THEN
              CALL WRIT40(' GXCHKI - Cannot do INLET with GRND Y-CO')
              CALL SET_ERR(270,
     *        'Error. GXCHKI - Cannot do INLET with GRND Y-CO.',1)
C              CALL EAROUT(1)
            ELSE
              PRES=PRESS0
C     inlet species defined in mol fractions
              IF(NPATCH(8:8)=='X') THEN
                CALL CKRHOX(PRES,RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,RHOMIX)
              ELSE
                CALL CKRHOY(PRES,RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,RHOMIX)
              ENDIF
            ENDIF
            RHOIN=REAL(RHOMIX)
          ELSE
            RHOIN=RHO1
          ENDIF
C
C--- Now look to see if we have a Poisseulle flow PATCH; this is
C    indicated by the 6th letter in the PATCH name, by having
C    a non-zero value of the relevant PROFx variable, and by
C    having the right PATCH-type etc. It is assumed that if the
C    velocity is set as Poisseulle flow, then so also is the
C    mass-flux (through P1). The work is all done in the subroutine
C    POISSL
C
          IF(BFC.AND.ISC==13) THEN
            IF(INDVAR==P1.OR.(INDVAR>=U1.AND.INDVAR<=W1)) THEN
              CALL FN1(VAL,RHOIN)
              CALL GXBFC
            ENDIF
          ELSE
            CALL POISSL(KVEL,IMUL)
C
            IF(INDVAR==P1) THEN
C
C--- If this is the pressure, we are determining the mass-flux from
C    the magnitude of the velocity multiplied by the density of the
C    inflowing gas (NB.: the calculation is incorrect for BFC grids)
C
              IF(KVEL/=W1) THEN
                CALL FN1(VAL,0.0)
                DO 13221 IVEL=U1,W1,2
                  CALL GETCOV(NPATCH,IVEL,VCO,VVALA)
                  IF(GRN(VVALA)) THEN
                  ELSEIF(NEZ(VVALA)) THEN
                    CALL FN33(VAL,VVALA*VVALA)
                  ENDIF
13221           CONTINUE
                CALL FN30(VAL)
              ENDIF
              CALL FN25(VAL,RHOIN)
            ELSEIF(.NOT.GRN(PVAL)) THEN
C
C--- If this is a velocity, we are determining the inflow-velocity
C    from the mass-flux and the inflow density (note that this assumes
C    only 1 non-zero inflow velocity)
C
              CALL FN1(VAL,PVAL/RHOIN)
            ELSE
              CALL FN25(VAL,FLOAT(IMUL))
            ENDIF
          ENDIF
        ELSEIF(INDVAR==ITIND) THEN
C
C--- Set the value of the incoming enthalpy NOCP PATCHes; get the
C    incoming temperature form TINLET (always have GRND set enthapy
C    flux for this option)
C
          LCKERR=.FALSE.
          CALL GETSPD(NPATCH,'TINLET',1,TINLET,IDUM,.FALSE.,
     1                '   ',IERR)
          IF(IERR==1) TINLET=TMP1A
          RTMP=TEMP0+TINLET
          SUMY=0.0
          DO 13212 K=1,KK-1
            CALL GETCOV(NPATCH,LBCKF+K-1,YCO,YVAL)
            LCKERR=LCKERR.OR.GRN(YVAL)
            RCKMC(K0BUF+K)=YVAL
            SUMY=SUMY+YVAL
13212     CONTINUE
          RCKMC(K0BUF+KK)=1.-SUMY
          IF(LCKERR) THEN
            CALL WRIT40(' GXCHKI - Cannot do INLET with GRND Y-CO')
              CALL SET_ERR(271,
     *        'Error. GXCHKI - Cannot do INLET with GRND Y-CO.',1)
C            CALL EAROUT(1)
          ELSE
            CALL CKHBMS(RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,CPMIX)
C
C--- Must subtract off the heat of formation part of the enthalpy
C    (Note: CPMIX is used to store the enthalpy here)
C
            DO 13213 K=1,KK
              CPMIX=CPMIX-RCKMC(K0BUF+K)*RCKMC(K0H0K+K)
13213       CONTINUE
            CALL FN1(VAL,REAL(CPMIX))
          ENDIF
        ENDIF
C
      ELSEIF(NPATCH(1:6)=='CHEMID') THEN
C
C--- Inlet diffusion PATCHes - the enthalpy flux due to mass-diffusion
C    is summed and the added into the enthalpy source
C
        IF(INDVAR>=LBCKF.AND.INDVAR<=LBCKF+KK-2.AND.
     1                                          ISC==10) THEN
          IPV=L0PVAR(PVCKID,IREG,IZSTEP)
          IF(INDVAR==LBCKF) CALL FNSCPA(-IPV,0.0)
          CALL GETCOV(NPATCH,INDVAR,YCO,YVAL)
          CALL GETSPD(NPATCH,'TINLET',1,TINLET,IDUM,.FALSE.,
     1                '   ',IERR)
          IF(IERR==1) TINLET=TMP1A
          CALL L0F2(CO,INDVAR,ICO,IFIC,IADD,'GXCHKI')
          DO 13301 IX=IXF,IXL
            ICO=ICO+IADD
            DO 13301 IY=IYF,IYL
              ICO=ICO+1
              IPV=IPV+1
              RTMP=TEMP0+TINLET
              CALL CKHMS(RTMP,ICKMC,RCKMC,RCKMC(K0BUF+1))
              F(IPV)=F(IPV)+MAX(F(ICO)*(YVAL-F(ICO+IFIC)),0.0)
     1                            *(RCKMC(K0BUF+K)-RCKMC(K0BUF+KK))
13301     CONTINUE
        ELSEIF(INDVAR==ITIND.AND.ISC==21) THEN
          CALL L0F1(VAL,IVAL,IADD,'GXCHKI')
          IPV=L0PVAR(PVCKID,IREG,IZSTEP)
          DO 13302 IX=IXF,IXL
            IVAL=IVAL+IADD
            DO 13302 IY=IYF,IYL
              IVAL=IVAL+1
              IPV=IPV+1
              F(IVAL)=F(IPV)
13302     CONTINUE
        ENDIF
C
      ELSEIF(MODDIF.AND.NPATCH(1:5)=='CKWAL'.AND.ISC==10.AND.
     1        INDVAR>=LBCKF.AND.INDVAR=LBCKF.AND.INDVAR
      SUBROUTINE PHINDX(NGO,LINK,LOUT,PRECIS,NTDA,K0TD)
C
C---------------------------------------------------------------
C  PHINDX - Checks that the precision of GXCHKI matches that of
C           the CHEMKIN and TRANLIB data-bases; also sets various
C           indices for GXCHKI
C
C  returns .TRUE. if there are thermal diffusion
C           coeff.s
C---------------------------------------------------------------
C
C*****double precision
      DOUBLE PRECISION
C*****END double precision
C*****single precision
C      REAL
C*****END single precision
     1     RU, PATMOS, SMALL
      CHARACTER*40 PRECIS
      CHARACTER*16 PREC, VERS, PRECMC, VERSMC
      LOGICAL ERROR, KERR, KERRMC
      COMMON /MCCONS/ VERSMC, PRECMC, KERR, LENI, LENR
      COMMON /CKCONS/ VERS, PREC, KERRMC, LENIMC, LENRMC
C
      COMMON /MCMCMC/ RU, PATMOS, SMALL, NKK, NO, NLITE, INLIN, IKTDIF,
     1                IPVT, NWT, NEPS, NSIG, NDIP, NPOL, NZROT, NLAM,
     2                NETA, NDIF, NTDIF, NXX, NVIS, NXI, NCP,
     3                NCROT, NCINT, NPARK, NBIND, NEOK, NSGM,
     4                NAST, NBST, NCST, NXL, NR, NWRK, K3
C
      IF(NGO==1) THEN
C
C--- Check CHEMKIN data-base
C
        ERROR = PREC(1:6)/=PRECIS(2:7)
C
      ELSEIF(NGO==2) THEN
C
C--- Check the TRANLIB data-base; get the number of species for which
C    there are thermal diffusion ratios/coefficients
C
        ERROR = PREC(1:6)/=PRECIS(2:7)
C
        NTDA=NLITE
        K0TD=IKTDIF-1
C
      ENDIF
C
      END
C
C===============================================================
C
      SUBROUTINE POISSL(KVEL,IMUL)
C
C-----------------------------------------------------------------
C SUBROUTINE POISSL - puts a POISSL flow into the VAL slot for
C group 13 or group 11 PATCHes
C-----------------------------------------------------------------
C Arguments:
C      KVEL - returns a code that indicates whether or not a
C             setting has been made
C-----------------------------------------------------------------
C
      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
      LOGICAL GRN
C
      IMUL=0
      KVEL=0
      IF((INDVAR==P1.AND.IGR==13).OR.INDVAR==W1) THEN
        VPOISS=GRND
        IF(NPATCH(7:7)=='A') THEN
          VPOISS=PROFA
        ELSEIF(NPATCH(7:7)=='B') THEN
          VPOISS=PROFB
        ELSEIF(NPATCH(7:7)=='C') THEN
          VPOISS=PROFC
        ENDIF
        KVEL=0
        IF(.NOT.GRN(VPOISS)) THEN
          IF(IGR==13.AND.(INTTYP==6.OR.INTTYP==7).OR.
     1       IGR==11.AND.IPATNF(2,IREG)==26) THEN
            KVEL=7
          ENDIF
          IF(IGR==13) THEN
            IGEO=PATGEO
          ELSEIF(IGR==11) THEN
            IGEO=ANORTH
          ENDIF
          IF(VPOISS>0.0) THEN
            IDVP=YG2D
            IF(IYF==1) THEN
              XYF=0.0
            ELSE
              XYF=F(L0F(YV2D)+IYF-1)
            ENDIF
            XYL=F(L0F(YV2D)+IYL)
            IF(.NOT.CARTES) THEN
              XYF=XYF+RINNER
              XYL=XYL+RINNER
            ENDIF
          ELSE
            VPOISS=-VPOISS
            IDVP=XG2D
            IF(IXF==1) THEN
              XYF=0.0
            ELSE
              XYF=F(L0F(XU2D)+(IXF-2)*NY)
            ENDIF
            XYL=F(L0F(XU2D)+(IXL-1)*NY)
          ENDIF
          IF(IGR==13) IMUL=2*MOD(INTTYP,2)-1
          SUM=0.0
          ASUM=0.0
          CALL L0F3(VAL,IDVP,IGEO,I,IDMI,IGMI,IADD,'GXCHKI')
          IF(XYF>0.0) RLRAT=(XYL**2-XYF**2)/LOG(XYL/XYF)
          DO 13215 IX=IXF,IXL
            I=I+IADD
            DO 13215 IY=IYF,IYL
              I=I+1
              F(I)=XYL**2-F(I+IDMI)**2
              IF(XYF>0.0) F(I)=F(I)+RLRAT*LOG(F(I+IDMI)/XYL)
              SUM=SUM+F(I)*F(I+IGMI)
              ASUM=ASUM+F(I+IGMI)
13215     CONTINUE
          CALL FN25(VAL,VPOISS*ASUM/SUM)
        ENDIF
      ENDIF
      END
c