gxcvd

cgxcvd

C FILE NAME GXCVD.FOR--------------------------------------------- 110121
      SUBROUTINE GXCVD
C --------------------------------------------------------------------
      INCLUDE 'patcmn'
      INCLUDE 'farray'
      INCLUDE 'cmncvd'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE 'grdbfc'
      COMMON/SLDPRP/SOLPRP,PORPRP,VACPRP
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION STARTS:
C
C 1   Set dimensions of data-for-GROUND arrays here. WARNING: the
C     corresponding arrays in the MAIN program of the satellite
C     and EARTH must have the same dimensions.
      PARAMETER (NLG=100, NIG=200, NRG=200, NCG=100)
C
      COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
      LOGICAL LG,WHF,EQZ,QEQ
      CHARACTER*4 CG
C     rname: used in calls for special data from Q1
      CHARACTER*15 RNAME
C
C 2   User dimensions own arrays here, for example:
C     DIMENSION GUH(10,10),GUC(10,10),GUX(10,10),GUZ(10)
C --------------------------------------------------------------------
C PARAMETERS:
      REAL          GR
      PARAMETER     (GR=8.314413)
C     diffusion coefficient
      REAL          GDIF
C     mean molecular weight
      REAL          GMMW
      REAL          GSUM
C     temperature
      REAL          GTEMP
C     pressure
      REAL          GPRES
C zero f-array indices
      INTEGER       L0TEM
      INTEGER       L0P
      INTEGER       L0C1
      INTEGER       L0VPOR
      INTEGER       L0DEN
      INTEGER       L0CO,L0VAL
      INTEGER       L0CP
      INTEGER       L0DEPO
C.........................................................
C      character*2 gel
C     Chemistry data:
      INTEGER       JHOMAX
      PARAMETER     (JHOMAX=100)
C     number of gas phase reactions:
      INTEGER       NGREAC
C     integer array of gas phase reaction numbers
      INTEGER       GREAC(JHOMAX)
C     stoichiometric coefficients:
      INTEGER       IGAMMA(JMAXS,JHOMAX)
C     st. coeff. summed over all species:
      INTEGER       JGAMMA(JHOMAX)
C...........................................................
C     Chemistry:
      REAL          GC1
      REAL          GHKMOL
      REAL          GHKMAS
      REAL          GCP
      REAL          GCDOT(JMAXS)
      REAL          GSOURCE
      REAL          FPRT(JMAXS)
      REAL          RF(JHOMAX)
      REAL          RR(JHOMAX)
      INTEGER       JR
      REAL          GFTS
      REAL          GSP(JMAXS)
      REAL          GSP1
      REAL          GSN(JMAXS)
      REAL          GSN1
      INTEGER       JDUM(JMAXS)
      INTEGER       JNUM(JMAXS)
      INTEGER       JGAS(JMAXS)
      INTEGER       JADS(JMAXS)
      CHARACTER*2   JPHSE(JMAXS)
      CHARACTER*2   JPHSV(JIND)
C...........................................................
C     Surface Chemistry:
      INTEGER       JSUMAX,JTT
      PARAMETER     (JSUMAX=100)
      INTEGER       NSREAC
      INTEGER       SREAC(JSUMAX)
      REAL          GDSOL
      REAL          GMSOL
      REAL          GSRCN
      CHARACTER*15  GHSOL
      INTEGER       GALFA(JMAXS,JSUMAX)
      INTEGER       GSALFA(JSUMAX)
      INTEGER       JDEP(JSUMAX)
      REAL          GFAC
C.................................................................
C     Surface chemistry sourceterm:
      REAL          GRKIN(JSUMAX)
      REAL          GDIST
      REAL          GRDIF(JMAXS)
      REAL          GDIFW(JMAXS)
      REAL          GREFF(JSUMAX)
      REAL          GSURFL
      REAL          GZ(1000)
C
C Additional storage for data-base code
      PARAMETER(NIWS=30,NRWS=90,NCWS=30)
      INTEGER       IWS(NIWS),KGSP(JSUMAX),IGRT(JHOMAX),ISRT(JSUMAX)
      REAL          RWS(NRWS),GRPAR(4,JHOMAX),GAPAR(90,JHOMAX),
     1              GFTHB(JMAXS,JHOMAX),SRPAR(4,JHOMAX),
     2              SAPAR(90,JHOMAX),SFTHB(JSUMAX,JHOMAX)
      CHARACTER*16  CWS(NCWS)
      LOGICAL       ERROR
      LOGICAL LINIL
C.... Set dimension of patch-name array here.  WARNING: the array
C     NAMPAT in the MAIN program of the satellite, and EARTH must have
C     the same dimension.
c      PARAMETER (NPNAM=5000)
c      COMMON/NPAT/NAMPAT(NPNAM)
      CHARACTER*8 PNAME
      COMMON/ICOVL/M04,I0PHI
C     ggnms : names of gas phase reactions
      CHARACTER*32  GGNMS(JHOMAX)
C     gsnms : names of surface reactions
      CHARACTER*32  GSNMS(JSUMAX)
      CHARACTER*15  CDUM
C     gmafr : mass fractions
      REAL          GMAFR(JMAXS)
C     gsocc : number of bonds per adsorbate
      INTEGER       GSOCC(JMAXS)
C     gsfra : site fractions of adsorbates
      REAL          GSFRA(JMAXS)
C     gdepr : deposition rate
      REAL          GDEPR
C     gsflux: : fluxes of gas species at surface
      REAL          GSFLUX(JMAXS)
C     l0fsu : zero f-array index for surface fluxes
      INTEGER       L0FSU
C     l0socc : zero f-array index for surface fractions
      INTEGER       L0SOCC
C     gmolw : molecular weights
      REAL          GMOLW(JMAXS)
C     is a reaction from the Arora model used?
      LOGICAL       LMODEL
C     are all reactions from the Arora?
      LOGICAL       LMODEL2
C     index of vacancy fraction
      INTEGER       KVAC
C     jtot: total number of species
      INTEGER       JTOT
C     gtrial: number of Newton iterations
      INTEGER       GTRIAL
C     gsuc, gapc: phoenics solver coefficients
      REAL          GSUC(JMAXS)
      REAL          GAPC(JMAXS)
C     antol, rntol: tolerance parameters for newton solver
      REAL          ANTOL(JMAXS)
      REAL          RNTOL(JMAXS)
C     gxvol: cell volume
      REAL          GXVOL
C     gabsol, grelat : machine accuracy parameters for Newton solver
      REAL          GABSOL, GRELAT
C     grho : density
      REAL          GRHO
C     istiff : integer flag : 0=phoenics solver; 1=Newton solver PBP
      INTEGER       ISTIFF
C     anjac: logical flag : f=numerical jacobian; t=analytical jacobian
      LOGICAL       ANJAC
C     thmdif: soret t/f
      LOGICAL       THMDIF
C     thmopt: option for soret calculation
      INTEGER       THMOPT
C     thmfrq: calculation frequency of soret
      INTEGER       THMFRQ
C     thmrlx: relaxation for soret calculation
      REAL          THMRLX
C     smrlx: relaxation applied to stefan-maxwell calculation
      REAL          SMRLX
C     mcdopt: multicomponent diffusion option
      INTEGER       MCDOPT
C     mcprop: multicomponent properties selection
      INTEGER       MCPROP
      COMMON/GENI/NXNY,IFIL1(10),IGSH,IFIL2(48)
      COMMON/GENL/LDUM1(15),SWPRIN,LDUM2(44)
      LOGICAL LDUM1,SWPRIN,LDUM2
      DIMENSION GLOWER(2*JMAXS,2*JMAXS),GUPPER(2*JMAXS+1,2*JMAXS+1),
     1          GDUMMY(2*JMAXS,2*JMAXS)
      INTEGER IORD(JMAXS)
      CHARACTER*15 CTEMP(JMAXS)
C --------------------------------------------------------------------
C
C 3   User places his data statements here, for example:
C     DATA NXDIM,NYDIM/10,10/
C
C 4   Insert own coding below as desired, guided by GREX examples.
C     Note that the satellite-to-GREX special data in the labelled
C     COMMONs /RSG/, /ISG/, /LSG/ and /CSG/ can be included and
C     used below but the user must check GREX for any conflicting
C     uses. The same comment applies to the EARTH-spare working
C     arrays EASP1, EASP2,....EASP20. In addition to the EASPs,
C     there are 10 GRound-earth SPare arrays, GRSP1,...,GRSP10,
C     supplied solely for the user, which are not used by GREX. If
C     the call to GREX has been deactivated then all of the arrays
C     may be used without reservation.
C
c***********************************************************************
c
      SAVE L0GC,L0CSP3,L0CSN3,L0CSPT,L0CSNT,L0ARL,L0TR,L0ST,L0JTE,
     1     L0JTN,L0JTH,L0JE,L0TE,L0UE,L0JN,L0TN,L0UN,L0JH,L0TH,L0UH,
     1     L0UL,L0SHWR,L0SURF,L0SUK,L0APK,L0RF,L0RR,L0CO
      SAVE JZL,JZH,L0TEM,L0PRP,L0PRPH,L0PRPL,L0VPOR,L0VPOH,L0VPOL,
     1     L0DEN,L0P,L0CP,JK1,JK2,JT1,JT2,L0SU,L0AP,L0T0,L0ION
      SAVE JTEM,JPRPS,JSPHT,JCOND,JT0,JNE,JION,JVPOR,JDEPO,JEPOR,JNPOR,
     1     JHPOR
      SAVE SMRLX,THMDIF,THMOPT,THMFRQ,THMRLX,NGREAC,GREAC,NSREAC,SREAC,
     1     CHMRLX,ANJAC
      SAVE GALFA,GSALFA,ISRT,KGSP,SAPAR,SRPAR,GAPAR,GRPAR,JDEP,GFTHB,
     1     IGAMMA,JGAMMA,GDSOL,GMSOL,GFAC,GSNMS,GSRCN,ISTIFF,FPRT,IGRT,
     1     JCHECK,LINIL,ISHWR,ISURF,KDEP,KVAC
C Redefinitions to cope with Group 9 index shifting
C L0TEM added to avoid problems with thermal diffusion
      IF(IGR.GT.1) THEN
        DO 10001 J1=1,JSP
          JLC (J1) = L0F(JSOL(J1))
10001   CONTINUE
        L0TEM    = L0F(JTEM)
      ENDIF
      IXL=IABS(IXL)
      IF(IGR.EQ.13) GO TO 13
      IF(IGR.EQ.19) GO TO 19
      GO TO (1,2,3,4,5,6,25,8,9,10,11,12,13,14,25,25,25,25,19,20,25,
     125,23,24),IGR
   25 CONTINUE
      RETURN
C*****************************************************************
C
C--- GROUP 1. Run title and other preliminaries
C
    1 GO TO (1001,1002,1003),ISC
 2061 FORMAT('(CVD) Soret model is ACTIVE.')
 1001 CONTINUE
 
C Default values for diffusion, properties, chemistry, solvers
      MCDOPT = 2
      SMRLX  = 0.5
      BINOPT = 4
      THMDIF = .FALSE.
      MCPROP = 3
      CHMRLX = 0.5
      NGREAC = 0
      NSREAC = 0
      ANJAC  = .TRUE.
      GSPNAM=' '; CTEMP=' '
C Calls to GETSPD for SPEDAT settings in Q1
      CALL GETSPD('CVD','MCDOPT',2,RDUM,MCDOPT,.FALSE.,' ',IERR)
      IF(MCDOPT.EQ.3) THEN
        CALL GETSPD('CVD','SMRLX',1,SMRLX,IDUM,.FALSE.,' ',IERR)
      ENDIF
      CALL GETSPD('CVD','BINOPT',2,RDUM,BINOPT,.FALSE.,' ',IERR)
      CALL GETSPD('CVD','THMDIF',3,RDUM,IDUM,THMDIF,' ',IERR)
      IF(THMDIF) THEN
        THMOPT = 1
        THMFRQ = 1
        THMRLX = 1.0
        CALL GETSPD('CVD','THMOPT',2,RDUM,THMOPT,.FALSE.,' ',IERR)
        CALL GETSPD('CVD','THMFRQ',2,RDUM,THMFRQ,.FALSE.,' ',IERR)
        CALL GETSPD('CVD','THMRLX',1,THMRLX,IDUM,.FALSE.,' ',IERR)
      ENDIF
      CALL GETSPD('CVD','NGREAC',2,RDUM,NGREAC,.FALSE.,' ',IERR)
      IF(NGREAC.GT.0) THEN
        RNAME(1:6) = 'GREAC('
        DO 998 I = 1,NGREAC
          IF(I.LT.10) THEN
            WRITE(RNAME(7:7),'(I1)') I
            RNAME(8:15) = ')       '
          ELSE
            WRITE(RNAME(7:8),'(I2)') I
            RNAME(9:15) = ')      '
          ENDIF
          CALL GETSPD('CVD',RNAME,2,RDUM,GREAC(I),.FALSE.,' ',IERR)
998     CONTINUE
      ENDIF
      CALL GETSPD('CVD','NSREAC',2,RDUM,NSREAC,.FALSE.,' ',IERR)
      IF(NSREAC.GT.0) THEN
        RNAME(1:6) = 'SREAC('
        DO 999 I = 1,NSREAC
          IF(I.LT.10) THEN
            WRITE(RNAME(7:7),'(I1)') I
            RNAME(8:15) = ')       '
          ELSE
            WRITE(RNAME(7:8),'(I2)') I
            RNAME(9:15) = ')      '
          ENDIF
          CALL GETSPD('CVD',RNAME,2,RDUM,SREAC(I),.FALSE.,' ',IERR)
999     CONTINUE
      ENDIF
      CALL GETSPD('CVD','CHMRLX',1,CHMRLX,IDUM,.FALSE.,' ',IERR)
      CALL GETSPD('CVD','MCPROP',2,RDUM,MCPROP,.FALSE.,' ',IERR)
      CALL GETSPD('CVD','ANJAC',3,RDUM,IDUM,ANJAC,' ',IERR)
      IF(MCPROP.EQ.1.OR.BINOPT.EQ.2) THEN
        IF(EQZ(TMP1A)) THEN
          CALL WRIT40('You need to set TMP1A for options chosen')
          CALL WRYT40('You need to set TMP1A for options chosen')
          CALL SET_ERR(461,
     *          'Error. You need to set TMP1A for options chosen.',1)
C          STOP
        ENDIF
      ENDIF
 
 
C --------------------------------------------------------------------
C Construction of:  JSOL(j)       = Contiguous array of specie indicies.
C                   JSOLV(indvar) = Inverse array of JSOL(j).
C                   JSP           = Total number of species.
C                   JST           = Index of STOREd-only specie.
C                   JFST          = Index of first SOLVEd specie.
C                   JLST          = Index of last SOLVEd specie.
C --------------------------------------------------------------------
      JTOT = 0
      DO 10043 JVAR = C1, C1+JMAXS-1
        IF(STORE(JVAR).OR.SOLVE(JVAR)) THEN
          IF(NAME(JVAR)(1:1)/='S') CYCLE ! species name starts 'S'
          READ(NAME(JVAR)(2:4),'(I3)',IOSTAT=IOS) JD
          IF(IOS/=0) CYCLE ! species name must be Snnn where nnn is a number
          JTOT = JTOT + 1
          JDUM(JTOT) = JVAR
        ENDIF
10043 CONTINUE
      DO 10042 J1 = 1, JTOT
        READ(NAME(JDUM(J1))(2:4),'(I3)') JD
        JNUM(J1) = JD
10042 CONTINUE
      GSPED = 'SPECIDAT'
      CALL GETSPC(GSPED,JTOT,JNUM,JPHSE,GSPNAM)
      DO 10044 J1 = 1, JTOT
        JPHSV(JDUM(J1)) = JPHSE(J1)
10044 CONTINUE
      JSU =    0
      JSP =    0
      JST = -999
!C#### JCL 14.08.18 increase number of gas species to 100 (CHAM-J 08/012/18)
      DO 1004 JVAR = C1, C1+JMAXS-1
!      DO 1004 JVAR = C1, MIN(NPHI,C1+JMAXS-1)
        IF(SOLVE(JVAR).AND.JPHSV(JVAR).EQ.'G ') THEN
          JSP = JSP + 1
          JGAS(JSP) = JVAR
        ELSEIF(SOLVE(JVAR).AND.JPHSV(JVAR).EQ.'A ') THEN
          CALL WRIT40('ERROR: You cannot solve adsorbed species')
          CALL WRYT40('ERROR: You cannot solve adsorbed species')
          CALL SET_ERR(462,
     *          'Error. You cannot solve adsorbed species.',1)
C          STOP
        ELSEIF(STORE(JVAR).AND.JPHSV(JVAR).EQ.'G ') THEN
          JSP = JSP + 1
          JGAS(JSP) = JVAR
          IF(JST.EQ.-999) THEN
            JST = JSP
          ELSE
            CALL WRIT40('Only allowed one STOREd gas species     ')
            CALL WRYT40('Only allowed one STOREd gas species     ')
            CALL SET_ERR(463,
     *          'Error. Only allowed one STOREd gas species.',1)
C            STOP
          ENDIF
        ELSEIF(STORE(JVAR).AND.JPHSV(JVAR).EQ.'A ') THEN
          JSU = JSU + 1
          JADS(JSU) = JVAR
        ENDIF
1004  CONTINUE
 
      DO 10051 J1 = 1, JSP
10051   JSOL(J1) = JGAS(J1)
      DO 1006 J1 = 1, JSU
       J = JSP + J1
1006   JSOL(J) = JADS(J1)
      IF(JST.NE.-999) THEN
C Reorder species indices and names
        JFST              = 2
        IF(JST.NE.1) JFST = JST
        CDUM              = GSPNAM(1)
        JJ                = JSOL(JST)
        JK                = JSOL(1  )
        JSOL(1  )         = JJ
        JSOL(JST)         = JK
        GSPNAM(1)         = GSPNAM(JST)
        GSPNAM(JST)       = CDUM
        JST               = 1
      ELSE
        CALL WRIT40('You must STORE one gas species          ')
        CALL WRYT40('You must STORE one gas species          ')
        CALL SET_ERR(464,
     *          'Error. You must STORE one gas species.',1)
      ENDIF
C
      IF(NSREAC==0.AND.JSU>0) THEN
        CALL WRIT40('Adsorbed species present but no Surface Reactions')
        CALL WRYT40('Adsorbed species present but no Surface Reactions')
        CALL SET_ERR(428,'Error. See result file',1)
      ENDIF
