c

c          filename      GXCONVEC.HTM             070405
      SUBROUTINE GXCONV
C**********************************************************************
C    This subroutine applies corrections to the convective terms of the
C    momentum equation. The corrections are required if shocks or
C    density gradients due to other causes are present. The effect of
C    these corrections is to make the convective fluxes used in the
C    momentum equations consistent with those used in the continuity
C    equation. It is activated when UCONV=T and NAMGRD=CONV.
C
C    If the static enthalpy or static temperature equations are
C    solved, the correction applied to the momentum equations means
C    that total enthalpy is no longer preserved. This is corrected
C    by deactivating the built-in pressure term in the H1 equation,
C    and replacing it by a modifed version supplied in this routine.
C    The modified source is introduced via the following Q1 statements:
C
C        PATCH(ENTSOR,CELL,1,NX,1,NY,1,NZ,1,LSTEP)
C        COVAL(ENTSOR,H1 [or TEM1] ,FIXFLU,GRND1)
C
C    Transient cases require
C        PATCH(DPDT,VOLUME,1,NX,1,NY,1,NZ,1,LSTEP)
C        COVAL(DPDT, H1[or TEM1], FIXFLU, GRND1)
C
C    If very steep density gradients are present, convergence may be
C    improved by setting DENPCO=T, which fully accounts for density
C    gradients in the pressure correction equation.
C
C    The routine is called from GREX, and contains active code in
C    Group 1, Section 1, and Group 8, Section 8.
C
C    The present coding makes the following assumptions:
C    1. The density is STOREd as STORE(RHO1) in Q1.
C    2. The two phase option requires STORE(RHO2) in Q1.
C    3. BFC cases require STORE(UCMP,VCMP,WCMP) in Q1.
C    4. Two phase BFC cases require STORE(U2CM,V2CM,W2CM)
C
C    The modfications are described in detail in CHAM TR/147, and in
C    the Enyclopaedia entry 'Compressible Flow Corrections'.
C
C    Authors: Mike Malin, John Ludwig, CHAM UK
C    Date   : 24/12/91
C    Modification date: 1/07/94
C**********************************************************************
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/grdear'
      LOGICAL ERRO, HSOR, SLD
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB,ERRMES*40
      SAVE IU1,IU2,IV1,IV2,IW1,IW2,ITEM
      SAVE L0DUP,L0RUP,L0CON,L0CH,L0CL,L0AL,L0FE,L0FW,L0FN,L0FS,
     1 L0FH,L0FL,L0CLL
      SAVE HSOR
C
      NAMSUB = 'GXCONV'
C*****************************************************************
C
C---- GROUP 1. Run title and other preliminaries
C
      IF(IGR.EQ.1.AND.ISC.EQ.1) THEN
        IF(UCONV) THEN
C     User may here change message transmitted to the VDU screen or
C     batch-run log file.
          IF(.NOT.NULLPR) THEN
            CALL WRYT40('GXCONVEC OF 070405 HAS BEEN CALLED   ')
            CALL WRIT40('GXCONVEC OF 070405 HAS BEEN CALLED   ')
          ENDIF
C
C--- Reserve local storage
          CALL GXMAKE(L0DUP,NX*NY,'RHOUP')
          CALL GXMAKE(L0CON,NX*NY,'CON ')
          CALL GXMAKE(L0CH,NX*NY,'CH  ')
          CALL GXMAKE(L0CL,NX*NY,'CL  ')
          IF(NZ.GT.1) CALL GXMAKE(L0AL,NX*NY,'ALOW')
          HSOR=.FALSE.
          ITEM=LBNAME('TEM1')
C--- Stores for enthalpy pressure-term correction
          IF(SOLVE(H1).OR.SOLVE(ITEM)) THEN
            HSOR=.TRUE.
            IF(NX.GT.1) THEN
              CALL GXMAKE(L0FE,NX*NY*NZ,'FE  ')
              CALL GXMAKE(L0FW,NX*NY*NZ,'FW  ')
              CALL ZERNUM(L0FE,NX*NY*NZ)
              CALL ZERNUM(L0FW,NX*NY*NZ)
            ENDIF
            IF(NY.GT.1) THEN
              CALL GXMAKE(L0FN,NX*NY*NZ,'FN  ')
              CALL GXMAKE(L0FS,NX*NY*NZ,'FS  ')
              CALL ZERNUM(L0FN,NX*NY*NZ)
              CALL ZERNUM(L0FS,NX*NY*NZ)
            ENDIF
            IF(NZ.GT.1) THEN
              CALL GXMAKE(L0FH,NX*NY*NZ,'FH  ')
              CALL GXMAKE(L0FL,NX*NY*NZ,'FL  ')
              CALL ZERNUM(L0FH,NX*NY*NZ)
              CALL ZERNUM(L0FL,NX*NY*NZ)
              CALL GXMAKE(L0CLL,NX*NY,'CLL ')
              CALL FN1(-(L0FL+(NZ-1)*NX*NY),1.0)
            ENDIF
          ENDIF