C
C... Ensure that all 'A' species come after 'G' species
      III=0
      DO II=1,JTOT
        IF(JPHSE(II)=='G') THEN
          III=III+1; IORD(III)=II
        ENDIF
      ENDDO
      DO II=1,JTOT
        IF(JPHSE(II)=='A') THEN
          III=III+1; IORD(III)=II
        ENDIF
      ENDDO
      DO II=1,JTOT
        III=IORD(II)
        CTEMP(II)=GSPNAM(III)
      ENDDO
      GSPNAM=CTEMP
 
      DO 2012 J2      = 1 , JIND
 2012 JSOLV(J2)       = -999
      DO 2011 J1      = 1 , JTOT
 2011 JSOLV(JSOL(J1)) = J1
      JLST                  = -999
      IF(JSP.GT.1) THEN
        DO 2014 J1            = 2 , JSP
        JVAR                  = JSOL(J1)
 2014   IF(JVAR.GT.JLST) JLST = JVAR
        JLST                  = JSOLV(JLST)
      ENDIF
 
C Set identifying integers for non-generic variable STOREs.
      JTEM  = LBNAME('TEM1')
      JPRPS = LBNAME('PRPS')
      JSPHT = LBNAME('SPHT')
      JCOND = LBNAME('CNDT')
      JT0   = LBNAME('T0')
      JNE   = LBNAME('NE')
      JION  = LBNAME('GION')
      JVPOR = LBNAME('VPOR')
      JDEPO = LBNAME('DEPO')
      JEPOR = LBNAME('EPOR')
      JNPOR = LBNAME('NPOR')
      JHPOR = LBNAME('HPOR')
      IF(JTEM.EQ.0) THEN
        CALL WRYT40('No storage provided for TEM1            ')
        CALL WRIT40('No storage provided for TEM1            ')
        CALL SET_ERR(465,
     *          'Error. No storage provided for TEM1.',1)
C        STOP
      ENDIF
      IF(JDEPO.EQ.0) THEN
        CALL WRYT40('No storage provided for DEPO            ')
        CALL WRIT40('No storage provided for DEPO            ')
        CALL SET_ERR(466,
     *          'Error. No storage provided for DEPO.',1)
C        STOP
      ENDIF
C  Call GXPLAS for Plasma variables
      IF(JT0.GT.0) CALL GXPLAS
C --------------------------------------------------------------------
C Make user part of f-array:
C --------------------------------------------------------------------
C     L0MW   = store for species molar mass, (kg/mol).
C     L0MF   = store for molar fraction of species, (non-dim).
C     L0BD   = store for species binary diffusion coefficients when
C              BINOPT is equal to 2.
C     L0SG   = store for standard Gibbs energy of formation, (J/mol).
C     L0EK   = store for Lennard Jones Parameters, (K). ??check
C     L0HO   = store for standard heat of formation, (J/mol).
C     L0SO   = store for standard entropy, (J/mol.K).
C     L0SP   = store for species conservation check, (kg).
C     L0AC   = store for species average field mass concentration.
C  The next 10 stores all changed to include face area
C     L0JE   = stores for old EAST  face diffusion fluxes, (kg/m2/s).
C     L0JN   = stores for old NORTH face diffusion fluxes, (kg/m2/s).
C     L0JH   = stores for old HIGH  face diffusion fluxes, (kg/m2/s).
C     L0TE   = stores for new EAST  face diffusion fluxes, (kg/m2/s).
C     L0TN   = stores for new NORTH face diffusion fluxes, (kg/m2/s).
C     L0TH   = stores for new HIGH  face diffusion fluxes, (kg/m2/s).
C     L0UE   = store  for EAST  face diffusion velocity, (kg/m2.s).
C     L0UN   = store  for NORTH face diffusion velocity, (kg/m2.s).
C     L0UH   = store  for HIGH  face diffusion velocity, (kg/m2.s).
C     L0UL   = store  for LOW   face diffusion velocity, (kg/m2.s).
C     L0DH   = store for multicomponent diffusivity of HIGH slab.
C     L0SC   = store for diffusion source (SU).
C     L0ARL  = store for cell LOW face area, (m2).
C     L0TR   = store for thermal diffusion coefficient.
C     L0ST   = store for thermal diffusion source.
C     L0JTE  = EAST thermal diffusion flux (kg/s)
C     L0JTN  = NORTH thermal diffusion flux (kg/s)
C     L0JTH  = HIGH thermal diffusion flux (kg/s)
C     L0CSP3 = store for positive chemical reaction rates in field
C     L0CSN3 = store for negative chemical reaction rates in field
C     L0CSPT = store for positive reaction sources (gas and surface)
C     L0CSNT = store for negative reaction sources (gas and surface)
C     L0DIF  = store for species diffusion coefficients in field
C     L0GC   = store for mass fractions of species
C     L0SUK  = source term (numerator)
C     L0APK  = source term (denominator)
C     L0RF   = forward reaction rate
C     L0RR   = reverse reaction rate
C --------------------------------------------------------------------
      CALL GXMAKE(L0GC,JSP,'MAFR')
      CALL GXMAKE(L0DIF,JSP*NXNY*NZ,'WDIF')
      DO 2088 J=1,JSP*NXNY*NZ
        F(L0DIF+J) = 1.0e-10
2088  CONTINUE
      CALL GXMAKE(L0CSP3,JSP*NXNY*NZ,'CS_P')
      CALL GXMAKE(L0CSN3,JSP*NXNY*NZ,'CS_N')
      CALL GXMAKE(L0CSPT,JSP*NXNY*NZ,'CSTP')
      CALL GXMAKE(L0CSNT,JSP*NXNY*NZ,'CSTN')
      DO 2089 J=1,JSP*NXNY*NZ
        F(L0CSP3+J) = 1.0e-10
        F(L0CSN3+J) = -1.0e-10
        F(L0CSPT+J) = 1.0e-10
        F(L0CSNT+J) = -1.0e-10
2089  CONTINUE
      IF(JSU.GT.0) THEN
        CALL GXMAKE(L0MW,JTOT+1,'MOLW')
      ELSE
        CALL GXMAKE(L0MW,JSP,'MOLW')
      ENDIF
      CALL GXMAKE(L0MF,JSP,'MOLF')
      IF(BINOPT.EQ.2)CALL GXMAKE(L0BD,JSP*JSP,'BIND')
      CALL GXMAKE(L0SG,JSP,'STGB')
      CALL GXMAKE(L0EK,JSP,'EPSK')
      CALL GXMAKE(L0HO,JSP,'STHF')
      CALL GXMAKE(L0SO,JSP,'STEN')
      IF(USOURC) THEN
        CALL GXMAKE(L0ARL,NXNY,'AREL')
        DO 2049 J   = L0ARL+1 , L0ARL+NXNY
 2049   F(J)        = 0.0
        IF(THMDIF) THEN
          IF(.NOT.SOLVE(JTEM)) THEN
            CALL WRYT40('SORET effect is active but TEM1 is not ')
            CALL WRIT40('SORET effect is active but TEM1 is not ')
          ENDIF
          WRITE(LUPR1,2061)
          CALL GXMAKE(L0TR,(JSP-1)*NXNY*NZ,'THDI')
          DO 2059 J = L0TR+1 , L0TR+(JSP-1)*NXNY*NZ
 2059     F(J)      = 0.0
          CALL GXMAKE(L0ST,(JSP-1)*NXNY*NZ,'SORT')
          DO 2076 J = L0ST+1 , L0ST+(JSP-1)*NXNY*NZ
 2076     F(J)      = 0.0
          CALL GXMAKE(L0JTE,JSP*NXNY*NZ,'JTE ')
          CALL GXMAKE(L0JTN,JSP*NXNY*NZ,'JTN ')
          CALL GXMAKE(L0JTH,JSP*NXNY*NZ,'JTH ')
          DO 20591 J=1,JSP*NXNY*NZ
            F(L0JTE+J) = 0.0
            F(L0JTN+J) = 0.0
            F(L0JTH+J) = 0.0
20591     CONTINUE
        ENDIF
          IF(NX.GT.1) THEN
            CALL GXMAKE(L0JE,JSP*NXNY*NZ,'JFLE')
            CALL GXMAKE(L0TE,JSP*NXNY   ,'TFLE')
            DO 1992 J   = L0JE+1 , L0JE+JSP*NXNY*NZ
 1992       F(J)        = 0.0
            DO 2022 J   = L0TE+1 , L0TE+JSP*NXNY
 2022       F(J)        = 0.0
            CALL GXMAKE(L0UE,NXNY,'UDFE')
            DO 2040 J = L0UE+1 , L0UE+NXNY
 2040       F(J)      = 0.0
          ENDIF
          IF(NY.GT.1) THEN
            CALL GXMAKE(L0JN,JSP*NXNY*NZ,'JFLN')
            CALL GXMAKE(L0TN,JSP*NXNY   ,'TFLN')
            DO 1993 J   = L0JN+1 , L0JN+JSP*NXNY*NZ
 1993       F(J)        = 0.0
            DO 2023 J   = L0TN+1 , L0TN+JSP*NXNY
 2023       F(J)        = 0.0
            CALL GXMAKE(L0UN,NXNY,'UDFN')
            DO 2041 J = L0UN+1 , L0UN+NXNY
 2041       F(J)      = 0.0
          ENDIF
          IF(NZ.GT.1) THEN
            CALL GXMAKE(L0JH,JSP*NXNY*NZ,'JFLH')
            CALL GXMAKE(L0TH,JSP*NXNY   ,'TFLH')
            DO 1994 J = L0JH+1 , L0JH+JSP*NXNY*NZ
 1994       F(J)      = 0.0
            DO 2024 J = L0TH+1 , L0TH+JSP*NXNY
 2024       F(J)      = 0.0
            CALL GXMAKE(L0UH,(JSP-1)*NXNY,'VELH')
            CALL GXMAKE(L0UL,(JSP-1)*NXNY,'VELL')
            DO 2025 J = L0UH+1 , L0UH+(JSP-1)*NXNY
 2025       F(J)      = 0.0
            DO 2026 J = L0UL+1 , L0UL+(JSP-1)*NXNY
 2026       F(J)      = 0.0
            CALL GXMAKE(L0DH,(JSP-1)*NXNY,'MDHI')
            DO 2027 J = L0DH+1 , L0DH+(JSP-1)*NXNY
 2027       F(J)      = 0.0
          ENDIF
      ENDIF
      IF(JSU.GT.0) THEN
        CALL GXMAKE(L0FSU,JSP*NXNY*NZ,'FSUR')
        DO 2101 J = L0FSU+1, L0FSU+JSP*NXNY*NZ
           F(J)   =0.0
 2101   CONTINUE
      ENDIF
      IF(JSU.GT.0) THEN
        CALL GXMAKE(L0SOCC,JSU,'SOCC')
        DO 2102 J = L0SOCC+1,L0SOCC+JSU
          F(J) = 0.0
 2102   CONTINUE
      ENDIF
C
      LINIL= .TRUE.
      ISHWR = 0
      ISURF = 0
      DO 1016 IPN=1,NUMPAT
        IF(NAMPAT(IPN)(1:4).EQ.'SHWR') ISHWR = ISHWR + 1
        IF(NAMPAT(IPN)(1:4).EQ.'SHWX') ISHWR = ISHWR + 1
        IF(NAMPAT(IPN)(1:4).EQ.'SURF') ISURF = ISURF + 1
1016  CONTINUE
      IF(ISHWR.GT.0) THEN
        CALL GXMAKE(L0SHWR,3*ISHWR,'SHWR')
        ISH = 0
      ENDIF
      IF(ISURF.GT.0) THEN
        CALL GXMAKE(L0SURF,2*ISURF,'SURF')
        ISR = 0
      ENDIF
      IF(ISHWR.GT.0.OR.ISURF.GT.0) THEN
        DO 1017 IPN=1,NUMPAT
          IF(NAMPAT(IPN)(1:4).EQ.'SHWR'.OR.NAMPAT(IPN)(1:4).EQ.
     1                           'SHWX') THEN
            ISH = ISH + 1
            F(L0SHWR+ISH*3-2) = REAL(IPN)
            L0PTCH = I0PAT + 10*(IPN-1)
            ITYPE = F(L0PTCH+2) + 0.1
            IF(NAMPAT(IPN)(1:4).EQ.'SHWR') THEN
              F(L0PTCH+5)=ABS(F(L0PTCH+5))
              IF(ITYPE.EQ.2.OR.ITYPE.EQ.3) THEN
                F(L0PTCH+3) = F(L0PTCH+3) - 1
                F(L0PTCH+4) = F(L0PTCH+4) + 1
              ELSEIF(ITYPE.EQ.4.OR.ITYPE.EQ.5) THEN
                F(L0PTCH+5) = F(L0PTCH+5) - 1
                F(L0PTCH+6) = F(L0PTCH+6) + 1
              ELSEIF(ITYPE.EQ.6.OR.ITYPE.EQ.7) THEN
                F(L0PTCH+7) = F(L0PTCH+7) - 1
                F(L0PTCH+8) = F(L0PTCH+8) + 1
              ENDIF
              IF(F(L0PTCH+3).LE.0.0.OR.F(L0PTCH+4).GT.REAL(NX).
     1          OR.F(L0PTCH+5).LE.0.0.OR.F(L0PTCH+6).GT.REAL(NY).
     2          OR.F(L0PTCH+7).LE.0.0.OR.F(L0PTCH+8).GT.REAL(NZ)) THEN
                CALL WRYT40('Shower patch too near edge of domain    ')
                CALL WRIT40('Shower patch too near edge of domain    ')
                CALL SET_ERR(467,
     *               'Error. Shower patch too near edge of domain.',1)
C                STOP
              ENDIF
              F(L0PTCH+5)=-(F(L0PTCH+5))
            ENDIF
          ENDIF
          IF(NAMPAT(IPN)(1:4).EQ.'SURF') THEN
            ISR = ISR + 1
            F(L0SURF+ISR*2-1) = REAL(IPN)
          ENDIF
1017    CONTINUE
      ENDIF
C Change SHWR patch name so that RESULT file echoes new name and size
      DO 1018 IPN = 1, NUMPAT
        IF(NAMPAT(IPN)(1:4).EQ.'SHWR') THEN
          PNAME       = NAMPAT(IPN)
          PNAME(1:4)  = 'SHWX'
          NAMPAT(IPN) = PNAME
        ENDIF
1018  CONTINUE
C --------------------------------------------------------------------
C Pick up data from databases:
C --------------------------------------------------------------------
      GCHED = 'CHEMIDAT'
      IF(NGREAC.GT.0) THEN
        ERROR=.FALSE.
        JSP0=JSP
        CALL GETCDB(3,'GAS',LUPR1,4,GCHED,GSPNAM,JMAXS,JSP,
     1              MIN(JHOMAX,NGREAC),NGREAC,GREAC,IGRT,IGAMMA,JGAMMA,
     2              GRPAR,90,GAPAR,GFTHB,GALFA,GSALFA,' ',0,KGSP,
     3              JDEP,0,0,0,0,0.0,0.0,0.0,NIWS,IWS,NRWS,RWS,NCWS,CWS,
     4              ERROR,GGNMS,GSRCN)
        IF(ERROR) THEN
          CALL WRIT40(' ERROR REPORTED BY GETCDB FOR GAS-PHASE ')
          CALL WRIT40(' REACTIONS                              ')
          CALL WRYT40(' ERROR REPORTED BY GETCDB FOR GAS-PHASE ')
          CALL WRYT40(' REACTIONS                              ')
      CALL SET_ERR(310,
     *   'ERROR REPORTED BY GETCDB FOR GAS-PHASE REACTIONS.',1)
C          CALL EAROUT(1)
        ELSEIF(JSP.NE.JSP0) THEN
          JSP=JSP0
        ENDIF
      ENDIF
      IF(NSREAC.GT.0) THEN
        ERROR=.FALSE.
        JTT=JTOT
        JSP0=JSP
        CALL WRIT40(' SURFACE CHEMISTRY READ                 ')
        CALL GETCDB(3,'SURFACE',LUPR1,4,GCHED,GSPNAM,JMAXS,JTT,
     1              MIN(JSUMAX,NSREAC),NSREAC,SREAC,ISRT,
     2              IGAMMA,JGAMMA,SRPAR,90,SAPAR,SFTHB,GALFA,GSALFA,
     3              GHSOL,KDEP,KGSP,JDEP,0,0,0,1,GDSOL,GMSOL,GFAC,NIWS,
     4              IWS,NRWS,RWS,NCWS,CWS,ERROR,GSNMS,GSRCN)
        IF(ERROR) THEN
          CALL WRIT40(' ERROR REPORTED BY GETCDB FOR SURFACE-  ')
          CALL WRIT40(' PHASE REACTIONS                        ')
          CALL WRYT40(' ERROR REPORTED BY GETCDB FOR SURFACE-  ')
          CALL WRYT40(' PHASE REACTIONS                        ')
      CALL SET_ERR(311,
     *   'ERROR REPORTED BY GETCDB FOR SURFACE-PHASE REACTIONS.',1)
C          CALL EAROUT(1)
        ELSEIF(JSP.NE.JSP0) THEN
          JSP=JSP0
        ENDIF
C If using complex chemistry, check that reactions are compatible
        LMODEL = .FALSE.
        LMODEL2 = .TRUE.
C Arora & Pollard model
        DO 3020 J = 1, NSREAC
          IF(ISRT(J).EQ.401) THEN
            LMODEL = .TRUE.
          ELSE
            LMODEL2 = .FALSE.
          ENDIF
 3020   CONTINUE
        IF(LMODEL.AND..NOT.LMODEL2)
     1        CALL WRIT40(' Error in combination of surface react. ')
        LMODEL = .FALSE.
        LMODEL2 = .TRUE.
C Wang model
        DO 3023 J = 1, NSREAC
          IF(ISRT(J).EQ.402) THEN
            LMODEL = .TRUE.
          ELSE
            LMODEL2 = .FALSE.
          ENDIF
 3023   CONTINUE
        IF(LMODEL.AND..NOT.LMODEL2)
     1        CALL WRIT40(' Error in combination of surface react. ')
        KVAC = 0
C Identify the vacancy species in the Arora and Wang model
        DO 3021 J = JSP+1, JSP+JSU
          IF(GSPNAM(J).EQ.'1V(A)') KVAC = J
 3021   CONTINUE
      ENDIF
      GTRAD = 'TRANSDAT'
      GTHED = 'THERMDAT'
      CALL WRIT40('TRANSPORT DATABASE WILL BE READ         ')
      CALL GETTRA(GTRAD,GSPNAM,JSP,JSP,JNDATA,GDATA)
      CALL WRIT40('THERMODYNAMIC DATABASE WILL BE READ     ')
      IF(JSU.EQ.0) THEN
        CALL GETTHM(GTHED,GSPNAM,JSP,JSP,JCOEFF,GCOMP,
     1             GEL,GMEL,F(L0MW+1),GTLOW,GTHIGH,GTAV,GLOW,GHIGH)
      ELSE
        CALL GETTHM(GTHED,GSPNAM,JTOT+1,JTOT+1,JCOEFF,GCOMP,
     1             GEL,GMEL,F(L0MW+1),GTLOW,GTHIGH,GTAV,GLOW,GHIGH)
      ENDIF
C Open file for detailed surface chemistry output
      IF(JSU.GT.0) THEN
        OPEN(UNIT=3,FILE='srfdat',ACCESS='SEQUENTIAL',
     1      FORM='FORMATTED',STATUS='UNKNOWN')
        REWIND(3)
        WRITE(3,'(''SURFACE CHEMISTRY RESULTS:'')')
      ENDIF
      DO 10181 J = 1, JSU
        CDUM = GSPNAM(JSP+J)
        IF(CDUM(1:1).EQ.'1') THEN
          GSOCC(J) = 1
        ELSEIF(CDUM(1:1).EQ.'2') THEN
          GSOCC(J) = 2
        ELSEIF(CDUM(1:1).EQ.'4') THEN
          GSOCC(J) = 4
        ELSE
          GSOCC(J) = 1
        ENDIF
        F(L0SOCC+J) = GSOCC(J)
10181  CONTINUE
C All low temperatures are set to 100.
      DO 10041 I=1,JTOT+1
        GTLOW(I)=100.0
10041 CONTINUE
      DO 1011 J1   = 1 , JSP
        F(L0SG+J1) = GDATA(3,J1)
        F(L0EK+J1) = GDATA(2,J1)
        F(L0HO+J1) = GDATA(4,J1)
 1011 F(L0SO+J1)   = GDATA(5,J1)
C --------------------------------------------------------------------
C Calculation of Binary Diffusion Coefficients, for option 2:
C
C Two options - 2: Empirical constant, reference temperature.
C   (BINOPT)    4: Empirical constant, actual    temperature.
C --------------------------------------------------------------------
      IF(BINOPT.EQ.2) THEN
        DO 1009 J1                  = 1 , JSP
          DO 1009 J2                = 1 , JSP
            JJ                      = J2+JSP*(J1-1)
            JK                      = J1+JSP*(J2-1)
            IF(J1.EQ.J2) F(L0BD+JJ) = 0.0
            IF(J1.GT.J2) F(L0BD+JJ) = F(L0BD+JK)
 1009   IF(J1.LT.J2) F(L0BD+JJ) = FBDIF(F(L0MW+J1),F(L0SG+J1),
     1          F(L0EK+J1),F(L0MW+J2),F(L0SG+J2),F(L0EK+J2),TMP1A,
     1          TMP1A,PRESS0,BINOPT)
      ENDIF
C --------------------------------------------------------------------
C Temperatures for use in calculating the in-cell and neighbour
C binary diffusion coefficients using reference temperature option,
C in Group 8 Section 12.
C --------------------------------------------------------------------
      GTEMPP = TMP1A
      GTEMPE = TMP1A
      GTEMPW = TMP1A
      GTEMPN = TMP1A
      GTEMPS = TMP1A
      GTEMPH = TMP1A
      GTEMPL = TMP1A
C --------------------------------------------------------------------
C
C     User may here change message transmitted to the VDU screen
      IF(IGR.EQ.1.AND.ISC.EQ.1.AND..NOT.NULLPR) THEN
        CALL WRYT40('-------------(PHOENICS-CVD)-------------')
        CALL WRYT40('GXCVD  file is gxcvd      of:  09/05/02 ')
        CALL WRYT40('Associated files are:  cvdchm           ')
        CALL WRYT40('                       cvddat           ')
        CALL WRYT40('                       cvdmth           ')
        CALL WRYT40('                       cvdpr1           ')
        CALL WRYT40('                       cvdpr2           ')
        CALL WRYT40('                       cvdsur           ')
        CALL WRYT40('                       cvdutl           ')
        CALL WRYT40('                       gxs2sr           ')
        CALL WRYT40('                       gxplas           ')
        CALL WRYT40('                       gxview           ')
        CALL WRYT40('                       view             ')
        CALL WRYT40('Data files are:        SPECIDAT         ')
        CALL WRYT40('                       TRANSDAT         ')
        CALL WRYT40('                       THERMDAT         ')
        CALL WRYT40('                       CHEMIDAT         ')
        CALL WRYT40('                       OPTICDAT         ')
        CALL WRYT40('----------------------------------------')
      ENDIF
C.................................................................
      CALL WRIT40(' CHEMISTRY DATA WRITTEN TO FILE         ')
      CALL GOUTP(JSP,GSPNAM,F(L0EK+1),F(L0SG+1),F(L0MW+1),GTLOW,GTHIGH,
     1           GTAV,GLOW,GHIGH,JCAR,JOPMU,JOPLAB,JOPCP,JOPDIF,
     1           GREFT,GTMIN,GTMAX,GERMU,GERLAB,GERDIF,PRESS0,IGAMMA,
     1           JCOEFF,JGAMMA,JMAXS,JHOMAX,NSREAC,GHSOL,GMSOL,GDSOL,
     1           GALFA,GSALFA,JSUMAX,JDEP,NGREAC,FPRT,IGRT,
     1           GRPAR,90,GAPAR,GFTHB,JTOT,JSU,GSRCN,GSOCC)
C
      IF(USOLVE) THEN
        ISTIFF=1
        DO 3100 K = 1, JSP
          K1 = JSOL(K)
          ISLN(K1)=546
3100    CONTINUE
        CALL GXMAKE(L0SUK,NXNY*JSP,'SUKT')
        CALL GXMAKE(L0APK,NXNY*JSP,'APKT')
        CALL GXMAKE(L0RF,NXNY*NGREAC,'RF  ')
        CALL GXMAKE(L0RR,NXNY*NGREAC,'RR  ')
        DO 3110 J = 1, NXNY*JSP
          F(L0SUK+J) = 0.0
          F(L0APK+J) = 1.0E+10
3110    CONTINUE
C  Determine machine accuracy for newton solver:
        CALL MCHACC(GABSOL,GRELAT)
      ELSE
        ISTIFF=0
      ENDIF
C
      RETURN
 1002 CONTINUE
C  Call GXPLAS for Plasma variables
      IF(JT0.GT.0) CALL GXPLAS
      RETURN
 1003 CONTINUE
      RETURN
C*****************************************************************
C
C--- GROUP 2. Transience; time-step specification
C
    2 CONTINUE
      RETURN
C*****************************************************************
C
C--- GROUP 3. X-direction grid specification
C
    3 CONTINUE
      RETURN
C*****************************************************************
C
C--- GROUP 4. Y-direction grid specification
C
    4 CONTINUE
      RETURN
C*****************************************************************
C
C--- GROUP 5. Z-direction grid specification
C
    5 CONTINUE
      RETURN
C*****************************************************************
C
C--- GROUP 6. Body-fitted coordinates or grid distortion
C
    6 CONTINUE
      RETURN
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.LE.GRND--- phase 1 additional velocity. Index VELAD
      RETURN
   82 CONTINUE
C   * ------------------- SECTION  2 ---------------------------
C    For U2AD.LE.GRND--- phase 2 additional velocity. Index VELAD
      RETURN
   83 CONTINUE
C   * ------------------- SECTION  3 ---------------------------
C    For V1AD.LE.GRND--- phase 1 additional velocity. Index VELAD
      RETURN
   84 CONTINUE
C   * ------------------- SECTION  4 ---------------------------
C    For V2AD.LE.GRND--- phase 2 additional velocity. Index VELAD
      RETURN
   85 CONTINUE
C   * ------------------- SECTION  5 ---------------------------
C    For W1AD.LE.GRND--- phase 1 additional velocity. Index VELAD
      RETURN
   86 CONTINUE
C   * ------------------- SECTION  6 ---------------------------
C    For W2AD.LE.GRND--- phase 2 additional velocity. Index VELAD
      RETURN
   87 CONTINUE
C   * ------------------- SECTION 7 ---- Volumetric source for gala
      RETURN
   88 CONTINUE
C   * ------------------- SECTION 8 ---- Convection fluxes
      RETURN
   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
      RETURN
  810 CONTINUE
C   * ------------------- SECTION 10 --- Convection neighbours
      RETURN
  811 CONTINUE
C   * ------------------- SECTION 11 --- Diffusion neighbours
      IF(INDVAR.GE.C1.AND.INDVAR.LE.C1+JMAXS-1.AND.INDVAR/=JTEM) THEN
        L0VPON            = 0
        L0PRPN            = 0
        JZL               = IZ-1
        JZH               = IZ+1
        IF(IZ.EQ.1 ) JZL  = IZ
        IF(IZ.EQ.NZ) JZH  = IZ
        IF(NZ.GT.1) THEN
          L0DENH          = L0F(ANYZ(DEN1,JZH)    )
          L0TEMH          = L0F(ANYZ(JTEM,JZH)    )
          L0TEML          = L0F(ANYZ(JTEM,JZL)    )
        ENDIF
        IF(NDIREC.EQ.1) THEN
          LDIFF = LAN
          LNEI  = LD8
          L0T   = L0TN + JT1
        ELSEIF(NDIREC.EQ.3) THEN
          LDIFF = LAE
          LNEI  = LD8
          L0T   = L0TE + JT1
        ELSEIF(NDIREC.EQ.5) THEN
          LDIFF = LAH
C          LDIFF = LD11
          LNEI  = LD7
          L0T   = L0TH + JT1
        ENDIF
        CALL FN10(-L0T,INDVAR,LNEI,0.0,1.0,-1.0)
        CALL FN26(-L0T,LDIFF)
C  Stefan-Maxwell coding starts
        IF(GRNDNO(8,-ABS(PRNDTL(JSOL(JFST)))).AND.MCDOPT.EQ.3) THEN
          IF(NDIREC.EQ.1) THEN
            CALL SUB3(L0DXX,L0F(LXYDY),L0DXG,L0F(LXYDYG),
     1                L0AREA,L0F(ANORTH))
            CALL SUB3(L0U,L0UN,L0T,L0TN,L0J,L0JN)
            CALL SUB2(L0DENN,L0DEN,L0TEMN,L0TEM)
            IF(JPRPS.NE.0) L0PRPN = L0PRP
            IF(JVPOR.NE.0) L0VPON = L0VPOR
          ELSEIF(NDIREC.EQ.3) THEN
            CALL SUB3(L0DXX,L0F(LXYDX),L0DXG,L0F(LXYDXG),
     1                L0AREA,L0F(AEAST))
            CALL SUB3(L0U,L0UE,L0T,L0TE,L0J,L0JE)
            CALL SUB2(L0DENN,L0DEN,L0TEMN,L0TEM)
            IF(JPRPS.NE.0) L0PRPN = L0PRP
            IF(JVPOR.NE.0) L0VPON = L0VPOR
          ELSEIF(NDIREC.EQ.5) THEN
            CALL SUB3(L0DXX,L0F(LXYDZ),L0DXG,L0F(LXYDZG),
     1                L0AREA,L0F(AHIGH))
            CALL SUB3(L0U,L0UH+JT2,L0T,L0TH,L0J,L0JH)
            CALL SUB2(L0DENN,L0DENH,L0TEMN,L0TEMH)
            IF(JPRPS.NE.0) L0PRPN = L0PRPH
            IF(JVPOR.NE.0) L0VPON = L0VPOH
          ENDIF
          DO 8120 JX=1,NX
          DO 8120 JY=1,NY
            JC = JY + NY*(JX-1)
            JCZ = (IZ-1)*NXNY + JC
            IF(NDIREC.EQ.1) THEN
              IF(JY.EQ.NY) GO TO 8120
              JCN = JC + 1
              JCZN = JCZ + 1
            ELSEIF(NDIREC.EQ.3) THEN
              IF(JX.EQ.NX.AND..NOT.XCYCLE) GO TO 8120
              JCN = JC + NY
              IF(JX.EQ.NX) JCN = JY
              JCZN = (IZ-1)*NXNY + JCN
            ELSEIF(NDIREC.EQ.5) THEN
              IF(IZ.EQ.NZ) GO TO 8120
              JCN = JC
              JCZN = JCZ + NXNY
            ENDIF
            GAREA = F(L0AREA+JC)
            FACEBL = .FALSE.
            IF(GAREA.LT.TINY) FACEBL=.TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRP+JC).GT.99.0) FACEBL=.TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRPN+JCN).GT.99.0) FACEBL=.TRUE.
            IF(JVPOR.NE.0.AND.F(L0VPOR+JC).LT.TINY) FACEBL=.TRUE.
            IF(JVPOR.NE.0.AND.F(L0VPON+JCN).LT.TINY) FACEBL=.TRUE.
            IF(FACEBL) THEN
              IF(NDIREC.EQ.5) THEN
                IF(IZ.EQ.1) THEN
                  F(L0UL+JC+JT2) = 0.0
                ELSE
                  F(L0UL+JC+JT2) = F(L0U+JC)
                ENDIF
              ENDIF
              F(L0U+JC) = 0.0
              F(L0T+JC+JT1) = 0.0
              GO TO 8120
            ENDIF
            IF(BINOPT.EQ.4) THEN
              GTEMPN = F(L0TEMN+JCN) + TEMP0
              GTEMPP = F(L0TEM+JC) + TEMP0
            ENDIF
            GDENP = F(L0DEN+JC)
            GDENN = F(L0DENN+JCN)
            GSUM  = 0.0
            GMMWP = 0.0
            GMMWN = 0.0
            DO 8123 J1=1,JSP
              IF(J1.NE.JSUB) THEN
                JKJ1 = (J1-1)*NXNY*NZ
                IF(BINOPT.EQ.4) THEN
                  GTEMP = (GTEMPP+GTEMPN)/2.0
                  GBIN = FBDIF(F(L0MW+JSUB),F(L0SG+JSUB),
     1                          F(L0EK+JSUB),F(L0MW+J1),F(L0SG+J1),
     2                          F(L0EK+J1),TMP1A,GTEMP,PRESS0,BINOPT)
                ELSEIF(BINOPT.EQ.2) THEN
                  GBIN = F(L0BD+JSUB+JSP*(J1-1))
                ENDIF
                GSUM = GSUM + F(L0J+JCZ+JKJ1)/((F(L0MW+J1)+TINY)*GBIN)
              ENDIF
              GMMWP = GMMWP + F(JLC(J1)+JC)/(F(L0MW+J1)+TINY)
              IF(NDIREC.EQ.5) THEN
                GMMWN = GMMWN + F(JLCH(J1)+JCN)/(F(L0MW+J1)+TINY)
              ELSE
                GMMWN = GMMWN + F(JLC(J1)+JCN)/(F(L0MW+J1)+TINY)
              ENDIF
8123        CONTINUE
            GMMWP = 1.0/(GMMWP+TINY)
            GMMWN = 1.0/(GMMWN+TINY)
            GMMW = (GMMWP+GMMWN)/2.0
            DX = F(L0DXX+JC)
            DXG = F(L0DXG+JC)
            IF(NDIREC.EQ.5) THEN
              DXNEI = 2.*F(L0DXG+JC)-F(L0DXX+JC)
            ELSE
              DXNEI = F(L0DXX+JCN)
            ENDIF
            DIFP = F(L0DIF+(JSUB-1)*NXNY*NZ+JCZ)
            DIFN = F(L0DIF+(JSUB-1)*NXNY*NZ+JCZN)
            GDIF=(DX+DXNEI)/(DX/DIFP+DXNEI/DIFN)
            GDFDEN =
     1                 (DX+DXNEI)/(DX/DIFP/GDENP+DXNEI/DIFN/GDENN)
            GTERM1 = -GDFDEN*(GMMWN-GMMWP)/GMMW/DXG*GAREA
C Second term of velocity calculation.
            GTERM2 = -GDIF*GMMW*GSUM
            IF(NDIREC.EQ.5) THEN
              IF(IZ.EQ.1) THEN
                F(L0UL+JC+JT2) = 0.0
              ELSE
                F(L0UL+JC+JT2) = F(L0U+JC)
              ENDIF
            ENDIF
            F(L0U+JC) = GTERM1 - GTERM2
            F(L0T+JC+JT1) = F(L0T+JC+JT1)+
     1               AMAX1(0.0,F(L0U+JC)*F(JLC(JSUB)+JC))
            IF(NDIREC.EQ.5) THEN
              F(L0T+JC+JT1) = F(L0T+JC+JT1)
     1                -AMAX1(0.0,-F(L0U+JC)*F(JLCH(JSUB)+JCN))
            ELSE
              F(L0T+JC+JT1) = F(L0T+JC+JT1)
     1                -AMAX1(0.0,-F(L0U+JC)*F(JLC(JSUB)+JCN))
            ENDIF
8120      CONTINUE
        ENDIF