C--- Find density stores
          ERRO=.NOT.STORE(DEN1)
          IF(ERRO) THEN
            ERRMES='STORE(RHO1);RHO1=GRNDn in Q1.           '
            GO TO 10
          ENDIF
          IF(.NOT.ONEPHS) THEN
            CALL GXMAKE(L0RUP,NX*NY,'RUP ')
            ERRO=.NOT.STORE(DEN2)
            IF(ERRO) THEN
              ERRMES='STORE(RHO2);RHO2=GRNDn in Q1.           '
              GO TO 10
            ENDIF
          ENDIF
C--- Select velocity stores - use components for BFC
          IF(BFC) THEN
            IU1=LBNAME('UCMP')
            IV1=LBNAME('VCMP')
            IW1=LBNAME('WCMP')
            ERRO=(NX.GT.1.AND.IU1.EQ.0).OR.(NY.GT.1.AND.IV1.EQ.0).OR.
     1         (NZ.GT.1.AND.IW1.EQ.0)
            IF(ERRO) THEN
              ERRMES='STORE(UCMP,VCMP,WCMP) in Q1.            '
              GO TO 10
            ENDIF
            IF(.NOT.ONEPHS) THEN
              IU2=LBNAME('U2CM')
              IV2=LBNAME('V2CM')
              IW2=LBNAME('W2CM')
              ERRO=(NX.GT.1.AND.IU2.EQ.0).OR.(NY.GT.1.AND.IV2.EQ.0).OR.
     1           (NZ.GT.1.AND.IW2.EQ.0)
              IF(ERRO) THEN
                ERRMES='STORE(U2CM,V2CM,W2CM) in Q1.            '
                GO TO 10
              ENDIF
            ENDIF
          ELSE
            IU1=U1
            IV1=V1
            IW1=W1
            IF(.NOT.ONEPHS) THEN
              IU2=U2
              IV2=V2
              IW2=W2
            ENDIF
          ENDIF
          RETURN
C--- Error exit
   10     CONTINUE
          CALL WRIT40('Shock correction requires:              ')
          CALL WRIT40(ERRMES)
          CALL WRIT40('Please correct Q1 and run again.        ')
          CALL SET_ERR(214,'Error. See result file.',2)
C          CALL EAROUT(2)
        ENDIF
C*****************************************************************
C
C--- GROUP 8. Terms (in differential equations) & devices
C
      ELSEIF(IGR.EQ.8.AND.ISC.EQ.8) THEN
C
C   * -----------GROUP 8  SECTION 8 --- CONVECTION COEFFICIENTS
C--- Entered when UCONV =.TRUE.
C
C--- Corrections applicable to velocities in compressible flow
        IF(INDVAR.EQ.W1.AND.(NDIREC.EQ.5.OR.NDIREC.EQ.6)) THEN
C--- Correct W1 fluxes
          CALL GXCNWH(NDIREC,IW1,DEN1,R1,L0DUP,L0RUP,L0CH,L0CL,L0AL,
     1                                      L0FH,L0FL,L0CLL,HSOR)
        ELSEIF(INDVAR.EQ.W2.AND.(NDIREC.EQ.5.OR.NDIREC.EQ.6)) THEN
C--- Correct W2 fluxes
          CALL GXCNWH(NDIREC,IW2,DEN2,R2,L0DUP,L0RUP,L0CH,L0CL,L0AL,
     1                                      L0FH,L0FL,L0CLL,HSOR)
        ELSEIF(INDVAR.EQ.V1.AND.NDIREC.EQ.1) THEN
C--- Correct V1 fluxes
          CALL GXCNVN(IV1,DEN1,R1,L0DUP,L0RUP,L0CH,L0CL,L0CON,
     1                                            L0FN,L0FS,HSOR)
        ELSEIF(INDVAR.EQ.V2.AND.NDIREC.EQ.1) THEN
C--- Correct V2 fluxes
          CALL GXCNVN(IV2,DEN2,R2,L0DUP,L0RUP,L0CH,L0CL,L0CON,
     1                                            L0FN,L0FL,HSOR)
        ELSEIF(INDVAR.EQ.U1.AND.NDIREC.EQ.3) THEN
C--- Correct U1 fluxes
          CALL GXCNUE(IU1,DEN1,R1,L0DUP,L0RUP,L0CH,L0CL,L0CON,
     1                                            L0FE,L0FW,HSOR)
        ELSEIF(INDVAR.EQ.U2.AND.NDIREC.EQ.3) THEN
C--- Correct U2 fluxes
          CALL GXCNUE(IU2,DEN2,R2,L0DUP,L0RUP,L0CH,L0CL,L0CON,
     1                                            L0FE,L0FW,HSOR)
        ENDIF
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
C                                       Index for Coefficient - CO
C                                       Index for Value       - VAL
      ELSEIF(IGR.EQ.13.AND.ISC.EQ.13) THEN
        IF(INDVAR.EQ.H1.OR.INDVAR.EQ.ITEM) THEN
C---Enthalpy/Temperature pressure source term consistent with momentum
C   modification.
          IF(NPATCH(1:4).EQ.'DPDT') THEN