C  Stefan-Maxwell coding ends
C  Soret coding begins
        IF(THMDIF) THEN
          IF(JCHECK.NE.0.AND.ISWEEP.GT.THMFRQ) GO TO 2047
          IF(NDIREC.EQ.1) THEN
            CALL SUB3(L0DXG,L0F(LXYDYG),L0AREA,L0F(ANORTH),L0JT,L0JTN)
            L0TEMN = L0TEM
            IF(JPRPS.NE.0) L0PRPN = L0PRP
            IF(JVPOR.NE.0) L0VPON = L0VPOR
          ELSEIF(NDIREC.EQ.3) THEN
            CALL SUB3(L0DXG,L0F(LXYDXG),L0AREA,L0F(AEAST),L0JT,L0JTE)
            L0TEMN = L0TEM
            IF(JPRPS.NE.0) L0PRPN = L0PRP
            IF(JVPOR.NE.0) L0VPON = L0VPOR
          ELSEIF(NDIREC.EQ.5) THEN
            CALL SUB3(L0DXG,L0F(LXYDZG),L0AREA,L0F(AHIGH),L0JT,L0JTH)
            L0TEMN = L0TEMH
            IF(JPRPS.NE.0) L0PRPN = L0PRPH
            IF(JVPOR.NE.0) L0VPON = L0VPOH
          ENDIF
          DO 8110 JX=1,NX
          DO 8110 JY=1,NY
            JC = JY + NY*(JX-1)
            JCZ = (IZ-1)*NXNY + JC
            IF(NDIREC.EQ.1) THEN
              IF(JY.EQ.NY) GO TO 8110
              JCN = JC + 1
              JCZN = JCZ + 1
            ELSEIF(NDIREC.EQ.3) THEN
              IF(JX.EQ.NX.AND..NOT.XCYCLE) GO TO 8110
              JCN = JC + NY
              IF(JX.EQ.NX) JCN = JY
              JCZN = (IZ-1)*NXNY + JCN
            ELSEIF(NDIREC.EQ.5) THEN
              IF(IZ.EQ.NZ) GO TO 8110
              JCN = JC
              JCZN = JCZ + NXNY
            ENDIF
            GAREA = F(L0AREA+JC)
            FACEBL = .FALSE.
            IF(GAREA.LT.TINY) FACEBL=.TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRP+JC).GT.99.0) FACEBL=.TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRPN+JCN).GT.99.0) FACEBL=.TRUE.
            IF(JVPOR.NE.0.AND.F(L0VPOR+JC).LT.TINY) FACEBL=.TRUE.
            IF(JVPOR.NE.0.AND.F(L0VPON+JCN).LT.TINY) FACEBL=.TRUE.
            IF(FACEBL) THEN
              F(L0JT+JK1+JCZ) = 0.0
              GO TO 8110
            ENDIF
            GTEMPN = F(L0TEMN+JCN) + TEMP0
            GTEMPP = F(L0TEM+JC) + TEMP0
            GTMEAN = 0.5*(GTEMPP+GTEMPN)
            DIFP = F(L0TR+JK2+JCZ)
            DIFN = F(L0TR+JK2+JCZN)
            DIFAV = 0.5*(DIFP+DIFN)
            DXG = F(L0DXG+JC)
            F(L0JT+JK1+JCZ) = DIFAV*(GTEMPP-GTEMPN)/DXG/GTMEAN*GAREA
            IF(JSUB.EQ.JLST) THEN
              GSUM = 0.0
              DO 81101 J1=2,JSP
                GSUM = GSUM + F(L0JT+JK1+JCZ)
81101         CONTINUE
              F(L0JT+JCZ) = -GSUM
            ENDIF
8110      CONTINUE
        ENDIF
C  Soret coding ends
2047  ENDIF
      RETURN
  812 CONTINUE
C   * ------------------- SECTION 12 --- Linearised sources
C --------------------------------------------------------------------
C Only access this section when the Stefan-Maxwell option for
C multicomponent diffusion or the Soret effect has been selected.
C --------------------------------------------------------------------
      IF(USOURC.AND.ISWEEP.GE.FSWEEP) THEN
        IF(INDVAR.GE.C1.AND.INDVAR.LE.C1+JMAXS-1.AND.INDVAR/=JTEM) THEN
          IF(WHF(INDVAR)) THEN
            L0SU = L0F(L3SU)
            L0AP = L0F(L3AP)
          ELSE
            L0SU = L0F(LSU )
            L0AP = L0F(LAP )
          ENDIF
          IF(JSP.EQ.1) GO TO 10104
          IF(GRNDNO(8,-ABS(PRNDTL(JSOL(JFST)))).AND.MCDOPT.EQ.3) THEN
C --------------------------------------------------------------------
C Source term for the Stefan-Maxwell Equation.
C --------------------------------------------------------------------
            L0L          = L0F(LAMPR)
            DO 8121 JX   = 1 , NX
              DO 8121 JY = 1 , NY
                JC       = JY+NY*(JX-1)
                JCZ      = JC+NXNY*(IZ-1)
                JCO      = JC
                IF(WHF(INDVAR))JCO = JCZ
C --------------------------------------------------------------------
C Source Su:
C --------------------------------------------------------------------
                JCN = JC + 1
                JCS = JC - 1
                JCE = JC + NY
                JCW = JC - NY
                IF(JX.EQ.1.AND.XCYCLE) JCW = (NX-1)*NY + JC
                IF(JX.EQ.NX.AND.XCYCLE) JCE = JY
                GSOURC = 0.0
                IF(XCYCLE.OR.(JX.LT.NX)) GSOURC = GSOURC-
     1            (AMAX1(0.0,F(L0UE+JC)*F(JLC(JSUB)+JC))-
     2            AMAX1(0.0,-F(L0UE+JC)*F(JLC(JSUB)+JCE)))
                IF(XCYCLE.OR.(JX.GT.1 )) GSOURC = GSOURC-
     1            (AMAX1(0.0,-F(L0UE+JCW)*F(JLC(JSUB)+JC))-
     2            AMAX1(0.0,F(L0UE+JCW)*F(JLC(JSUB)+JCW)))
                IF(JY.LT.NY) GSOURC = GSOURC-
     1            (AMAX1(0.0,F(L0UN+JC)*F(JLC(JSUB)+JC))-
     2            AMAX1(0.0,-F(L0UN+JC)*F(JLC(JSUB)+JCN)))
                IF(JY.GT.1 ) GSOURC = GSOURC-
     1            (AMAX1(0.0,-F(L0UN+JCS)*F(JLC(JSUB)+JC))-
     2            AMAX1(0.0,F(L0UN+JCS)*F(JLC(JSUB)+JCS)))
                IF(IZ.LT.NZ) GSOURC = GSOURC-
     1             (AMAX1(0.0,F(L0UH+JC+JT2)*F(JLC(JSUB)+JC))-
     2             AMAX1(0.0,-F(L0UH+JC+JT2)*F(JLCH(JSUB)+JC)))
                IF(IZ.GT.1 ) GSOURC = GSOURC-
     1             (AMAX1(0.0,-F(L0UL+JC+JT2)*F(JLC(JSUB)+JC))-
     2             AMAX1(0.0,F(L0UL+JC+JT2)*F(JLCL(JSUB)+JC)))
C Add zero to the coefficient Ap if source is negative.
                GCOEFF = AMAX1(0.0,-GSOURC/(F(JLC(JSUB)+JC)+TINY))
C --------------------------------------------------------------------
C Coefficient Ap:
C --------------------------------------------------------------------
                GADD        = GSOURC
                F(L0SU+JCO) = F(L0SU+JCO)+GADD
                F(L0AP+JCO)  = F(L0AP+JCO)+ABS(GADD)/
     1                        (F(JLC(JSUB)+JC)+TINY)
 8121       CONTINUE
C --------------------------------------------------------------------
          ENDIF
C --------------------------------------------------------------------
C Transfer species diffusion fluxes from temporary (2d) storage
C to (3d) storage at end of the last specie counter.
C --------------------------------------------------------------------
            DO 81271 JC            = 1 , NXNY
              JCZ                  = JC+NXNY*(IZ-1)
              IF(NX.GT.1) THEN
                  IF(GRNDNO(8,-ABS(PRNDTL(JSOL(JFST)))).AND.
     1                                           MCDOPT.EQ.3) THEN
                    F(L0JE+JCZ+JK1) = SMRLX*F(L0TE+JC+JT1)+
     1                                  (1.-SMRLX)*F(L0JE+JCZ+JK1)
                  ELSE
                    F(L0JE+JCZ+JK1) = F(L0TE+JC+JT1)
                  ENDIF
              ENDIF
              IF(NY.GT.1) THEN
                  IF(GRNDNO(8,-ABS(PRNDTL(JSOL(JFST)))).AND.
     1                                           MCDOPT.EQ.3) THEN
                    F(L0JN+JCZ+JK1) = SMRLX*F(L0TN+JC+JT1)+
     1                                  (1.-SMRLX)*F(L0JN+JCZ+JK1)
                  ELSE
                    F(L0JN+JCZ+JK1) = F(L0TN+JC+JT1)
                  ENDIF
              ENDIF
              IF(NZ.GT.1) THEN
                  IF(GRNDNO(8,-ABS(PRNDTL(JSOL(JFST)))).AND.
     1                                           MCDOPT.EQ.3) THEN
                    F(L0JH+JCZ+JK1) = SMRLX*F(L0TH+JC+JT1)+
     1                                  (1.-SMRLX)*F(L0JH+JCZ+JK1)
                  ELSE
                    F(L0JH+JCZ+JK1) = F(L0TH+JC+JT1)
                  ENDIF
              ENDIF
81271       CONTINUE
          IF(JSUB.EQ.JLST) THEN
            DO 8127 JC             = 1 , NXNY
              JCZ                  = JC+NXNY*(IZ-1)
              IF(NX.GT.1) THEN
                GFLUXE             = 0.0
                DO 8128 J1         = 2 , JSP
                  JKJ1             = (J1-1)*NXNY*NZ
 8128           GFLUXE             = GFLUXE+F(L0JE+JCZ+JKJ1)
                F(L0JE+JCZ)        = -GFLUXE
              ENDIF
              IF(NY.GT.1) THEN
                GFLUXN             = 0.0
                DO 8129 J1         = 2 , JSP
                  JKJ1             = (J1-1)*NXNY*NZ
 8129           GFLUXN             = GFLUXN+F(L0JN+JCZ+JKJ1)
                F(L0JN+JCZ)        = -GFLUXN
              ENDIF
              IF(NZ.GT.1) THEN
                GFLUXH             = 0.0
                DO 8130 J1         = 2 , JSP
                  JKJ1             = (J1-1)*NXNY*NZ
 8130           GFLUXH             = GFLUXH+F(L0JH+JCZ+JKJ1)
                F(L0JH+JCZ)        = -GFLUXH
              ENDIF
 8127       CONTINUE
          ENDIF
10104     CONTINUE
C --------------------------------------------------------------------
C Source term for the Soret Effect:
C --------------------------------------------------------------------
          IF(THMDIF) THEN
            DO 2053   JX = 1 , NX
            DO 2053 JY = 1 , NY
              JC       = JY+NY  *(JX-1)
              BLCKGE = .FALSE.
              IF(JPRPS.NE.0.AND.F(L0PRP+JC).GE.100) BLCKGE = .TRUE.
              IF(JVPOR.NE.0.AND.F(L0VPOR+JC).LE.TINY) BLCKGE = .TRUE.
              IF(BLCKGE) GO TO 2053
              JCZ      = JC+NXNY*(IZ-1)
              JCO      = JC
              IF(WHF(INDVAR))JCO = JCZ
              IF(JCHECK.NE.0.AND.ISWEEP.GT.THMFRQ) GO TO 20471
              JCZS = JCZ - 1
              JCZW = JCZ - NY
              JCZL = JCZ - NXNY
              IF(JX.EQ.1.AND.XCYCLE) JCZW = (IZ-1)*NXNY+(NX-1)*NY+JC
C --------------------------------------------------------------------
C Source Su:
C --------------------------------------------------------------------
              GSOURT          = 0.0
              IF(XCYCLE.OR.(JX.LT.NX)) GSOURT = GSOURT -
     1                                          F(L0JTE+JK1+JCZ)
              IF(XCYCLE.OR.(JX.GT.1)) GSOURT = GSOURT +
     1                                          F(L0JTE+JK1+JCZW)
              IF(JY.LT.NY) GSOURT = GSOURT - F(L0JTN+JK1+JCZ)
              IF(JY.GT.1) GSOURT = GSOURT + F(L0JTN+JK1+JCZS)
              IF(IZ.LT.NZ) GSOURT = GSOURT - F(L0JTH+JK1+JCZ)
              IF(IZ.GT.1) GSOURT = GSOURT + F(L0JTH+JK1+JCZL)
              F(L0ST+JCZ+JK2) = GSOURT
20471         IF(THMFRQ.GT.1) THEN
                GADD = F(L0ST+JCZ+JK2)
              ELSE
                GADD = GSOURT
              ENDIF
              F(L0SU+JCO) = F(L0SU+JCO) + GADD
 2053       CONTINUE
C --------------------------------------------------------------------
C Dumping total flux due to thermal diffusion.
C --------------------------------------------------------------------
          ENDIF
        ENDIF
C  Energy transport by multicomponent diffusion
        IF(INDVAR.EQ.JTEM) THEN
          L0HPOR            = 0
          L0HPOL            = 0
          L0NPOR            = 0
          L0EPOR            = 0
          IF(NX.GT.1.AND.JEPOR.NE.0) L0EPOR = L0F(JEPOR)
          IF(NY.GT.1.AND.JNPOR.NE.0) L0NPOR = L0F(JNPOR)
          JZL               = IZ-1
          JZH               = IZ+1
          IF(IZ.EQ.1 ) JZL  = IZ
          IF(IZ.EQ.NZ) JZH  = IZ
          IF(NZ.GT.1) THEN
            L0TEMH          = L0F(ANYZ(JTEM,JZH)    )
            L0TEML          = L0F(ANYZ(JTEM,JZL)    )
            IF(JHPOR.NE.0) THEN
              L0HPOR        = L0F(JHPOR             )
              L0HPOL        = L0F(ANYZ(JHPOR,JZL)   )
            ENDIF
          ENDIF
          IF(WHF(INDVAR)) THEN
            L0SU = L0F(L3SU)+(IZ-1)*NXNY
            L0AP = L0F(L3AP)+(IZ-1)*NXNY
          ELSE
            L0SU = L0F(LSU )
            L0AP = L0F(LAP )
          ENDIF
          DO 81272 JY=1,NY
          DO 81272 JX=1,NX
            JC = (JX-1)*NY+JY
            JCZ = JC+NXNY*(IZ-1)
            IF(JVPOR.NE.0.AND.F(L0VPOR+JC).LT.TINY) GO TO 81272
            IF(JPRPS.NE.0.AND.F(L0PRP+JC).GT.99) GO TO 81272
            GSOURC = 0.0
            GAP = 0.0
            DO 81273 J1=1,JSP
              JKJ1 = (J1-1)*NXNY*NZ
              IF(NX.GT.1) THEN
                JCW = JC-NY
                JCE = JC+NY
                IF(XCYCLE.AND.JX.EQ.1) JCW = JY+NY*(NX-1)
                IF(XCYCLE.AND.JX.EQ.NX) JCE = JY
                JCZW = JCW+NXNY*(IZ-1)
                IF(XCYCLE.OR.(JX.LT.NX)) THEN
                  FACEBL = .FALSE.
                  IF(JVPOR.NE.0.AND.F(L0VPOR+JCE).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(JPRPS.NE.0.AND.F(L0PRP+JCE).GT.99)
     1                                  FACEBL = .TRUE.
                  IF(JEPOR.NE.0.AND.F(L0EPOR+JC).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(FACEBL) GO TO 81274
                  FLUX = F(L0JE+JCZ+JKJ1)
                  IF(THMDIF) FLUX = FLUX + F(L0JTE+JCZ+JKJ1)
                  IF(FLUX.GE.-TINY) THEN
                    TEMP = F(L0TEM+JC) + TEMP0
                    CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                    CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                    CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC - FLUX*CPT
                    GAP = GAP + FLUX*CPT/TEMP
                  ELSE
                    TEMP = F(L0TEM+JCE) + TEMP0
                    CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                    CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                    CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC - FLUX*CPT
                  ENDIF
                ENDIF
81274           IF(XCYCLE.OR.(JX.GT.1)) THEN
                  FACEBL = .FALSE.
                  IF(JVPOR.NE.0.AND.F(L0VPOR+JCW).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(JPRPS.NE.0.AND.F(L0PRP+JCW).GT.99)
     1                                  FACEBL = .TRUE.
                  IF(JEPOR.NE.0.AND.F(L0EPOR+JCW).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(FACEBL) GO TO 81275
                  FLUX = F(L0JE+JCZW+JKJ1)
                  IF(THMDIF) FLUX = FLUX + F(L0JTE+JCZW+JKJ1)
                  IF(FLUX.GT.TINY) THEN
                    TEMP = F(L0TEM+JCW) + TEMP0
                      CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                   GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                      CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                   GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                      CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC + FLUX*CPT
                  ELSE
                    TEMP = F(L0TEM+JC) + TEMP0
                      CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                   GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                      CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                   GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                      CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC + FLUX*CPT
                    GAP = GAP - FLUX*CPT/TEMP
                  ENDIF
                ENDIF
              ENDIF
81275         IF(NY.GT.1) THEN
                JCS = JC-1
                JCN = JC+1
                JCZS = JCS+NXNY*(IZ-1)
                IF(JY.LT.NY) THEN
                  FACEBL = .FALSE.
                  IF(JVPOR.NE.0.AND.F(L0VPOR+JCN).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(JPRPS.NE.0.AND.F(L0PRP+JCN).GT.99)
     1                                  FACEBL = .TRUE.
                  IF(JEPOR.NE.0.AND.F(L0NPOR+JC).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(FACEBL) GO TO 81276
                  FLUX = F(L0JN+JCZ+JKJ1)
                  IF(THMDIF) FLUX = FLUX + F(L0JTN+JCZ+JKJ1)
                  IF(FLUX.GE.-TINY) THEN
                    TEMP = F(L0TEM+JC) + TEMP0
                    CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                    CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                    CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC - FLUX*CPT
                    GAP = GAP + FLUX*CPT/TEMP
                  ELSE
                    TEMP = F(L0TEM+JCN) + TEMP0
                    CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                    CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                    CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC - FLUX*CPT
                  ENDIF
                ENDIF
81276           IF(JY.GT.1) THEN
                  FACEBL = .FALSE.
                  IF(JVPOR.NE.0.AND.F(L0VPOR+JCS).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(JPRPS.NE.0.AND.F(L0PRP+JCS).GT.99)
     1                                  FACEBL = .TRUE.
                  IF(JEPOR.NE.0.AND.F(L0NPOR+JCS).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(FACEBL) GO TO 81277
                  FLUX = F(L0JN+JCZS+JKJ1)
                  IF(THMDIF) FLUX = FLUX + F(L0JTN+JCZS+JKJ1)
                  IF(FLUX.GE.TINY) THEN
                    TEMP = F(L0TEM+JCS) + TEMP0
                    CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                    CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                    CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC + FLUX*CPT
                  ELSE
                    TEMP = F(L0TEM+JC) + TEMP0
                    CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                    CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                    CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC + FLUX*CPT
                    GAP = GAP - FLUX*CPT/TEMP
                  ENDIF
                ENDIF
              ENDIF
81277         IF(NZ.GT.1) THEN
                JCZL = JCZ - NXNY
                JCZH = JCZ + NXNY
                IF(IZ.LT.NZ) THEN
                  FACEBL = .FALSE.
                  IF(JVPOR.NE.0.AND.F(L0VPOH+JC).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(JPRPS.NE.0.AND.F(L0PRPH+JC).GT.99)
     1                                  FACEBL = .TRUE.
                  IF(JEPOR.NE.0.AND.F(L0HPOR+JC).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(FACEBL) GO TO 81278
                  FLUX = F(L0JH+JCZ+JKJ1)
                  IF(THMDIF) FLUX = FLUX + F(L0JTH+JCZ+JKJ1)
                  IF(FLUX.GE.-TINY) THEN
                    TEMP = F(L0TEM+JC) + TEMP0
                    CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                    CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                    CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC - FLUX*CPT
                    GAP = GAP + FLUX*CPT/TEMP
                  ELSE
                    TEMP = F(L0TEMH+JC) + TEMP0
                    CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                    CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                    CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC - FLUX*CPT
                  ENDIF
                ENDIF
81278           IF(IZ.GT.1) THEN
                  FACEBL = .FALSE.
                  IF(JVPOR.NE.0.AND.F(L0VPOL+JC).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(JPRPS.NE.0.AND.F(L0PRPL+JC).GT.99)
     1                                  FACEBL = .TRUE.
                  IF(JEPOR.NE.0.AND.F(L0HPOL+JC).LT.TINY)
     1                                  FACEBL = .TRUE.
                  IF(FACEBL) GO TO 81273
                  FLUX = F(L0JH+JCZL+JKJ1)
                  IF(THMDIF) FLUX = FLUX + F(L0JTH+JCZL+JKJ1)
                  IF(FLUX.GE.TINY) THEN
                    TEMP = F(L0TEML+JC) + TEMP0
                    CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                    CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                    CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC + FLUX*CPT
                  ELSE
                    TEMP = F(L0TEM+JC) + TEMP0
                    CPT = H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),TEMP)
                    CPT = CPT - H0(JCOEFF,GLOW(1,J1),GHIGH(1,J1),
     1                 GTLOW(J1),GTAV(J1),GTHIGH(J1),250.0)
                    CPT = CPT/(F(L0MW+J1)+TINY)
                    GSOURC = GSOURC + FLUX*CPT
                    GAP = GAP - FLUX*CPT/TEMP
                  ENDIF
                ENDIF
              ENDIF
81273       CONTINUE
            F(L0SU+JC) = F(L0SU+JC) + GSOURC
            F(L0AP+JC) = F(L0AP+JC) + GAP
81272     CONTINUE
        ENDIF
      ENDIF
C --------------------------------------------------------------------
      RETURN
  813 CONTINUE
C   * ------------------- SECTION 13 --- Correction coefficients
      RETURN
  814 CONTINUE
C   * ------------------- SECTION 14 --- User's own solver
      USER=.FALSE.
      IF(ISTIFF.EQ.1) THEN
        IF(INDVAR.GE.JSOL(1).AND.INDVAR.LE.JSOL(JLST)) THEN
          USER=.TRUE.
          CALL FN1(LD7,0.0)
        ENDIF
      ENDIF
      RETURN
  815 CONTINUE
C   * ------------------- SECTION 15 --- Change solution
      RETURN
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 GO TO (91,92,93,94,95,96,97,98,99,900,901,902,903,904,905),ISC
C*****************************************************************
  900 CONTINUE
C   * ------------------- SECTION 10 ---------------------------
C    For TMP1.LE.GRND--------- phase-1 temperature Index TEMP1
      RETURN
  901 CONTINUE
C   * ------------------- SECTION 11 ---------------------------
C    For TMP2.LE.GRND--------- phase-2 temperature Index TEMP2
      RETURN
  902 CONTINUE
C   * ------------------- SECTION 12 ---------------------------
C    For EL1.LE.GRND--------- phase-1 length scale Index LEN1
      RETURN
  903 CONTINUE
C   * ------------------- SECTION 13 ---------------------------
C    For EL2.LE.GRND--------- phase-2 length scale Index LEN2
      RETURN
  904 CONTINUE
C   * ------------------- SECTION 14 ---------------------------
C    For SOLVE(TEMP1)-------------------------- phase-1 specic heat
C --------------------------------------------------------------------
      RETURN
  905 CONTINUE
C   * ------------------- SECTION 15 ---------------------------
C    For SOLVE(TEMP2)-------------------------- phase-2 specic heat
      RETURN
   91 CONTINUE
C   * ------------------- SECTION  1 ---------------------------
C    For RHO1.LE.GRND--- density for phase 1       Index DEN1
      RETURN
   92 CONTINUE
C   * ------------------- SECTION  2 ---------------------------
C    For DRH1DP.LE.GRND--- D(LN(DEN))/DP for phase 1
C                                                  Index D1DP
      RETURN
   93 CONTINUE
C   * ------------------- SECTION  3 ---------------------------
C    For RHO2.LE.GRND--- density for phase 2       Index DEN2
      RETURN
   94 CONTINUE
C   * ------------------- SECTION  4 ---------------------------
C    For DRH2DP.LE.GRND--- D(LN(DEN))/DP for phase 2
C                                                  Index D2DP
      RETURN
   95 CONTINUE
C   * ------------------- SECTION  5 ---------------------------
C    For ENUT.LE.GRND--- reference turbulent kinematic viscosity
C                                                  Index VIST
      RETURN
   96 CONTINUE
C   * ------------------- SECTION  6 ---------------------------
C    For ENUL.LE.GRND--- reference laminar kinematic viscosity
C                                                  Index VISL
      RETURN
   97 CONTINUE
C   * ------------------- SECTION  7 ---------------------------
C    For PRNDTL( ).LE.GRND--- laminar PRANDTL nos., or diffusivity
C                                                  Index LAMPR
C  Call GXPLAS for Plasma variables
      IF(JT0.GT.0) CALL GXPLAS
C --------------------------------------------------------------------
C Dumping of thermal conductivity field to optional CNDT STORE.
C --------------------------------------------------------------------
      IF(ISWEEP.EQ.LSWEEP.OR.SWPRIN) THEN
        IF(INDVAR.EQ.JTEM) THEN
crj          IF((GRNDNO(10,-ABS(PRNDTL(JTEM)))).AND.STORE(JCOND)) THEN
          IF(JPRPS.NE.0.AND.STORE(JCOND)) THEN
            CALL FN2(JCOND,LAMPR,0.0,-1.0)
          ENDIF
C --------------------------------------------------------------------
C Dumping species multicomponent diffusivity fields to optional
C MD02 - MD20 STOREs.
C --------------------------------------------------------------------
        ELSEIF(INDVAR.GE.C1.AND.INDVAR.LE.C1+JMAXS-1.AND.INDVAR/=JTEM)
     1                                                              THEN
          JNAME(1:2)                = 'MD'
          WRITE(JNAME(3:4),'(I2)') JSUB
          IF(JSUB.LT.10) JNAME(3:3) = '0'
          JDIFF                     = LBNAME(JNAME)
          IF(STORE(JDIFF)) THEN
            CALL FN0(JDIFF,LAMPR)
          ENDIF
          IF(JSUB.EQ.2) THEN
            JNAME(3:4) = '01'
            JDIFF = LBNAME(JNAME)
            IF(STORE(JDIFF)) CALL FN0(JDIFF,-(L0DIF+(IZ-1)*NX*NY))
          ENDIF
        ENDIF
      ENDIF
C --------------------------------------------------------------------
C Calculation of thermal diffusion coefficient for use in Soret effect.
C --------------------------------------------------------------------
      IF(USOURC.AND.THMDIF.AND.JSUB.GE.1.AND.JSUB.LE.JSP) THEN
C Minimum molar concentration set for GFMIN
        GFMIN        = 1.0E-06
C IGSH used to put thermal diffusion coefficient in correct place in
C L0TR storage (note that L0TEM and JLC already take account of IGSH)
        DO 2062 JC   = 1 , NXNY
          JCZ        = JC+NXNY*(IZ-1+IGSH)
          GSUM       = 0.0
          DO 2063 J1 = 1 , JSP
 2063     GSUM       = GSUM+F(JLC(J1)+JC)/(F(L0MW+J1)+TINY)
          IF(GSUM.LT.GFMIN) GO TO 2062
          BLCKGE = .FALSE.
          IF(IGSH.EQ.0) THEN
            IF(JVPOR.NE.0.AND.F(L0VPOR+JC).LE.TINY) BLCKGE = .TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRP+JC).GE.100) BLCKGE = .TRUE.
          ELSEIF(IGSH.EQ.1) THEN
            IF(JVPOR.NE.0.AND.F(L0VPOH+JC).LE.TINY) BLCKGE = .TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRPH+JC).GE.100) BLCKGE = .TRUE.
          ENDIF
          IF(BLCKGE) GO TO 2062
          DO 2064 J1 = 1 , JSP
 2064     F(L0MF+J1) = F(JLC(J1)+JC)/(F(L0MW+J1)+TINY)/(GSUM+TINY)
          IF(THMOPT.EQ.1.OR.THMOPT.EQ.2) THEN
C THMOPT = 1 for Clark-Jones formulation using rigid spheres assumption.
C THMOPT = 2 for Clark-Jones formulation using Lennard-Jones potential.
            F(L0TR+JCZ+JK2) = FMDFT1(JSUB,F(L0MF+1),JSP,JMAXS,JST,
     1         F(L0MW+1),F(L0SG+1),F(L0EK+1),F(L0TEM+JC)+TEMP0
     2         ,GFMIN,THMOPT)
          ELSEIF(THMOPT.EQ.3.OR.THMOPT.EQ.4) THEN
C THMOPT = 3 for exact formulation using rigid spheres assumption.
C THMOPT = 4 for exact formulation using Lennard-Jones potential.
            CALL DETLOW(F(L0MF+1),JSP,JMAXS,F(L0MW+1),F(L0EK+1),
     1        F(L0SG+1),F(L0TEM+JC)+TEMP0,GFMIN,GLOWER,GDUMMY
     2        ,GDET,THMOPT-2)
            F(L0TR+JCZ+JK2) = FMDFT2(JSUB,F(L0MF+1),GFMIN,JSP,JMAXS,
     1      F(L0MW+JSUB),GLOWER,GUPPER,GDET)
          ELSE
            CALL WRYT40('THMOPT incorrectly set                  ')
            CALL WRIT40('THMOPT incorrectly set                  ')
                CALL SET_ERR(468, 'Error. THMOPT incorrectly set.',1)
C            STOP
          ENDIF
 2062   CONTINUE
C --------------------------------------------------------------------
C Dumping species thermal diffusivity fields to optional
C TD02 - TD20 STOREs.
C --------------------------------------------------------------------
        IF(ISWEEP.EQ.LSWEEP) THEN
          JNAME(1:2)                = 'TD'
          WRITE(JNAME(3:4),'(I2)') JSUB
          IF(JSUB.LT.10) JNAME(3:3) = '0'
          JDIFF                     = LBNAME(JNAME)
          IF(STORE(JDIFF)) THEN
            L0DUMP       = L0F(JDIFF)
            DO 2065 JC   = 1 , NXNY
              JCZ        = JC+NXNY*(IZ-1)
 2065       F(L0DUMP+JC) = F(L0TR+JCZ+JK2)
          ENDIF
        ENDIF
      ENDIF
C --------------------------------------------------------------------
      RETURN
   98 CONTINUE
C   * ------------------- SECTION  8 ---------------------------
C    For PHINT( ).LE.GRND--- interface value of first phase
C                                                  Index FII1
      RETURN
   99 CONTINUE
C   * ------------------- SECTION  9 ---------------------------
C    For PHINT( ).LE.GRND--- interface value of second phase
C                                                  Index FII2
      RETURN
C*****************************************************************
C
C--- GROUP 10. Inter-phase-transfer processes and properties
C
   10 GO TO (101,102,103,104),ISC
  101 CONTINUE
C   * ------------------- SECTION  1 ---------------------------
C    For CFIPS.LE.GRND--- inter-phase friction coeff.
C                                                  Index INTFRC
      RETURN
  102 CONTINUE
C   * ------------------- SECTION  2 ---------------------------
C    For CMDOT.EQ.GRND- inter-phase mass transfer  Index INTMDT
      RETURN
  103 CONTINUE
C   * ------------------- SECTION  3 ---------------------------
C    For CINT( ).EQ.GRND--- phase1-to-interface transfer coefficients
C                                                  Index COI1
      RETURN
  104 CONTINUE
C   * ------------------- SECTION  4 ---------------------------
C    For CINT( ).EQ.GRND--- phase2-to-interface transfer coefficients
C                                                  Index COI2
      RETURN
C*****************************************************************
C
C--- GROUP 11. Initialization of variable or porosity fields
C                                                  Index VAL
   11 CONTINUE
      RETURN
C*****************************************************************
C
C--- GROUP 12. Convection and diffusion adjustments
C
   12 CONTINUE
      RETURN
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
C                                       Index for Coefficient - CO
C                                       Index for Value       - VAL
   13 CONTINUE
      GO TO (130,131,132,133,134,135,136,137,138,139,1310,
     11311,1312,1313,1314,1315,1316,1317,1318,1319,1320,1321),ISC
  130 CONTINUE
C------------------- SECTION  1 ------------- coefficient = GRND
      RETURN
  131 CONTINUE
C------------------- SECTION  2 ------------- coefficient = GRND1
C/////////////////////////////////////////////////////////////////
C Gas-phase chemistry
C/////////////////////////////////////////////////////////////////
      IF(NPATCH(1:4).EQ.'CHEM') THEN
        L0CO = L0F(CO)
        IF(INDVAR.EQ.JTEM) THEN
          DO 1301 IXY = 1, NX*NY
            F(L0CO+IXY) = TINY
1301      CONTINUE
        ELSEIF(INDVAR.GE.C1.AND.INDVAR.LE.(C1+JMAXS-1).AND.INDVAR/=JTEM)
     1                                                              THEN
          L0C1 = L0F(INDVAR)
          DO 1302 IXY = 1, NX*NY
            BLCKGE = .FALSE.
            IF(JVPOR.NE.0.AND.F(L0VPOR+IXY).LE.TINY) BLCKGE = .TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRP+IXY).GE.100) BLCKGE = .TRUE.
            IF(BLCKGE) THEN
              F(L0CO+IXY) = 0.0
            ELSE
              GC1  = F(L0C1  +IXY)
              GSP1 = F(L0CSP3+JK1+(IZ-1)*NXNY+IXY)
              GSN1 = F(L0CSN3+JK1+(IZ-1)*NXNY+IXY)
C The linearised source term has the form
C    S = SP + (SN/c)*phi
C which is written as
C    S = Coef*(Val-phi).
C Thus
C    Val = -c*SP/SN ; Coef = - SN/c
C  No linearisation if stiff solver (Jacobian used elsewhere)
              IF(ISTIFF.EQ.0) THEN
                F(L0CO+IXY) = -(GSN1-TINY)/(GC1+TINY)
              ELSE
                F(L0CO+IXY) = FIXFLU
              ENDIF
            ENDIF
1302      CONTINUE
        ENDIF
C/////////////////////////////////////////////////////////////////
C Surface-phase chemistry is not linearised: mass loss makes it
C more complicated than gas-phase, so FIXFLU is used
C/////////////////////////////////////////////////////////////////
C     Automated false time step underrelaxation:
      ELSEIF(NPATCH.EQ.'RELT') THEN
        IF(INDVAR.GE.C1.AND.INDVAR.LE.(C1+JMAXS-1).AND.INDVAR/=JTEM)
     1                                                             THEN
          L0CO = L0F(CO)
          L0C1 = L0F(INDVAR)
          DO 1303 IXY = 1, NX*NY
            BLCKGE = .FALSE.
            IF(JVPOR.NE.0.AND.F(L0VPOR+IXY).LE.TINY) BLCKGE = .TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRP+IXY).GE.100) BLCKGE = .TRUE.
            IF(BLCKGE) THEN
              F(L0CO+IXY) = 0.0
            ELSE
C An optimal false time step is estimated, based on
C   T = rho*c/Source
C for both the positive and negative sources.
C Of these two estimates, the smaller is chosen.
              GFTSP =
     1          ABS(F(L0DEN+IXY)*F(L0C1+IXY)/
     1          (F(L0CSPT+JK1+(IZ-1)*NXNY+IXY)+TINY))
              GFTSN =
     1          ABS(F(L0DEN+IXY)*F(L0C1+IXY)/
     1          (F(L0CSNT+JK1+(IZ-1)*NXNY+IXY)-TINY))
              GFTS = MIN(GFTSP,GFTSN)
C Multiplication by factor from Q1
              GFTS = GFTS*CHMRLX
              F(L0CO+IXY) = 1.0/(GFTS+TINY)
C  Fix values on first sweep because CHEM source is incorrect
C  until TEM1 section has been visited
              IF(ISWEEP.EQ.FSWEEP) F(L0CO+IXY)=GREAT
            ENDIF
1303      CONTINUE
        ENDIF
      ENDIF
      RETURN
  132 CONTINUE
C------------------- SECTION  3 ------------- coefficient = GRND2
      RETURN
  133 CONTINUE
C------------------- SECTION  4 ------------- coefficient = GRND3
      RETURN
  134 CONTINUE
C------------------- SECTION  5 ------------- coefficient = GRND4
      RETURN
  135 CONTINUE
C------------------- SECTION  6 ------------- coefficient = GRND5
      RETURN
  136 CONTINUE
C------------------- SECTION  7 ------------- coefficient = GRND6
      RETURN
  137 CONTINUE
C------------------- SECTION  8 ------------- coefficient = GRND7
C  Call GXPLAS for Plasma variables
      IF(JT0.GT.0) CALL GXPLAS
      IF(NPATCH(1:4).EQ.'SHWX') THEN
        L0CO = L0F(CO)
        CALL FN1(CO,0.0)
      ENDIF
      RETURN
  138 CONTINUE
C------------------- SECTION  9 ------------- coefficient = GRND8
      RETURN
  139 CONTINUE
C------------------- SECTION 10 ------------- coefficient = GRND9
      RETURN
 1310 CONTINUE
C------------------- SECTION 11 ------------- coefficient = GRND10
      RETURN
 1311 CONTINUE
C------------------- SECTION 12 ------------------- value = GRND
      RETURN
 1312 CONTINUE
C------------------- SECTION 13 ------------------- value = GRND1
C/////////////////////////////////////////////////////////////////
C Gas-phase chemistry
C/////////////////////////////////////////////////////////////////
      IF(NPATCH(1:4).EQ.'CHEM') THEN
        L0CO  = L0F(CO)
        L0VAL = L0F(VAL)
C
C Temperature
C
        IF(INDVAR.EQ.JTEM) THEN
          DO 1306 IXY = 1, NX*NY
            BLCKGE = .FALSE.
            IF(JVPOR.NE.0.AND.F(L0VPOR+IXY).LE.TINY) BLCKGE = .TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRP+IXY).GE.100) BLCKGE = .TRUE.
            IF(BLCKGE) THEN
              F(L0CO+IXY)  = TINY
              F(L0VAL+IXY) = TINY
            ELSE
              GTEMP = F(L0TEM+IXY) + TEMP0
              GPRES = F(L0P+IXY)+PRESS0
              IF(JT0.GT.0.AND.JION.GT.0) THEN
                L0T0 = L0F(JT0)
                L0ION = L0F(JION)
                PTEMP = F(L0T0+IXY)
                PION = F(L0ION+IXY)
              ENDIF
              DO 1307 JC = 1, JSP
                F(L0GC+JC) = F(JLC(JC) + IXY)