C--- Add (P - Pold)/dt
            CALL FN1(VAL,0.0)
            IF(.NOT.STEADY) THEN
              L0P1=L0F(P1)
              L0PO=L0F(OLD(P1))
              L0VAL=L0F(VAL)
              DO I = 1,NX*NY
                IF(.NOT.SLD(I)) THEN
                  F(L0VAL+I)=HUNIT*(F(L0P1+I)-F(L0PO+I))/DT
                ENDIF
              ENDDO
            ENDIF
          ELSEIF(NPATCH(1:6).EQ.'ENTSOR') THEN
            CALL FN1(VAL,0.0)
C--- Add spatial contributions
            IF(NX.GT.1) CALL GXENTX(L0FE,L0FW,IU1,DEN1)
            IF(NY.GT.1) CALL GXENTY(L0FN,L0FS,IV1,DEN1)
            IF(NZ.GT.1) CALL GXENTZ(L0FH,L0FL,L0AL,IW1,DEN1)
          ENDIF
        ENDIF
      ENDIF
      END
C**********************************************************************
c
      SUBROUTINE GXCNWH(NDIR,IWVEL,IRHO,IR12,L0DUP,L0RUP,L0CH,L0CL,L0AL,
     1                                       L0FH,L0FL,L0CLL,HSOR)
C
C    This subroutine operates on the High/Low fluxes for W1
C    Input arguments:
C    NDIR     -  Direction indicator; 5 = High, 6 = Low
C    IWVEL    -  W velocity index
C    IRHO     -  Phase density index
C    IR12     -  Phase volume fraction index
C    L0DUP    -  Store for upwinded density
C    L0RUP    -  Store for upwinded volume fraction
C    L0CH     -  Store for convection at HIGH face of scalar cell
C    L0CL     -  Store for convection at LOW face of scalar cell
C    L0AL     -  Store for LOW cell face area
C    HSOR     -  Logical switch for enthalpy pressure factors
C
C    LD2      -  Convection flux index for High and Low
C**********************************************************************
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      COMMON/GENI/IGNI1(9),NFM,IGNI11(50)
      INTEGER HPORH
      LOGICAL HSOR
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB = 'GXCNWH'
C
C     LOW CONVECTIONS
C     ===============
C
      IF(NDIR.EQ.6) THEN
C--- Convection through low face of scalar cell
C... Upwind density into L0DUP
        CALL FN105(-L0DUP,LOW(IRHO),LOW(IWVEL),NFM)
        IF(.NOT.ONEPHS) THEN
C... For IPSA, upwind volume fraction into L0RUP, then multiply density
          CALL FN105(-L0RUP,LOW(IR12),LOW(IWVEL),NFM)
          CALL FN26(-L0DUP,-L0RUP)
        ENDIF
C... Now get convection Cl as rho*vel*area
        IF(BFC) THEN
          CALL FNGBFC(-L0AL,3,IZ-1)
        ELSE
          IF(STORE(HPOR)) THEN
            CALL FN21(-L0AL,LOW(HPOR),LXYAL,0.0,1.0)
          ELSE
            CALL FN0(-L0AL,AHIGH)
          ENDIF
        ENDIF
        CALL FN55(-L0CL,-L0DUP,LOW(IWVEL),-L0AL,1.0)
C--- Convection through high face of scalar cell
C... Upwind density into L0DUP
        CALL FN105(-L0DUP,IRHO,IWVEL,NFM)
        IF(.NOT.ONEPHS) THEN
C... For IPSA, upwind volume fraction into L0RUP, then multiply density
          CALL FN105(-L0RUP,IR12,IWVEL,NFM)
          CALL FN26(-L0DUP,-L0RUP)
        ENDIF
C... Now get convection Ch as rho*vel*area
        CALL FN55(-L0CH,-L0DUP,IWVEL,AHIGH,1.0)
C--- Compute convection CP for low face of w-cell by averaging scalars
        CALL FN10(LD2,-L0CL,-L0CH,0.,0.5,0.5)
C... Enthalpy pressure source factors
        IF(HSOR) THEN
           LLFL=L0FL+(IZ-1)*NX*NY
           IF(IZ.EQ.1) THEN
             CALL FN1(-LLFL,0.0)
           ELSEIF(IZ.EQ.2) THEN
             CALL FN1(-LLFL,1.0)
           ELSE
             CALL FN15(-LLFL,-L0CL,-L0CLL,0.0,1.0)
           ENDIF
           CALL FN0(-L0CLL,LD2)
        ENDIF
C... Finally remove outflows
        CALL FN22(LD2,0.0)
C
C     HIGH CONVECTIONS
C     ================
C
      ELSEIF(NDIR.EQ.5) THEN
C--- Convection through high face of scalar cell
        IF(IZ.EQ.1) THEN
C... On first Z slab, upwind density into L0DUP
          CALL FN105(-L0DUP,IRHO,IWVEL,NFM)
          IF(.NOT.ONEPHS) THEN
C... For IPSA, upwind volume fraction into L0RUP, then multiply density
            CALL FN105(-L0RUP,IR12,IWVEL,NFM)
          ENDIF
C... Now get convection Ch as rho*vel*area
          CALL FN55(-L0CL,-L0DUP,IWVEL,AHIGH,1.0)
        ELSE
C... For other slabs, copy previous slab's Ch
          CALL FN0(-L0CL,-L0CH)
        ENDIF
C--- Convection through high face of the high scalar cell
C... Upwind density into L0DUP
        CALL FN105(-L0DUP,HIGH(IRHO),HIGH(IWVEL),NFM)
        IF(.NOT.ONEPHS) THEN
C... For IPSA, upwind volume fraction into L0RUP, then multiply density
          CALL FN105(-L0RUP,HIGH(IR12),HIGH(IWVEL),NFM)
          CALL FN26(-L0DUP,-L0RUP)
        ENDIF
C... Now get convection Chh as rho*vel*area
        IF(.NOT.BFC.AND.STORE(HPOR)) THEN
C--- CARTES or POLAR, HPOR stored - must recover High area at IZ+2
          HPORH=ANYZ(HPOR,IZ+2)
          CALL FN55(-L0CH,HPORH,DXU2D,DYV2D,1.0)
          CALL FN55(-L0CH,-L0DUP,HIGH(IWVEL),-L0CH,1.0)
        ELSEIF(BFC) THEN
C--- BFC, so can get High area at IZ=2 directly
          CALL FNGBFC(-L0CH,APROJH,IZ+2)
          CALL FN55(-L0CH,-L0DUP,HIGH(IWVEL),-L0CH,1.0)
        ELSE
C--- CARTES or POLAR, HPOR not stored, so can use current IZ area
          CALL FN55(-L0CH,-L0DUP,HIGH(IWVEL),AHIGH,1.0)
        ENDIF
C--- Compute convection CH for high face of w-cell by averaging scalars
        CALL FN10(LD2,-L0CL,-L0CH,0.,0.5,0.5)
C... Enthalpy pressure source factors
        IF(HSOR) THEN
          LLFH=L0FH+(IZ-1)*NX*NY
          CALL FN15(-LLFH,-L0CL,LD2,0.0,1.0)
        ENDIF
C... Finally remove outflows
        CALL FN25(LD2,-1.0)
        CALL FN22(LD2,0.0)
      ENDIF
      END
C**********************************************************************
c
      SUBROUTINE GXCNVN(IVVEL,IRHO,IR12,L0DUP,L0RUP,L0CN,L0CS,L0CON,
     1                                               L0FN,L0FS,HSOR)
C
C    This subroutine operates on the North/South fluxes for V1
C    Input arguments:
C    IVVEL    -  V velocity index
C    IRHO     -  Phase density index
C    IR12     -  Phase volume fraction index
C    L0DUP    -  Store for upwinded density
C    L0RUP    -  Store for upwinded volume fraction
C    L0CN     -  Store for convection at NORTH face of velocity cell
C    L0CS     -  Store for convection at SOUTH face of velocity cell
C    L0CON    -  Store for convection at NORTH face of scalar cell
C    L0FS     -  Store for SOUTH face enthalpy pressure factor
C    L0FN     -  Store for NORTH face enthalpy pressure factor
C    HSOR     -  Logical switch for enthalpy pressure factors
C
C    LD11     -  Convection flux index for North
C    LD12     -  Convection flux index for South
C**********************************************************************
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/grdear'
C
      LOGICAL HSOR
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB = 'GXCNVN'
C
C     NORTH CONVECTIONS
C     =================
C
C--- Convection through north face of scalar cell
C... Upwind density into L0DUP
      CALL FN105(-L0DUP,IRHO,IVVEL,1)
      IF(.NOT.ONEPHS) THEN
C... For IPSA, upwind volume fraction into L0RUP, then multiply density
        CALL FN105(-L0RUP,IR12,IVVEL,1)
        CALL FN26(-L0DUP,-L0RUP)
      ENDIF
C... Now get convection as rho*vel*area
      CALL FN55(-L0CON,-L0DUP,IVVEL,ANORTH,1.0)
C... Compute convection for v through north face by averaging
      DO 92 IX=1,NX
        DO 90 IY=1,NY-2
          I=(IX-1)*NY+IY
          F(L0CN+I)=0.5*(F(L0CON+I)+F(L0CON+I+1))
   90   CONTINUE
        DO 91 IY=NY-1,NY
          I=(IX-1)*NY+IY
          F(L0CN+I)=0.0
   91   CONTINUE
   92 CONTINUE
C... Finally remove outflows
      CALL FN66(LD11,-L0CN,0.0,-1.0)
C
C     SOUTH CONVECTIONS
C     =================
C
C--- Convection through south face of scalar cell
      DO 101 IX=1,NX
        II=1+NY*(IX-1)
        F(L0CS+II)=0.0
        DO 101 IY=2,NY
          I=IY+NY*(IX-1)
          F(L0CS+I)=F(L0CN+I-1)
 101  CONTINUE
      CALL FN66(LD12,-L0CS,0.0,1.0)