1307          CONTINUE
              CALL GCALF(JMAXS,JSP,F(L0GC+1),F(L0MW+1),F(L0MF+1))
C Calculate mole concentrations: mole/m^3
              DO 1323 JC=1,JSP
                FPRT(JC) = F(L0MF+JC)*GPRES/(GR*GTEMP+TINY)
1323          CONTINUE
              DO 1324 JR = 1, NGREAC
C Forward reaction rate
                RF(JR)=GAF(JR,GTEMP,GPRES,PTEMP,PION,JSP,FPRT,IGRT,
     1          GRPAR,90,GAPAR,GFTHB)
C Store raw rate for MNEWT
                IF(ISTIFF.EQ.1) F(L0RF+NGREAC*(IXY-1)+JR) = RF(JR)
                DO 1325 JC = 1, JSP
                  IF(IGAMMA(JC,JR) .LT. 0) THEN
                    RF(JR) =RF(JR)*(FPRT(JC)**(-IGAMMA(JC,JR)))
                  ENDIF
1325            CONTINUE
C Reverse reaction rate
                RR(JR)= GAR(JR,GTEMP,GPRES,PTEMP,PION,
     1                      IGAMMA,JSP,JCOEFF,GLOW,GHIGH,GTLOW,
     1                      GTAV,GTHIGH,JGAMMA,JMAXS,JHOMAX,FPRT,
     1                      IGRT,GRPAR,90,GAPAR,GFTHB)
C Store raw rate for MNEWT
                IF(ISTIFF.EQ.1) F(L0RR+NGREAC*(IXY-1)+JR) = RR(JR)
                DO 1326 JC = 1, JSP
                  IF(IGAMMA(JC,JR) .GT. 0) THEN
                    RR(JR) =RR(JR)*(FPRT(JC)**(+IGAMMA(JC,JR)))
                  ENDIF
1326            CONTINUE
1324          CONTINUE
C Linearisation
C Positive contributions are added to GSP, negative to GSN
              DO 1327 JC = 1, JSP
                GSP(JC) = 0.0
                GSN(JC) = 0.0
                DO 1328 JR = 1, NGREAC
                  IF(IGAMMA(JC,JR).GE.0) THEN
                    GSP(JC) = GSP(JC) + IGAMMA(JC,JR)*RF(JR)
                    GSN(JC) = GSN(JC) - IGAMMA(JC,JR)*RR(JR)
                  ELSE
                    GSP(JC) = GSP(JC) - IGAMMA(JC,JR)*RR(JR)
                    GSN(JC) = GSN(JC) + IGAMMA(JC,JR)*RF(JR)
                  ENDIF
1328            CONTINUE
                GSP(JC) = GSP(JC)*F(L0MW+JC)
                GSN(JC) = GSN(JC)*F(L0MW+JC)
                F(L0CSP3+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1            = GSP(JC)
                F(L0CSN3+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1            = GSN(JC)
                F(L0CSPT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1            = F(L0CSPT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)+GSP(JC)
                F(L0CSNT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1            = F(L0CSNT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)+GSN(JC)
                GCDOT(JC) = GSP(JC)+GSN(JC)
1327          CONTINUE
              GSUM = 0.0
              GCP = F(L0CP+IXY)
              DO 1329 JC = 1, JSP
C Formation enthalpy in molar units
                GHKMOL = H0(JCOEFF,GLOW(1,JC),GHIGH(1,JC),
     1                      GTLOW(JC),GTAV(JC),
     1                      GTHIGH(JC),250.0)
C Formation enthalpy in mass units
                GHKMAS = GHKMOL/(F(L0MW+JC)+TINY)
                GSUM = GSUM + GCDOT(JC)*GHKMAS
1329          CONTINUE
              GSOURCE = -GSUM
              F(L0VAL+IXY) = GSOURCE/TINY
            ENDIF
1306      CONTINUE
C Dumping of chemistry sources to optional stores
          IF(ISWEEP.EQ.LSWEEP) THEN
            DO 1308 JC = 1,JSP
              WRITE(JNAME(3:4),'(I2)') JC
              IF(JC.LT.10) JNAME(3:3) = '0'
              JNAME(1:2) = 'SP'
              JDIFF = LBNAME(JNAME)
              IF(STORE(JDIFF)) CALL FN0(JDIFF,
     1         -(L0CSP3+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY))
              JNAME(1:2) = 'SN'
              JDIFF = LBNAME(JNAME)
              IF(STORE(JDIFF)) CALL FN0(JDIFF,
     1         -(L0CSN3+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY))
1308        CONTINUE
          ENDIF
        ELSEIF(INDVAR.GE.C1.AND.INDVAR.LE.(C1+JMAXS-1).AND.INDVAR/=JTEM)
     1                                                              THEN
          L0C1 = L0F(INDVAR)
          DO 1330 IXY = 1, NX*NY
            BLCKGE = .FALSE.
            IF(JVPOR.NE.0.AND.F(L0VPOR+IXY).LE.TINY) BLCKGE = .TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRP+IXY).GE.100) BLCKGE = .TRUE.
            IF(BLCKGE) THEN
              F(L0CO+IXY) = 0.0
              F(L0VAL+IXY) = 0.0
              F(L0CSP3+JK1+(IZ-1)*NXNY+IXY) = 0.0
              F(L0CSN3+JK1+(IZ-1)*NXNY+IXY) = 0.0
              F(L0CSP3+(IZ-1)*NXNY+IXY) = 0.0
              F(L0CSN3+(IZ-1)*NXNY+IXY) = 0.0
            ELSE
              GC1 = F(L0C1+IXY)
              GSP1 = F(L0CSP3+JK1+(IZ-1)*NXNY+IXY)
              GSN1 = F(L0CSN3+JK1+(IZ-1)*NXNY+IXY)
C The linearised source term has the form
C    S = SP + (SN/c)*phi
C which is written as
C    S = Coef*(Val-phi).
C Thus
C    Val = -c*SP/SN ; Coef = - SN/c
C  No linearisation if stiff solver (Jacobian used elsewhere)
              IF(ISTIFF.EQ.0) THEN
                F(L0VAL+IXY) = -(GC1+TINY)*GSP1/(GSN1-TINY)
              ELSE
                F(L0VAL+IXY) = (GSP1+GSN1)/FIXFLU
              ENDIF
            ENDIF
1330      CONTINUE
        ENDIF
      ENDIF
C/////////////////////////////////////////////////////////////////
C Surface chemistry
C/////////////////////////////////////////////////////////////////
      IF(NPATCH(1:4).EQ.'SURF'.AND.JSP.GT.1) THEN
        FACTOR = 1.0
        DO 13375 ISR=1,ISURF
          IPN = F(L0SURF+ISR*2-1) + 0.1
          IF(IPN.EQ.IPNAME(NPATCH)) THEN
            FACTOR = F(L0SURF+ISR*2)
            GO TO 13376
          ENDIF
13375   CONTINUE
13376   CONTINUE
        IF(JNE.GT.0) L0NE = L0F(JNE)
        L0VAL = L0F(VAL)
        L0DEPO = L0F(JDEPO)
        L0GEOM = L0F(PATGEO)
        L0VOL = L0F(VOL)
        IF(BFC) THEN
          L0X = L0B(DHXPE)+NY*NX*(IZ-1)
          L0Y = L0B(DHYPN)+NY*NX*(IZ-1)
          L0Z = L0B(DHZPH)+NY*NX*(IZ-1)
        ELSE
          L0X=L0F(DXU2D)
          L0Y=L0F(DYV2D)
          CALL GETZ(ZGNZ,GZ,1000)
        ENDIF
        IF((INDVAR.GE.C1.AND.INDVAR.LE.(C1+JMAXS-1)).OR.
     1      (INDVAR.EQ.P1).OR.(INDVAR.EQ.JTEM)) THEN
          DO 1332 IY = IYF, IYL
          DO 1332 IX = IXF, IXL
            IXY = IY+NY*(IX-1)
            JCZ = (IZ-1)*NXNY + IXY
            BLCKGE = .FALSE.
            IF(JVPOR.NE.0.AND.F(L0VPOR+IXY).LE.TINY) BLCKGE = .TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRP+IXY).GE.100) BLCKGE = .TRUE.
            IF(BLCKGE) THEN
              F(L0VAL+IXY) = TINY
            ELSE
C Use solid temperature if any, otherwise gas
              IXYW = IXY
              L0TEMW = L0TEM
              GDIST = 1.0
C       east or west wall:
              IF(INTTYP.EQ.2.OR.INTTYP.EQ.3) THEN
                GDIST=0.5*F(L0X+IXY)
                IF(INTTYP.EQ.2) THEN
                  IF(XCYCLE.OR.IX.LT.NX) THEN
                    IXYW = IXY + NY
                    IF(IX.EQ.NX) IXYW = IY
                    IF(JVPOR.NE.0) THEN
                      IF(F(L0VPOR+IXYW).LT.TINY) IXYW = IXY
                    ENDIF
                    IF(JPRPS.NE.0) THEN
                      IF(QEQ(F(L0PRP+IXYW),PORPRP).OR.
     1                   QEQ(F(L0PRP+IXYW),VACPRP)) IXYW = IXY
                    ENDIF
                  ENDIF
                ELSE
                  IF(XCYCLE.OR.IX.GT.1) THEN
                    IXYW = IXY - NY
                    IF(IX.EQ.1) IXYW = (NX-1)*NY + IY
                    IF(JVPOR.NE.0) THEN
                      IF(F(L0VPOR+IXYW).LT.TINY) IXYW = IXY
                    ENDIF
                    IF(JPRPS.NE.0) THEN
                      IF(QEQ(F(L0PRP+IXYW),PORPRP).OR.
     1                   QEQ(F(L0PRP+IXYW),VACPRP)) IXYW = IXY
                    ENDIF
                  ENDIF
                ENDIF
C       north or south wall:
              ELSEIF(INTTYP.EQ.4.OR.INTTYP.EQ.5) THEN
                GDIST=0.5*F(L0Y+IXY)
                IF(INTTYP.EQ.4) THEN
                  IF(IY.LT.NY) THEN
                    IXYW = IXY + 1
                    IF(JVPOR.NE.0) THEN
                      IF(F(L0VPOR+IXYW).LT.TINY) IXYW = IXY
                    ENDIF
                    IF(JPRPS.NE.0) THEN
                      IF(QEQ(F(L0PRP+IXYW),PORPRP).OR.
     1                   QEQ(F(L0PRP+IXYW),VACPRP)) IXYW = IXY
                    ENDIF
                  ENDIF
                ELSE
                  IF(IY.GT.1) THEN
                    IXYW = IXY - 1
                    IF(JVPOR.NE.0) THEN
                      IF(F(L0VPOR+IXYW).LT.TINY) IXYW = IXY
                    ENDIF
                    IF(JPRPS.NE.0) THEN
                      IF(QEQ(F(L0PRP+IXYW),PORPRP).OR.
     1                   QEQ(F(L0PRP+IXYW),VACPRP)) IXYW = IXY
                    ENDIF
                  ENDIF
                ENDIF
C       high or low wall:
              ELSEIF(INTTYP.EQ.6.OR.INTTYP.EQ.7) THEN
                IF(BFC) THEN
                  GDIST=0.5*F(L0Z+IXY)
                ELSE
                  GDIST=0.5*DZ
                ENDIF
                IF(INTTYP.EQ.6) THEN
                  IF(IZ.LT.NZ) THEN
                    L0TEMW = L0F(ANYZ(JTEM,IZ+1))
                    IF(JVPOR.NE.0) THEN
                      L0VPOW = L0F(ANYZ(VPOR,IZ+1))
                      IF(F(L0VPOW+IXY).LT.TINY) L0TEMW = L0TEM
                    ENDIF
                    IF(JPRPS.NE.0) THEN
                      L0PRPW = L0F(ANYZ(PRPS,IZ+1))
                      IF(QEQ(F(L0PRPW+IXY),SOLPRP).OR.
     1                   QEQ(F(L0PRPW+IXY),VACPRP)) L0TEMW = L0TEM
                    ENDIF
                  ENDIF
                ELSE
                  IF(IZ.GT.1) THEN
                    L0TEMW = L0F(ANYZ(JTEM,IZ-1))
                    IF(JVPOR.NE.0) THEN
                      L0VPOW = L0F(ANYZ(VPOR,IZ-1))
                      IF(F(L0VPOW+IXY).LT.TINY) L0TEMW = L0TEM
                    ENDIF
                    IF(JPRPS.NE.0) THEN
                      L0PRPW = L0F(ANYZ(PRPS,IZ-1))
                      IF(QEQ(F(L0PRPW+IXY),SOLPRP).OR.
     1                   QEQ(F(L0PRPW+IXY),VACPRP)) L0TEMW = L0TEM
                    ENDIF
                  ENDIF
                ENDIF
              ENDIF
C-------------------------------------------------------------------------------
C           Test for the presence of adsorbed species.
C           No adsorbed species => use overall surface chemistry
              IF(JSU.EQ.0) THEN
C-------------------------------------------------------------------------------
              GTEMP = F(L0TEM+IXY) + TEMP0
              GTEMPW = F(L0TEMW+IXYW) + TEMP0
              GPRES = F(L0P+IXY)+PRESS0
              GRHO = F(L0DEN+IXY)
              IF(JNE.GT.0) PNE = F(L0NE+IXY)
              DO 1334 JC = 1,JSP
                GDIFW(JC) = F(L0DIF+(JC-1)*NX*NY*NZ
     1                          +(IZ-1)*NX*NY+IXY)
                F(L0GC+JC)= F(JLC(JC)+IXY)
C Diffusive flux of species JC
                GRDIF(JC) = 1./(F(L0MW+JC)+TINY)*GRHO*GDIFW(JC)*
     1                                            F(L0GC+JC)/GDIST
                GRDIF(JC) = MAX(TINY,GRDIF(JC))
1334          CONTINUE
C Function gflux calculates surface reaction rate
C GRKIN contains reaction rate in mol/m^2/s
              DO 1333 JR = 1,NSREAC
                GRKIN(JR) = GFLUX(JMAXS,JSP,JSUMAX,NSREAC,JR,
     1                      GTEMPW,GPRES,PNE,F(L0GC+1),F(L0MW+1),
     1                      ISRT,KGSP,SRPAR,90,SAPAR,GSPNAM,
     1                      GALFA)
1333          CONTINUE
              DO 1335 JR = 1,NSREAC
                IF(JDEP(JR).EQ.0.OR.JDEP(JR).EQ.4) THEN
                  GREFF(JR) = 0.0
                ELSEIF(JDEP(JR).EQ.1) THEN
                  GREFF(JR) = GRKIN(JR)
                ELSE
                  GREFF(JR) = 0.0
                  DO 1336 JC = 1,JSP
                    IF(GALFA(JC,JR).LT.-TINY) THEN
                      GREFF(JR) = GREFF(JR) - GALFA(JC,JR) /
     1                                      (GRDIF(JC) + TINY)
                    ENDIF
1336              CONTINUE
                  IF(JDEP(JR).EQ.2) THEN
                    GREFF(JR) = 1. / (GREFF(JR) + TINY)
                  ELSEIF(JDEP(JR).EQ.3) THEN
                    GREFF(JR) = GRKIN(JR) / (1. + GRKIN(JR)
     1                                      * GREFF(JR))
                  ELSE
                    WRITE(LUPR1,*) 'Invalid value of JDEP'
                    CALL SET_ERR(469, 'Error. Invalid value of JDEP.',1)
C                    STOP
                  ENDIF
                ENDIF
1335          CONTINUE
C
C Gas species
C
              IF(INDVAR.GE.C1.AND.INDVAR.LE.(C1+JMAXS-1).AND.
     1                                                INDVAR/=JTEM) THEN
                GSOURCE = 0.0
                DO 13371 JC=1,JSP
                DO 13371 JR=1,NSREAC
                  IF(GALFA(JC,JR).NE.0) GSOURCE = GSOURCE +
     1                           GALFA(JC,JR)*GREFF(JR)*F(L0MW+JC)
13371           CONTINUE
                GSOURCE = -GSOURCE*F(L0GC+JSUB)
                DO 1337 JR = 1, NSREAC
                  IF(GALFA(JSUB,JR).NE.0) GSOURCE = GSOURCE +
     1                           GALFA(JSUB,JR)*GREFF(JR)*F(L0MW+JSUB)
1337            CONTINUE
                F(L0VAL+IXY) = GSOURCE * FACTOR
C
C Mass
C
              ELSEIF(INDVAR.EQ.P1) THEN
                GSOURCE = 0.0
                DO 1338 JC = 1, JSP
                  DO 1339 JR = 1, NSREAC
                    IF(GALFA(JC,JR).NE.0) GSOURCE = GSOURCE+
     1                           GALFA(JC,JR)*GREFF(JR)*F(L0MW+JC)
1339              CONTINUE
1338            CONTINUE
                F(L0VAL+IXY) = GSOURCE * FACTOR
C Deposition rate
                GSURFL= 0.0
                DO 1341 JR = 1, NSREAC
                  IF(GSALFA(JR).NE.0) GSURFL = GSURFL +
     1                                      GSALFA(JR)*GREFF(JR)
1341            CONTINUE
                GSURFL= GSURFL*GMSOL*GFAC/GDSOL
                F(L0DEPO+IXY) = GSURFL
C Temperature
              ELSEIF(INDVAR.EQ.JTEM) THEN
                GSOURCE = 0.0
                DO 13373 JC=1,JSP
                  GSP(JC) = 0.0
                  GSN(JC) = 0.0
                  DSOURCE = 0.0
                  DO 13374 JR=1,NSREAC
                    IF(GALFA(JC,JR).NE.0) DSOURCE = DSOURCE +
     1                                        GALFA(JC,JR)*GREFF(JR)
                    IF(GALFA(JC,JR).GT.0) THEN
                      GSP(JC) = GSP(JC) +
     1                             GALFA(JC,JR)*GREFF(JR)*F(L0MW+JC)
                    ELSE
                      GSN(JC) = GSN(JC) +
     1                             GALFA(JC,JR)*GREFF(JR)*F(L0MW+JC)
                    ENDIF