C... Enthalpy pressure term factors
      IF(HSOR) THEN
        LLFN=L0FN+(IZ-1)*NX*NY
        LLFS=L0FS+(IZ-1)*NX*NY
        DO 102 IX=1,NX
          DO 102 IY=1,NY
            I=(IX-1)*NY+IY
            IF(IY.LE.2) THEN
              F(LLFS+I)=1.0
            ELSE
              F(LLFS+I)=F(L0CON+I-1)/(F(L0CN+I-2)+TINY)
            ENDIF
            IF(IY.GE.NY-1) THEN
              F(LLFN+I)=1.0
            ELSE
              F(LLFN+I)=F(L0CON+I)/(F(L0CN+I)+TINY)
            ENDIF
  102   CONTINUE
      ENDIF
      END
C**********************************************************************
c
      SUBROUTINE GXCNUE(IUVEL,IRHO,IR12,L0DUP,L0RUP,L0CE,L0CW,L0CON,
     1                                               L0FE,L0FW,HSOR)
C
C    This subroutine operates on the East/West fluxes for U1
C    Input arguments:
C    IUVEL    -  U velocity index
C    IRHO     -  Phase density index
C    IR12     -  Phase volume fraction index
C    L0DUP    -  Store for upwinded density
C    L0RUP    -  Store for upwinded volume fraction
C    L0CE     -  Store for convection at EAST face of velocity cell
C    L0CW     -  Store for convection at WEST face of velocity cell
C    L0CON    -  Store for convection at EAST face of scalar cell
C    L0FE     -  Store for EAST face enthalpy pressure factor
C    L0FW     -  Store for WEST face enthalpy pressure factor
C    HSOR     -  Logical switch for enthalpy pressure factors
C
C    LD11     -  Convection flux index for East
C    LD12     -  Convection flux index for West
C**********************************************************************
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/grdear'
C
      LOGICAL HSOR
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB = 'GXCNUE'
C
C     EAST CONVECTIONS
C     =================
C
C--- Convection through east face of scalar cell
C... Upwind density ito L0DUP
      CALL FN105(-L0DUP,IRHO,IUVEL,NY)
      IF(.NOT.ONEPHS) THEN
C... For IPSA, upwind volume fraction into L0RUP, then multiply density
        CALL FN105(-L0RUP,IR12,IUVEL,NY)
        CALL FN26(-L0DUP,-L0RUP)
      ENDIF
C... Now get convection as rho*vel*area
      CALL FN55(-L0CON,-L0DUP,IUVEL,AEAST,1.0)
C... Compute convection for U through east face by averaging
      LIM=ITWO(NX-1,NX-2,XCYCLE)
      DO 90 I=1,NY*LIM
        F(L0CE+I)=0.5*(F(L0CON+I)+F(L0CON+I+NY))
   90 CONTINUE
      IF(XCYCLE) THEN
        IADD=(NX-1)*NY
        DO 91 IY=1,NY
          F(L0CE+IY+IADD)=0.5*(F(L0CON+IY+IADD)+F(L0CON+IY))
   91   CONTINUE
      ELSE
        DO 92 IX=NX-1,NX
          DO 92 IY=1,NY
            I=(IX-1)*NY+IY
            F(L0CE+I)=0.0
   92   CONTINUE
      ENDIF
C... Finally remove outflows
      CALL FN66(LD11,-L0CE,0.0,-1.0)
C
C     WEST CONVECTIONS
C     =================
C
C--- Convection through west face of velocity cell
      IF(.NOT.XCYCLE) THEN
c... If not cyclic, set Cw(1) to 0.0
        DO 100 IY=1,NY
          F(L0CW+IY)=0.0
 100    CONTINUE
      ELSE
C... If cyclic, shift Ce(NX) to Cw(1)
        IADD=(NX-1)*NY
        DO 101 IY=1,NY
          F(L0CW+IY)=F(L0CE+IY+IADD)
 101    CONTINUE
      ENDIF
C... Now shift Ce(IX-1) to Cw(IX)
      DO 102 IX=2,NX
        DO 102 IY=1,NY
          I=IY+NY*(IX-1)
          F(L0CW+I)=F(L0CE+I-NY)
 102  CONTINUE
      CALL FN66(LD12,-L0CW,0.0,1.0)
C... Enthalpy pressure source term factors
      IF(HSOR) THEN
        LLFE=L0FE+(IZ-1)*NX*NY
        LLFW=L0FW+(IZ-1)*NX*NY
        DO 104 IX=1,NX
          DO 104 IY=1,NY
            I=(IX-1)*NY+IY
            IF(IX.LE.2) THEN
              F(LLFW+I)=1.0
            ELSE
              F(LLFW+I)=F(L0CON+I-NY)/(F(L0CE+I-NY*2)+TINY)
            ENDIF
            IF(IX.GE.NX-1) THEN
              F(LLFE+I)=1.0
            ELSE
              F(LLFE+I)=F(L0CON+I)/(F(L0CE+I)+TINY)
            ENDIF
  104   CONTINUE
      ENDIF
      END
C**********************************************************************
c
      SUBROUTINE GXENTX(L0FE,L0FW,IUVEL,IRHO)
C
C    This subroutine creates the East/West Enthalpy sources
C    Input arguments:
C    L0FE     -  Index for East face correction factor
C    L0FW     -  Index for West face correction factor
C    IUVEL    -  U velocity index
C    IRHO     -  Phase density index
C**********************************************************************
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      logical sld, ef, wf
C
      NAMSUB = 'GXENTX'
C
      L0VAL=L0F(VAL)
      L0U1=L0F(IUVEL)
      L0D1=L0F(IRHO)
      L0AE=L0F(AEAST)
      L0P1=L0F(P1)
      LLFE=L0FE+(IZ-1)*NX*NY
      LLFW=L0FW+(IZ-1)*NX*NY
      DO 131 I=1,NY*NX
        if(sld(i)) go to 131
        IX=(I-1)/NY+1
C... East face term
        IF(IX.LT.NX) THEN
          IF(IX.EQ.NX-1.AND..NOT.XCYCLE) THEN
            GE=F(L0D1+I+NY)*F(L0AE+I)*F(L0U1+I)
          ELSE
            GE=F(LLFE+I)*F(L0D1+I+NY)*F(L0AE+I)*0.5*(F(L0U1+I)+
     1                                            F(L0U1+I+NY))
          ENDIF
        ELSE
          IF(.NOT.XCYCLE) THEN
            GE=0.0
          ELSE
            IY=I-(IX-1)*NY
            GE=F(LLFE+I)*F(L0D1+IY)*F(L0AE+I)*0.5*(F(L0U1+I)+
     1                                            F(L0U1+IY))
          ENDIF
        ENDIF
C... West face term
        IF(IX.GT.1) THEN
          IF(IX.EQ.2) THEN
            IF(.NOT.XCYCLE) THEN
              GW=F(L0D1+I-NY)*F(L0AE+I-NY)*F(L0U1+I-NY)
            ELSE
              IADD=(NX-1)*NY
              IY=I-(IX-1)*NY
              GW=F(LLFW+I)*F(L0D1+I-NY)*F(L0AE+I-NY)*0.5*(F(L0U1+I-NY)
     1                                             +F(L0U1+IADD+IY))
            ENDIF
          ELSE
            GW=F(LLFW+I)*F(L0D1+I-NY)*F(L0AE+I-NY)*0.5*(F(L0U1+I-NY)
     1                                             +F(L0U1+I-2*NY))
          ENDIF
        ELSE
          IF(.NOT.XCYCLE) THEN
            GW=0.0
          ELSE
            IADD=(NX-1)*NY
            GW=F(LLFW+I)*F(L0D1+IADD)*F(L0AE+IADD)*0.5*(F(L0U1+IADD)
     1                                             +F(L0U1+IADD-NY))
          ENDIF
        ENDIF
C... Assemble E/W source
        IF(.NOT.WF(I)) THEN
          IF(IX.GT.1) THEN
            F(L0VAL+I)=F(L0VAL+I) +
     1     AMAX1( GW,0.0)*(F(L0P1+I)-F(L0P1+I-NY))*HUNIT/F(L0D1+I-NY)
          ELSEIF(XCYCLE) THEN
            IADD=(NX-1)*NY
            F(L0VAL+I)=F(L0VAL+I) +
     1     AMAX1( GW,0.0)*(F(L0P1+I)-F(L0P1+I+IADD))*HUNIT/
     2                                                 F(L0D1+I+IADD)
          ENDIF
        ENDIF
        IF(.NOT.EF(I)) THEN
          IF(IX.LT.NX) THEN
            F(L0VAL+I)=F(L0VAL+I) +
     1     AMAX1(-GE,0.0)*(F(L0P1+I)-F(L0P1+I+NY))*HUNIT/F(L0D1+I+NY)
          ELSEIF(XCYCLE) THEN
            IY=I-(IX-1)*NY
            F(L0VAL+I)=F(L0VAL+I) +
     1     AMAX1(-GE,0.0)*(F(L0P1+I)-F(L0P1+IY))*HUNIT/F(L0D1+IY)
          ENDIF
        ENDIF
  131 CONTINUE
      END
C**********************************************************************
c
      SUBROUTINE GXENTY(L0FN,L0FS,IVVEL,IRHO)
C
C    This subroutine creates the North/South Enthalpy sources
C    Input arguments:
C    L0FN     -  Index for North face correction factor
C    L0FS     -  Index for South face correction factor
C    IVVEL    -  V velocity index
C    IRHO     -  Phase density index
C**********************************************************************
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      logical sld, nf, sf
C
      NAMSUB = 'GXENTY'
C
      L0VAL=L0F(VAL)
      L0V1=L0F(IVVEL)
      L0D1=L0F(IRHO)
      L0AN=L0F(ANORTH)
      L0P1=L0F(P1)
      LLFN=L0FN+(IZ-1)*NX*NY
      LLFS=L0FS+(IZ-1)*NX*NY
      DO 131 IX=1,NX
        DO 131 IY=1,NY
          I=(IX-1)*NY+IY
          if(sld(i)) go to 131