13374             CONTINUE
                  VOLAREA = F(L0GEOM+IXY)/(F(L0VOL+IXY)+TINY)
                  F(L0CSPT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1             = F(L0CSPT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1                + GSP(JC)*VOLAREA
                  F(L0CSNT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1             = F(L0CSNT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1                + GSN(JC)*VOLAREA
                  FACT1 = H0(JCOEFF,GLOW(1,JC),GHIGH(1,JC),
     1                            GTLOW(JC),GTAV(JC),GTHIGH(JC),GTEMP)
                  FACT1 = FACT1 - H0(JCOEFF,GLOW(1,JC),GHIGH(1,JC),
     1                            GTLOW(JC),GTAV(JC),GTHIGH(JC),250.0)
                  GSOURCE = GSOURCE - DSOURCE*
     1                      (F(L0MW+JC)*F(L0CP+IXY)*GTEMP-FACT1)
13373           CONTINUE
                F(L0VAL+IXY) = GSOURCE * FACTOR
              ENDIF
C-----------------------------------------------------------------------
C           Test for the presence of adsorbed species.
C           adsorbed species present=> use detailed surface chemistry
              ELSE
C-----------------------------------------------------------------------
              GTEMP = F(L0TEM+IXY) + TEMP0
              GTEMPW = F(L0TEMW+IXYW) + TEMP0
              GPRES = F(L0P+IXY)+PRESS0
C Calculate fluxes
              IF(INDVAR.EQ.JSOL(JFST)) THEN
C Mass fractions
                DO 13400 JC = 1, JSP
                  F(L0GC+JC)= F(JLC(JC)+IXY)
                  GMAFR(JC) = F(JLC(JC)+IXY)
13400           CONTINUE
C Surface coverages
                DO 13406 JC = 1, JSU
                  GSFRA(JC) = F(JLC(JC+JSP)+IXY)
                  GSOCC(JC) = F(L0SOCC+JC)+0.1
13406           CONTINUE
                DO 13407 JC = 1, JTOT
                  GMOLW(JC) = F(L0MW+JC)
13407           CONTINUE
                ATOL = 1E-8
                RTOL = 1E-7
                ENDTIM = 0.5
C Call ODE solver driver
                CALL SURODE(JMAXS,JSUMAX,JCOEFF,
     1                      GTEMPW,GPRES,GMAFR,
     1                      GSFRA,GSFLUX,GSURFL,
     1                      NSREAC,JSP,JSU,GSOCC,GMOLW,GMSOL,
     1                      GTLOW,GTHIGH,GTAV,GLOW,
     1                      GHIGH,ISRT,
     1                      SRPAR,SAPAR,GALFA,GSALFA,KDEP,KVAC,
     1                      LUPR1,GSRCN,RTOL,ATOL,ENDTIM)
C Store surface fluxes
                DO 13402 JC = 1, JSP
                  F(L0FSU+(JC-1)*NX*NY*NZ+JCZ)=GSFLUX(JC)
13402           CONTINUE
C Store surface coverages
                DO 13403 JC = 1, JSU
                  F(JLC(JC+JSP)+IXY) = GSFRA(JC)
13403           CONTINUE
C Deposition rate
                GDEPR= GSURFL*GMSOL*GFAC/GDSOL
                F(L0DEPO+IXY) = GDEPR
C Write surface chemistry results  file
                IF(ISWEEP.EQ.LSWEEP) THEN
                  WRITE(3,*)
                  WRITE(3,*) NPATCH
                  WRITE(3,'(''IX ='',I3,'' IY ='',I3,'' IZ ='',I3)')
     1                  IX,IY,IZ
                  WRITE(3,'(''-----------------------------'')')
                  WRITE(3,'('' DEPOSITION FLUX = '',E12.3)') GSURFL
                  WRITE(3,'(''-----------------------------'')')
                  DO 13501 JC = 1, JSP
                    WRITE(3,13502) JC, GSFLUX(JC), F(JLC(JC)+IXY)
13501             CONTINUE
13502             FORMAT('FLUX(',I2,')=',E10.3,' MASSFR=',E10.3)
                  WRITE(3,'(''-----------------------------'')')
                  DO 13503 JC = 1, JSU
                    WRITE(3,13504) JC, GSFRA(JC)
13503             CONTINUE
13504             FORMAT('FRAC(',I2,')=',E10.3)
                  WRITE(3,'(''-----------------------------'')')
                ENDIF
              ENDIF
C Gas species fluxes
              IF(INDVAR.GE.C1.AND.INDVAR.LE.(C1+JMAXS-1).AND.
     1                                                INDVAR/=JTEM) THEN
                GSOURCE = 0.0
                DO 13401 JC=1,JSP
                  GSFLUX(JC) = F(L0FSU+(JC-1)*NX*NY*NZ+JCZ)
                  GSOURCE = GSOURCE +  GSFLUX(JC)*F(L0MW+JC)
13401           CONTINUE
                GSOURCE = -GSOURCE*F(L0GC+JSUB)
                GSOURCE = GSOURCE + GSFLUX(JSUB)*F(L0MW+JSUB)
                F(L0VAL+IXY) = GSOURCE*FACTOR
C Mass source
              ELSEIF(INDVAR.EQ.P1) THEN
                GSOURCE = 0.0
                DO 13405 JC = 1, JSP
                    GSOURCE = GSOURCE+F(L0FSU+(JC-1)*NX*NY*NZ+JCZ)
     1                     *F(L0MW+JC)
13405           CONTINUE
                F(L0VAL+IXY) = GSOURCE*FACTOR
C Temperature
              ELSEIF(INDVAR.EQ.JTEM) THEN
                GSOURCE = 0.0
                DO 13408 JC=1,JSP
                  GSP(JC) = 0.0
                  GSN(JC) = 0.0
                  DSOURCE = F(L0FSU+(JC-1)*NX*NY*NZ+JCZ)
                  IF(DSOURCE.GE.0.0) THEN
                      GSP(JC) = GSP(JC) + DSOURCE*F(L0MW+JC)
                  ELSE
                      GSN(JC) = GSN(JC) + DSOURCE*F(L0MW+JC)
                  ENDIF
                  VOLAREA = F(L0GEOM+IXY)/(F(L0VOL+IXY)+TINY)
                  F(L0CSPT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1             = F(L0CSPT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1                + GSP(JC)*VOLAREA
                  F(L0CSNT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1             = F(L0CSNT+(JC-1)*NX*NY*NZ+(IZ-1)*NX*NY+IXY)
     1                + GSN(JC)*VOLAREA
                  FACT1 = H0(JCOEFF,GLOW(1,JC),GHIGH(1,JC),
     1                            GTLOW(JC),GTAV(JC),GTHIGH(JC),GTEMP)
                  FACT1 = FACT1 - H0(JCOEFF,GLOW(1,JC),GHIGH(1,JC),
     1                            GTLOW(JC),GTAV(JC),GTHIGH(JC),250.0)
                  GSOURCE = GSOURCE - DSOURCE*
     1                      (F(L0MW+JC)*F(L0CP+IXY)*GTEMP-FACT1)
13408           CONTINUE
                F(L0VAL+IXY) = GSOURCE*FACTOR
              ENDIF
              ENDIF
C-------------------------------------------------------------------------------
            ENDIF
1332      CONTINUE
        ENDIF
      ENDIF
      RETURN
 1313 CONTINUE
C------------------- SECTION 14 ------------------- value = GRND2
      RETURN
 1314 CONTINUE
C------------------- SECTION 15 ------------------- value = GRND3
      RETURN
 1315 CONTINUE
C------------------- SECTION 16 ------------------- value = GRND4
      RETURN
 1316 CONTINUE
C------------------- SECTION 17 ------------------- value = GRND5
      RETURN
 1317 CONTINUE
C------------------- SECTION 18 ------------------- value = GRND6
      RETURN
 1318 CONTINUE
C------------------- SECTION 19 ------------------- value = GRND7
C  Call GXPLAS for Plasma variables
      IF(JT0.GT.0) CALL GXPLAS
      IF(NPATCH(1:4).EQ.'SHWX') THEN
        IF(.NOT.STORE(VISL)) THEN
          CALL WRYT40('STORE(ENUL) required for SHOWER         ')
          CALL WRIT40('STORE(ENUL) required for SHOWER         ')
          CALL SET_ERR(470, 'Error. STORE(ENUL) required for SHOWER.',1)
C          STOP
        ENDIF
        IF(.NOT.STORE(DEN1)) THEN
          CALL WRYT40('STORE(RHO1) required for SHOWER         ')
          CALL WRIT40('STORE(RHO1) required for SHOWER         ')
          CALL SET_ERR(471, 'Error. STORE(RHO1) required for SHOWER.',1)
C          STOP
        ENDIF
        DO 13701 ISH=1,ISHWR
          IPN = F(L0SHWR+ISH*3-2) + 0.1
          IF(IPN.EQ.IPNAME(NPATCH)) GO TO 13702
13701   CONTINUE
13702   CONTINUE
        CALL GETPTC(NPATCH,TYPE,IDUM,IDUM,IDUM,IDUM,IZF,IZL,IDUM,IDUM)
        ITYPE = TYPE + 0.1
        IF(ITYPE.EQ.1.OR.ITYPE.GT.7) THEN
          CALL WRYT40('SHOWER patch must have an AREA type     ')
          CALL WRIT40('SHOWER patch must have an AREA type     ')
          CALL SET_ERR(472,
     *          'Error. SHOWER patch must have an AREA type.',1)
C          STOP
        ENDIF
        L0VIS = L0F(VISL)
        L0VAR = L0F(INDVAR)
        L0VAL = L0F(VAL)
        CALL FN1(VAL,0.0)
        IF(ITYPE.EQ.2.OR.ITYPE.EQ.3) THEN
          DO 1372 IY=IYF,IYL
            ICF = (IXF-1)*NY + IY
            ICL = (IXL-1)*NY + IY
            FACTOR = 0.5*(F(L0VIS+ICF)*F(L0DEN+ICF)
     1                   +F(L0VIS+ICL)*F(L0DEN+ICL))
            FACTOR = F(L0SHWR+ISH*3)/F(L0SHWR+ISH*3-1)/FACTOR
            FLUX = FACTOR*(F(L0P+ICF)-F(L0P+ICL))
            DO 1373 IX=IXF,IXL
              IC = (IX-1)*NY + IY
              IF(IX.EQ.IXF) THEN
                IF(INDVAR.EQ.P1) THEN
                  IF(FLUX.LT.0.0) THEN
                    F(L0CO+IC) = FACTOR*F(L0DEN+ICL)
                  ELSE
                    F(L0CO+IC) = FACTOR*F(L0DEN+ICF)
                  ENDIF
                  F(L0VAL+IC) = F(L0P+ICL)
                ELSEIF(INDVAR.EQ.U1) THEN
                  IF(FLUX.LT.0.0) THEN
                    F(L0VAL+IC) = FLUX/F(L0SHWR+ISH*3)
                  ENDIF
                ELSEIF(INDVAR.GE.C1) THEN
                  F(L0VAL+IC) = F(L0VAR+ICL)
                ENDIF
              ELSEIF(IX.EQ.IXL) THEN
                IF(INDVAR.EQ.P1) THEN
                  IF(FLUX.LT.0.0) THEN
                    F(L0CO+IC) = FACTOR*F(L0DEN+ICL)
                  ELSE
                    F(L0CO+IC) = FACTOR*F(L0DEN+ICF)
                  ENDIF
                  F(L0VAL+IC) = F(L0P+ICF)
                ELSEIF(INDVAR.EQ.U1) THEN
                  IF(FLUX.GT.0.0) THEN
                    F(L0VAL+IC) = FLUX/F(L0SHWR+ISH*3)
                  ENDIF
                ELSEIF(INDVAR.GE.C1) THEN
                  F(L0VAL+IC) = F(L0VAR+ICF)
                ENDIF
              ENDIF
1373        CONTINUE
1372      CONTINUE
        ELSEIF(ITYPE.EQ.4.OR.ITYPE.EQ.5) THEN
          DO 1374 IX = IXF,IXL
            ICF = (IX-1)*NY + IYF
            ICL = (IX-1)*NY + IYL
            FACTOR = 0.5*(F(L0VIS+ICF)*F(L0DEN+ICF)
     1                   +F(L0VIS+ICL)*F(L0DEN+ICL))
            FACTOR = F(L0SHWR+ISH*3)/F(L0SHWR+ISH*3-1)/FACTOR
            FLUX = FACTOR*(F(L0P+ICF)-F(L0P+ICL))
            DO 1375 IY=IYF,IYL
              IC = (IX-1)*NY + IY
              IF(IY.EQ.IYF) THEN
                IF(INDVAR.EQ.P1) THEN
                  IF(FLUX.LT.0.0) THEN
                    F(L0CO+IC) = FACTOR*F(L0DEN+ICL)
                  ELSE
                    F(L0CO+IC) = FACTOR*F(L0DEN+ICF)
                  ENDIF
                  F(L0VAL+IC) = F(L0P+ICL)
                ELSEIF(INDVAR.EQ.V1) THEN
                  IF(FLUX.LT.0.0) THEN
                    F(L0VAL+IC) = FLUX/F(L0SHWR+ISH*3)
                  ENDIF
                ELSEIF(INDVAR.GE.C1) THEN
                  F(L0VAL+IC) = F(L0VAR+ICL)
                ENDIF
              ELSEIF(IY.EQ.IYL) THEN
                IF(INDVAR.EQ.P1) THEN
                  IF(FLUX.LT.0.0) THEN
                    F(L0CO+IC) = FACTOR*F(L0DEN+ICL)
                  ELSE
                    F(L0CO+IC) = FACTOR*F(L0DEN+ICF)
                  ENDIF
                  F(L0VAL+IC) = F(L0P+ICF)
                ELSEIF(INDVAR.EQ.V1) THEN
                  IF(FLUX.GT.0.0) THEN
                    F(L0VAL+IC) = FLUX/F(L0SHWR+ISH*3)
                  ENDIF
                ELSEIF(INDVAR.GE.C1) THEN
                  F(L0VAL+IC) = F(L0VAR+ICF)
                ENDIF
              ENDIF
1375        CONTINUE
1374      CONTINUE
        ELSEIF(ITYPE.EQ.6.OR.ITYPE.EQ.7) THEN
          IF(IZ.GT.IZF.AND.IZ.LT.IZL)  GO TO 1377
          L0VISF = L0F(ANYZ(VISL,IZF))
          L0VISL = L0F(ANYZ(VISL,IZL))
          L0DENF = L0F(ANYZ(DEN1,IZF))
          L0DENL = L0F(ANYZ(DEN1,IZL))
          L0PF = L0F(ANYZ(P1,IZF))
          L0PL = L0F(ANYZ(P1,IZL))
          L0VARF = L0F(ANYZ(INDVAR,IZF))
          L0VARL = L0F(ANYZ(INDVAR,IZL))
          DO 1376 IX=IXF,IXL
          DO 1376 IY=IYF,IYL
            IC = (IX-1)*NY + IY
            FACTOR = 0.5*(F(L0VISF+IC)*F(L0DENF+IC)
     1                   +F(L0VISL+IC)*F(L0DENL+IC))
            FACTOR = F(L0SHWR+ISH*3)/F(L0SHWR+ISH*3-1)/FACTOR
            FLUX = FACTOR*(F(L0PF+IC)-F(L0PL+IC))
            IF(IZ.EQ.IZF) THEN
              IF(INDVAR.EQ.P1) THEN
                IF(FLUX.LT.0.0) THEN
                  F(L0CO+IC) = FACTOR*F(L0DENL+IC)
                ELSE
                  F(L0CO+IC) = FACTOR*F(L0DENF+IC)
                ENDIF
                F(L0VAL+IC) = F(L0PL+IC)
              ELSEIF(INDVAR.EQ.W1) THEN
                IF(FLUX.LT.0.0) THEN
                  F(L0VAL+IC) = FLUX/F(L0SHWR+ISH*3)
                ENDIF
              ELSEIF(INDVAR.GE.C1) THEN
                F(L0VAL+IC) = F(L0VARL+IC)
              ENDIF
            ELSEIF(IZ.EQ.IZL) THEN
              IF(INDVAR.EQ.P1) THEN
                IF(FLUX.LT.0.0) THEN
                  F(L0CO+IC) = FACTOR*F(L0DENL+IC)
                ELSE
                  F(L0CO+IC) = FACTOR*F(L0DENF+IC)
                ENDIF
                F(L0VAL+IC) = F(L0PF+IC)
              ELSEIF(INDVAR.EQ.W1) THEN
                IF(FLUX.GT.0.0) THEN
                  F(L0VAL+IC) = FLUX/F(L0SHWR+ISH*3)
                ENDIF
              ELSEIF(INDVAR.GE.C1) THEN
                F(L0VAL+IC) = F(L0VARF+IC)
              ENDIF
            ENDIF
1376      CONTINUE
1377      CONTINUE
        ENDIF
      ENDIF
      RETURN
 1319 CONTINUE
C------------------- SECTION 20 ------------------- value = GRND8
      RETURN
 1320 CONTINUE
C------------------- SECTION 21 ------------------- value = GRND9
      RETURN
 1321 CONTINUE
C------------------- SECTION 22 ------------------- value = GRND10
      RETURN
C***************************************************************
C
C--- GROUP 14. Downstream pressure for PARAB=.TRUE.
C
   14 CONTINUE
      RETURN
C***************************************************************
C* Make changes to data for GROUPS 15, 16, 17, 18  GROUP 19.
C***************************************************************
C
C--- GROUP 19. Special calls to GROUND from EARTH
C
   19 GO TO (191,192,193,194,195,196,197,198,199,1910),ISC
  191 CONTINUE
C   * ------------------- SECTION 1 ---- Start of time step.
      RETURN
  192 CONTINUE
C   * ------------------- SECTION 2 ---- Start of sweep.
      IF(THMFRQ.GT.0) THEN
        JCHECK = MOD(ISWEEP,THMFRQ)
      ELSE
        JCHECK=0
      ENDIF
C Inititalise patch data for SHWR, SURF & VIN
      IF(ISWEEP.EQ.FSWEEP.AND.LINIL) THEN
        LINIL = .FALSE.
        DO 1921 ISH=1,ISHWR
          IPN = F(L0SHWR+ISH*3-2) + 0.1
          L0PTCH = I0PAT+(IPN-1)*10
          ITYPE = F(L0PTCH+2) + 0.1
          CALL GETCOV(NAMPAT(IPN),P1,GCOF,GVAL)
          IF(GCOF.LE.0.0) THEN
            CALL WRYT40('Positive P1 CO needed for shower patch  ')
            CALL WRIT40('Positive P1 CO needed for shower patch  ')
            CALL SET_ERR(473,
     *          'Error. Positive P1 CO needed for shower patch.',1)
C            STOP
          ENDIF
          F(L0SHWR+ISH*3-1) = GCOF
          F(I0PHI+2) = GRND7
          IF(ITYPE.EQ.2.OR.ITYPE.EQ.3) THEN
            CALL GETCOV(NAMPAT(IPN),U1,GCOF,GVAL)
          ELSEIF(ITYPE.EQ.4.OR.ITYPE.EQ.5) THEN
            CALL GETCOV(NAMPAT(IPN),V1,GCOF,GVAL)
          ELSEIF(ITYPE.EQ.6.OR.ITYPE.EQ.7) THEN
            CALL GETCOV(NAMPAT(IPN),W1,GCOF,GVAL)
          ENDIF
          IF(GCOF.LE.0.0) THEN
            CALL WRYT40('Positive vel CO needed for shower patch ')
            CALL WRIT40('Positive vel CO needed for shower patch ')
            CALL SET_ERR(474,
     *          'Error. Positive vel CO needed for shower patch.',1)
C            STOP
          ENDIF
          F(L0SHWR+ISH*3) = GCOF
          F(I0PHI+2) = GRND7
1921    CONTINUE
        DO 1922 ISR=1,ISURF
          IPN = F(L0SURF+ISR*2-1) + 0.1
          CALL GETCOV(NAMPAT(IPN),P1,GCOF,GVAL)
          IF(GCOF.LE.0.0) THEN
            CALL WRYT40('P1 CO must be positive for SURF patch   ')
            CALL WRYT40(' - value of 1.0 assumed                 ')
            CALL WRYT40('P1 CO must be positive for SURF patch   ')
            CALL WRIT40(' - value of 1.0 assumed                 ')
            GCOF = 1.0
          ENDIF
          F(L0SURF+ISR*2) = GCOF
          F(I0PHI+2) = FIXFLU
1922    CONTINUE
        DO 1923 IPN=1,NUMPAT
          IF(NAMPAT(IPN)(1:3).EQ.'VIN') THEN
C P1 VAL gives volumetric flow in sccm
C normal velocity VAL gives corresponding area
C species VALs are in volume fractions rather than mass fractions
            GMAV = F(L0MW+1)
            DO 1924 I=2,JSP
              CALL GETCOV(NAMPAT(IPN),JSOL(I),GCOF,GCDOT(I))
              GMAV = GMAV + GCDOT(I)*(F(L0MW+I)-F(L0MW+1))
1924        CONTINUE
            GCDOT(1) = 1.0
            DO 1925 I=2,JSP
              GCDOT(I) = GCDOT(I)*F(L0MW+I)/(GMAV+TINY)
              GCDOT(1) = GCDOT(1)-GCDOT(I)
1925        CONTINUE
            RHOIN = GMAV*1.01325E5/GR/273.15
            CALL GETCOV(NAMPAT(IPN),P1,GCOF,GVAL)
            GDUM1 = GVAL*RHOIN/60.E6
            CALL GETCOV(NAMPAT(IPN),JTEM,GCOF,GVAL)
            IF(PRESS0.LT.1.0) THEN
              CALL WRYT40('PRESS0 too low for VIN inlet to be used ')
              CALL WRIT40('PRESS0 too low for VIN inlet to be used ')
              CALL SET_ERR(475,
     *          'Error. PRESS0 too low for VIN inlet to be used.',1)
C              STOP
            ENDIF
            GDUM2 = GDUM1*GVAL*GR/GMAV/PRESS0
            L0PTCH = I0PAT+(IPN-1)*10
            ITYPE = F(L0PTCH+2) + 0.1
            IF(ITYPE.EQ.2.OR.ITYPE.EQ.3) THEN
              CALL GETCOV(NAMPAT(IPN),U1,GCOF,GVAL)
            ELSEIF(ITYPE.EQ.4.OR.ITYPE.EQ.5) THEN
              CALL GETCOV(NAMPAT(IPN),V1,GCOF,GVAL)
            ELSEIF(ITYPE.EQ.6.OR.ITYPE.EQ.7) THEN
              CALL GETCOV(NAMPAT(IPN),W1,GCOF,GVAL)
            ENDIF
            IF(GVAL.LE.0.0) THEN
              CALL WRYT40('Velocity VAL for VIN inlet should be +ve')
              CALL WRIT40('Velocity VAL for VIN inlet should be +ve')
              CALL SET_ERR(476,
     *          'Error. Velocity VAL for VIN inlet should be +ve.',1)
C              STOP
            ENDIF
            ITYPE = ITYPE - (ITYPE/2)*2
            FACTOR = (REAL(ITYPE)-0.5)/0.5
            GDUM1 = GDUM1/GVAL
            GDUM2 = GDUM2/GVAL
            F(I0PHI+2) = 0.0
            F(I0PHI+3) = GDUM2 * FACTOR
            CALL GETCOV(NAMPAT(IPN),P1,GCOF,GVAL)
            F(I0PHI+2) = FIXFLU
            F(I0PHI+3) = GDUM1
            DO 1926 I=2,JSP
              CALL GETCOV(NAMPAT(IPN),JSOL(I),GCOF,GVAL)
              F(I0PHI+2) = 0.0
              F(I0PHI+3) = GCDOT(I)
1926        CONTINUE
          ENDIF
1923    CONTINUE
      ENDIF
      RETURN
  193 CONTINUE
C   * ------------------- SECTION 3 ---- Start of iz slab.
C --------------------------------------------------------------------
C Set-up f-array zero location indicies for use in other groups.
C --------------------------------------------------------------------
      JZL               = IZ-1
      JZH               = IZ+1
      IF(IZ.EQ.1 ) JZL  = IZ
      IF(IZ.EQ.NZ) JZH  = IZ
      DO 1931 J1        = 1 , JSP+JSU
        JLC (J1)        = L0F(JSOL(J1)          )
        IF(NZ.GT.1) THEN
          JLCL(J1)      = L0F(ANYZ(JSOL(J1),JZL))
          JLCH(J1)      = L0F(ANYZ(JSOL(J1),JZH))
        ENDIF
 1931 CONTINUE
      L0TEM             = L0F(JTEM              )
      L0PRP             = L0F(JPRPS             )
      L0PRP             = 0
      L0PRPH            = 0
      L0PRPL            = 0
      IF(JPRPS.NE.0) THEN
        L0PRP           = L0F(JPRPS             )
        L0PRPH          = L0F(ANYZ(JPRPS,JZH)   )
        L0PRPL          = L0F(ANYZ(JPRPS,JZL)   )
      ENDIF
      L0VPOR            = 0
      L0VPOH            = 0
      L0VPOL            = 0
      IF(JVPOR.NE.0) THEN
        L0VPOR          = L0F(JVPOR             )
        L0VPOH          = L0F(ANYZ(JVPOR,JZH)   )
        L0VPOL          = L0F(ANYZ(JVPOR,JZL)   )
      ENDIF
      L0DEN             = L0F(DEN1              )
      L0P               = L0F(P1                )
      L0CP              = L0F(LCP1              )
C --------------------------------------------------------------------
C When STOREing one species, its mass concentration will be calculated
C from the difference between 1.0 and the sum of the other species.
C This is done at the (beginning) and end of each z-slab.
C --------------------------------------------------------------------
      RETURN
  194 CONTINUE
C   * ------------------- SECTION 4 ---- Start of iterations over slab.
C  Call GXPLAS for Plasma variables
      IF(JT0.GT.0) CALL GXPLAS
      RETURN
  199 CONTINUE
C   * ------------------- SECTION 9 ---- Start of solution sequence for
C                                                          a variable
C --------------------------------------------------------------------
C Subject species index.
C --------------------------------------------------------------------
      JSUB   = -999
      IF(INDVAR.GE.C1.AND.INDVAR.LE.C1+JMAXS-1.AND.INDVAR/=JTEM) THEN
C Species identifier from contiguous species counter.
        JSUB = JSOLV(INDVAR)
        IF(JSUB.GE.1.AND.JSUB.LE.JTOT) THEN
C Indices for accessing f-array for stores allocated by GXMAKE.
          JK1  = (JSUB-1)*NXNY*NZ
          JK2  = (JSUB-2)*NXNY*NZ
          JT1  = (JSUB-1)*NXNY
          JT2  = (JSUB-2)*NXNY
        ENDIF
      ELSEIF(INDVAR.EQ.JTEM) THEN
C Zero chemistry sources prior to calculation
        DO 19911 JC=1,JSP
        DO 19911 ICELL=1,NXNY
          J = (JC-1)*NXNY*NZ + (IZ-1)*NXNY + ICELL
          F(L0CSP3+J) = 0.0
          F(L0CSN3+J) = 0.0
          F(L0CSPT+J) = 0.0
          F(L0CSNT+J) = 0.0
19911   CONTINUE
      ENDIF
C --------------------------------------------------------------------
      RETURN
 1910 CONTINUE
C   * ------------------- SECTION 10---- Finish of solution sequence for
C                                                          a variable
      IF(ISTIFF.EQ.1) THEN
        IF(INDVAR.EQ.JSOL(JLST).AND.ISWEEP.GT.1) THEN
          L0SU = L0F(LSU)
          L0AP = L0F(LAP)
          IF(JVPOR.GT.0) L0VPOR=L0F(JVPOR)
          IF(JPRPS.GT.0) L0PRP=L0F(JPRPS)
          L0VOL = L0F(VOL)
C Loop over cells calling Newton-Raphson solver
          DO 91400 IXY = 1, NXNY
            IF(JVPOR.NE.0.AND.F(L0VPOR+IXY).GE.TINY) THEN
              IF(JPRPS.NE.0.AND.F(L0PRP+IXY).LT.100) THEN
                GTEMP = F(L0TEM+IXY) + TEMP0
                GPRES = F(L0P+IXY)+PRESS0
                GXVOL = F(L0VOL+IXY)
                GRHO = F(L0DEN+IXY)
                IF(JT0.GT.0.AND.JION.GT.0) THEN
                  L0T0 = L0F(JT0)
                  L0ION = L0F(JION)
                  PTEMP = F(L0T0+IXY)
                  PION = F(L0ION+IXY)
                ENDIF
C Read mass fractions
                DO 91408 JC = 1, JSP
                  F(L0GC+JC) = F(JLC(JC) + IXY)
91408           CONTINUE
C           Calculate mol fractions
                CALL GCALF(JMAXS,JSP,F(L0GC+1),F(L0MW+1),F(L0MF+1))
C Calculate molar concentrations fprt: mole/m^3
C ANTOL & RNTOL are termination criteria for solver
                DO 91407 JC=1,JSP
                  FPRT(JC) = F(L0MF+JC)*GPRES/(GR*GTEMP+TINY)
                  ANTOL(JC) = 1.0E-10*FPRT(JC)
                  RNTOL(JC) = 1.0E-3
91407           CONTINUE
C Read temporary coefficient storage for other species
                DO 91401 J = 1, JSP-1
                  GSUC(J) = F(L0SUK+(J-1)*NXNY+IXY)
                  GAPC(J) = F(L0APK+(J-1)*NXNY+IXY)
91401           CONTINUE
C Coefficients for last species
                GSUC(JSP)=F(L0SU+IXY)
                GAPC(JSP)=F(L0AP+IXY)
                DO 91406 JR = 1, NGREAC
                  RF(JR) = F(L0RF+NGREAC*(IXY-1)+JR)
                  RR(JR) = F(L0RR+NGREAC*(IXY-1)+JR)
91406           CONTINUE
C Number of Newton-Raphson iterations
                GTRIAL = 20
C Call Newton solver:
                CALL MNEWT(GTRIAL,F(L0GC+1),JSP,ANTOL,RNTOL,NDID,NGREAC,
     1                     GAPC,GSUC,RF,RR,IGAMMA,GXVOL,F(L0MW+1),
     1                     JMAXS,JHOMAX,GRHO,GABSOL,GRELAT,ANJAC)
C Calculate mass fraction of first species from the sum of the others
                GREL = 0.5
                SUM = 0.0
                DO 91402 J = 2, JSP
                  IF(F(L0GC+J).LT.VARMIN(JSOL(J)))
     1                   F(L0GC+J)=VARMIN(JSOL(J))
                  IF(F(L0GC+J).GT.VARMAX(JSOL(J)))
     1                   F(L0GC+J)=VARMAX(JSOL(J))
                  SUM = SUM + F(L0GC+J)
C Store mass fractions
                  F(JLC(J)+IXY) = GREL*F(L0GC+J)+
     1               (1.0-GREL)*F(JLC(J)+IXY)
91402           CONTINUE
                F(L0GC+1)= 1.0-SUM
                F(JLC(1)+IXY) = GREL*F(L0GC+1)+
     1               (1.0-GREL)*F(JLC(1)+IXY)
              ENDIF
            ENDIF
91400     CONTINUE
        ENDIF
        IF(INDVAR.GE.JSOL(1).AND.INDVAR.LT.JSOL(JLST)
     1         .AND.ISWEEP.GT.1) THEN
          L0SU = L0F(LSU)
          L0AP = L0F(LAP)
          JC = JSOLV(INDVAR)
          DO 91403 IXY = 1, NXNY
            BLCKGE = .FALSE.
            IF(JVPOR.NE.0.AND.F(L0VPOR+IXY).LT.TINY) BLCKGE = .TRUE.
            IF(JPRPS.NE.0.AND.F(L0PRP+IXY).GT.100) BLCKGE = .TRUE.
            IF(BLCKGE) THEN
C Put coefficients in temporary storage
                F(L0SUK+(JC-1)*NXNY+IXY) = 0.0
                F(L0APK+(JC-1)*NXNY+IXY) = 1.0E10
            ELSE
                F(L0SUK+(JC-1)*NXNY+IXY) = F(L0SU+IXY)
                F(L0APK+(JC-1)*NXNY+IXY) = F(L0AP+IXY)
            ENDIF
91403     CONTINUE
        ENDIF
      ENDIF
      RETURN
  195 CONTINUE
C   * ------------------- SECTION 5 ---- Finish of iterations over slab.
      RETURN
  196 CONTINUE
C   * ------------------- SECTION 6 ---- Finish of iz slab.
C  Call GXPLAS for Plasma variables
      IF(JT0.GT.0) CALL GXPLAS
C --------------------------------------------------------------------
C For non-BFC cases fill the LOW-area array for the next slab
C (declared by GXMAKE) with the HIGH-areas of this slab.
C --------------------------------------------------------------------
      IF(.NOT.BFC.AND.USOURC) THEN
        DO 2074 JC      = 1 , NXNY
          IF(IZ.EQ.NZ) THEN
            F(L0ARL+JC) = 0.0
          ELSE
            F(L0ARL+JC) = F(L0ARH+JC)
          ENDIF
 2074   CONTINUE
      ENDIF
C --------------------------------------------------------------------
C When STOREing one species, its mass concentration will be calculated
C from the difference between 1.0 and the sum of the other species.
C This is done at the beginning and (end) of each z-slab.
C --------------------------------------------------------------------
      IF(JST.NE.-999) THEN
        DO 1967 JC        = 1 , NXNY
          GSUM            = 0.0
          DO 1968 J1      = 2 , JSP
 1968     GSUM            = GSUM+F(JLC(J1)+JC)
          F(JLC(JST)+JC)  = MAX(0.0,1.0-GSUM)
C  Coding to zero mass fractions in solids (removing the need to
C  set CO and VAL in the CHEM patch to achieve the same effect)
          BLCKGE = .FALSE.
          IF(JVPOR.NE.0.AND.F(L0VPOR+JC).LE.TINY) BLCKGE = .TRUE.
          IF(JPRPS.NE.0.AND.F(L0PRP+JC).GE.100) BLCKGE = .TRUE.
          IF(BLCKGE) THEN
            DO 2081 J1    = 1 , JSP
2081        F(JLC(J1)+JC) = 0.0
          ENDIF
 1967   CONTINUE
      ENDIF
C
C     read mass fractions F-array indices
      IF(JVPOR.GT.0) L0VPOR = L0F(JVPOR)
      IF(JPRPS.GT.0) L0PRP = L0F(JPRPS)
C Mass fraction of unsolved species is calculated
      DO 1934 IXY = 1, NX*NY
        BLCKGE = .FALSE.
        IF(JVPOR.NE.0.AND.F(L0VPOR+IXY).LE.TINY) BLCKGE = .TRUE.
        IF(JPRPS.NE.0.AND.F(L0PRP+IXY).GE.100) BLCKGE = .TRUE.
        IF(BLCKGE) THEN
          F(JLC(1)+IXY) = 0.0
        ELSE
          GSUM = 0.0
          DO 1933 JC = 2, JSP
            GSUM = GSUM + F(JLC(JC)+IXY)
1933      CONTINUE
          F(JLC(1)+IXY) = MAX(0.0,1.0-GSUM)
        ENDIF
1934  CONTINUE
C Molar fraction output
      IF(ISWEEP.EQ.LSWEEP.OR.SWPRIN) THEN
        DO 2054 J1   = 1 , JSP
          JNAME(1:2) = 'FR'
C Check that molar fraction STORE exists.
C FR01 .... FR20.
          WRITE(JNAME(3:4),'(I2)') J1
          IF(J1.LT.10) JNAME(3:3) = '0'
          JDIFF                   = LBNAME(JNAME)
          IF(STORE(JDIFF)) THEN
            L0DUMP       = L0F(JDIFF)
            DO 2055 JC   = 1 , NXNY
              GSUM       = 0.0
              DO 2056 J2 = 1 , JSP
 2056         GSUM       = GSUM+ F(JLC(J2)+JC)/(F(L0MW+J2)+TINY)
              F(L0MF+J1) = (F(JLC(J1)+JC))/(F(L0MW+J1)+TINY)/(GSUM+TINY)
 2055       F(L0DUMP+JC) = F(L0MF+J1)
          ENDIF
 2054   CONTINUE
        IF(STORE(JSPHT)) THEN
          CALL FN0(JSPHT,LCP1)
        ENDIF
      ENDIF
C --------------------------------------------------------------------
      RETURN
  197 CONTINUE
C   * ------------------- SECTION 7 ---- Finish of sweep.
      RETURN
  198 CONTINUE
C   * ------------------- SECTION 8 ---- Finish of time step.
C
      RETURN
C***************************************************************
C
C--- GROUP 20. Preliminary print-out
C
   20 CONTINUE
      RETURN
C***************************************************************
C* Make changes to data for GROUPS 21 and 22 only in GROUP 19.
C***************************************************************
C
C--- GROUP 23. Field print-out and plot control
   23 CONTINUE
      RETURN
C***************************************************************
C
C--- GROUP 24. Dumps for restarts
C
   24 CONTINUE
      END
c