C... North face term
          IF(IY.LT.NY) THEN
            IF(NY.LT.NY-1) THEN
              GN=F(LLFN+I)*F(L0D1+I+1)*F(L0AN+I)*0.5*(F(L0V1+I)+
     1                                             F(L0V1+I+1))
            ELSE
              GN=F(L0D1+I+1)*F(L0AN+I)*F(L0V1+I)
            ENDIF
          ELSE
            GN=0.0
          ENDIF
C... South face term
          IF(IY.GT.1) THEN
            IF(IY.EQ.2) THEN
              GS=F(L0D1+I-1)*F(L0AN+I-1)*F(L0V1+I-1)
            ELSE
              GS=F(LLFS+I)*F(L0D1+I-1)*F(L0AN+I-1)*0.5*(F(L0V1+I-1)
     1                                               +F(L0V1+I-2))
            ENDIF
          ELSE
            GS=0.0
          ENDIF
C... Assemble source
          IF(.NOT.SF(I)) THEN
            IF(IY.GT.1) THEN
              F(L0VAL+I)=F(L0VAL+I) +
     1        AMAX1( GS,0.0)*(F(L0P1+I)-F(L0P1+I-1))*HUNIT/F(L0D1+I-1)
            ENDIF
          ENDIF
          IF(.NOT.NF(I)) THEN
            IF(IY.LT.NY) THEN
              F(L0VAL+I)=F(L0VAL+I) +
     1        AMAX1(-GN,0.0)*(F(L0P1+I)-F(L0P1+I+1))*HUNIT/F(L0D1+I+1)
            ENDIF
          ENDIF
  131 CONTINUE
      END
C**********************************************************************
c
      SUBROUTINE GXENTZ(L0FH,L0FL,L0AL,IWVEL,IRHO)
C
C    This subroutine creates the High/Low Enthalpy sources
C    Input arguments:
C    L0FH     -  Index for High face correction factor
C    L0FL     -  Index for Low face correction factor
C    L0AL     -  Index for Low face area
C    IWVEL    -  U velocity index
C    IRHO     -  Phase density index
C**********************************************************************
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      logical sld, hf, lf
C
      NAMSUB = 'GXENTZ'
C
      L0VAL=L0F(VAL)
      L0W1=L0F(IWVEL)
      L0WH=L0F(HIGH(IWVEL))
      L0WL=L0F(LOW(IWVEL))
      L0WLL=L0F(ANYZ(IWVEL,IZ-2))
      L0D1=L0F(IRHO)
      L0DH=L0F(HIGH(IRHO))
      L0DL=L0F(LOW(IRHO))
      L0AH=L0F(AHIGH)
      L0P1=L0F(P1)
      L0PH=L0F(HIGH(P1))
      L0PL=L0F(LOW(P1))
      IF(BFC) THEN
        CALL FNGBFC(-L0AL,3,IZ-1)
      ELSE
        IF(STORE(HPOR)) THEN
          CALL FN21(-L0AL,LOW(HPOR),LXYAL,0.0,1.0)
        ELSE
          CALL FN0(-L0AL,AHIGH)
        ENDIF
      ENDIF
      LLFH=L0FH+(IZ-1)*NX*NY
      LLFL=L0FL+(IZ-1)*NX*NY
      call fn1(-llfl,1.0)
      call fn1(-llfh,1.0)
      DO 131 I=1,NX*NY
        if(sld(i)) go to 131
C... High Face
        IF(IZ.LT.NZ) THEN
          IF(IZ.EQ.NZ-1) THEN
            GH=F(L0DH+I)*F(L0AH+I)*F(L0W1+I)
          ELSE
            GH=F(LLFH+I)*F(L0DH+I)*F(L0AH+I)*0.5*(F(L0W1+I)+
     1                                             F(L0WH+I))
          ENDIF
        ELSE
          GH=0.0
        ENDIF
C... Low Face
        IF(IZ.GT.1) THEN
          IF(IZ.EQ.2) THEN
            GL=F(L0DL+I)*F(L0AL+I)*F(L0WL+I)
          ELSE
            GL=F(LLFL+I)*F(L0DL+I)*F(L0AL+I)*0.5*(F(L0WL+I)
     1                                           +F(L0WLL+I))
          ENDIF
        ELSE
          GL=0.0
        ENDIF
C... Assemble source
        IF(.NOT.HF(I)) THEN
          IF(IZ.LT.NZ) THEN
            F(L0VAL+I)=F(L0VAL+I) +
     1      AMAX1(-GH,0.0)*(F(L0P1+I)-F(L0PH+I))*HUNIT/F(L0DH+I)
          ENDIF
        ENDIF
        IF(.NOT.LF(I)) THEN
          IF(IZ.GT.1) THEN
            F(L0VAL+I)=F(L0VAL+I) +
     1      AMAX1( GL,0.0)*(F(L0P1+I)-F(L0PL+I))*HUNIT/F(L0DL+I)
          ENDIF
        ENDIF
  131 CONTINUE
      END
C**********************************************************************
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
      SUBROUTINE GXMACH(IMACH,IPTOT,IVABS,GAMA,IPH)
C.... This subroutine calculates the Mach Number, the total pressure
C     and the absolute velocity. These are stored in the indices
C     IMACH,IPTOT and IVABS respectively.
C
C     The ratio of specific heats, gamma, is passed in via the
C     argument GAMA.
C
C     For the built in options RHO1/RHO2 = GRND3/GRND5, it is set
C     automatically. Air is assumed if GAMA is zero.
C
C     The argument IPH sets the phase to 1 or 2.
C
C     It is called from Group 19.6 of GROUND.
C
C     Author : J C Ludwig, CHAM
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
      INCLUDE '/phoenics/d_includ/cmndmn'
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/fgemcb'
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      LOGICAL QEQ,EQZ,SLD
      SAVE IRUN0,L0VSQ,L0CSQ
      SAVE GAMMA
      DATA IRUN0 /0/
C
      NAMSUB = 'GXMACH'
C
      IF(STORE(IMACH)) THEN
        IF(IRUN0.NE.IRUN) THEN
          CALL GXMAKE(L0VSQ,NXYMAX,'VLSQ')
          CALL GXMAKE(L0CSQ,NXYMAX,'SNSQ')
          IRUN0=IRUN
        ENDIF
        GAMMA=GAMA
C... Compute VABS squared and store in L0VSQ
        CALL FNVLSQ(-L0VSQ,IPH)
C... Get VABS if store present
        IF(STORE(IVABS)) CALL FN71(IVABS,-L0VSQ,1.0,0.5)
C... Select density store
        IF(IPH.EQ.1) THEN
          IRHO=DEN1
          IF(QEQ(RHO1,GRND3)) GAMMA=1/RHO1B
          IF(QEQ(RHO1,GRND5)) GAMMA=1/RHO1C
        ELSE
          IRHO=DEN2
          IF(QEQ(RHO2,GRND3)) GAMMA=1/RHO2B
          IF(QEQ(RHO2,GRND5)) GAMMA=1/RHO2C
        ENDIF
C... If GAMMA not set, assume air
        IF(EQZ(GAMMA)) GAMMA=1.4
C... Calculate square of sonic velocity from Csq = GAMMA.P1/RHO
        CALL FN15(-L0CSQ,P1,IRHO,GAMMA*PRESS0,GAMMA)
C... Get Mach squared - Msq = VABSsq/Csq
        CALL FN15(IMACH,-L0VSQ,-L0CSQ,0.0,1.0)
C... Compute total pressure if P0 is STOREd
        IF(STORE(IPTOT)) THEN
          L0P0=L0F(IPTOT)
          L0P=L0F(P1)
          L0M=L0F(IMACH)
            ACON=GAMMA/(GAMMA-1.0)
            BCON=(GAMMA-1.0)/2.0
            DO 100 I=1,NX*NY
              IF(SLD(I)) GO TO 100
              AMACH=AMAX1(0.0,AMIN1(100.0,F(L0M+I)))
              F(L0P0+I)=ABSPRES(I) *
     1                  ABS((1.0 + BCON*AMACH)) ** ACON
  100       CONTINUE
        ENDIF
C... Now get square root of Mach Number
        CALL FN30(IMACH)
      ENDIF
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C     SUBROUTINE GXCMPR is called from section 8, group 8 of GREX3, and
C     is entered, for velocities only, when UCONV = T and COMPRS = T.
C
C.... The argument L0CNV is the F-array index for the convection
C     coefficient; L0P and L0NEXP are P1 and P1(next), where 'next'
C     is north, east or high; GA is a constant, set in GREX3 .
C
      SUBROUTINE GXCMPR(LCNV,LP,LNEXTP,GA)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      COMMON /IDATA/NX,NY,NZ,LUPR1,IGFILL(116)
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB = 'GXCMPR'
      G1 = (GA-1.0)/GA
      G2 = (1.0+1.0/GA)*0.3333333
      G3 = 2.0*GA
      L0CNV = L0F(LCNV)
      L0P = L0F(LP)
      L0NEXP = L0F(LNEXTP)
      IXLAST=IXL
      IF(INDVAR.EQ.3.OR.INDVAR.EQ.4) IXLAST=MIN0(IXL,NX-1)
      IYLAST=IYL
      IF(INDVAR.EQ.5.OR.INDVAR.EQ.6) IYLAST=MIN0(IYL,NY-1)
      DO 10 IX = IXF,IXLAST
        IADD = NY * (IX-1)
        DO 10 IY = IYF,IYLAST
          I = IY + IADD
          INCONV = L0CNV + I
          IF(F(INCONV).GT.0.0) THEN
            PRAT = (F(L0NEXP+I) + PRESS0)/ABSPRES(I)
            PRATM1 = 1.0 - PRAT
            IF(ABS(PRATM1).LT.0.03) THEN
              TOP = PRATM1*(1.0+G2*PRATM1)
              F(INCONV) = F(INCONV)*(1.0-TOP/(G3+TOP))
            ELSE
              PRATG1=PRAT**G1
              F(INCONV) = F(INCONV)*G1*PRATM1/(1.0-PRATG1)
            ENDIF
          ENDIF
   10 CONTINUE
      NAMSUB = 'gxcmpr'
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c