cGxhocs.htm c c

c
C     File name ............. GXHOCS.FOR ....................... 110817
C     File includes the following subroutines and functions:
C          GXHOCS, GXHOEN, GXHOHL, BLKSLD, GXFLPS
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
      SUBROUTINE GXHOCS
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C     SUBROUTINE GXHOCS introduces
c 
c     higher-order convection schemes for
C     use with the staggered-mesh option in PHOENICS. These are
C     discretization schemes for the convection terms which have
C     truncation errors of order greater than one (i.e. higher than
C     First-Order Upwind). These schemes allow more accurate resolution
C     of convection-dominated flows, especially when these contain
C     gradients of the flow variables which are not aligned with the
C     grid directions (e.g. recirculating flows).
C
C     Two classes of schemes have been implemented:
C
C     LINEAR SCHEMES: - These include the Central-Differencing scheme
C     (CDS), the Linear-Upwind scheme (LUS), the Quadratic-Upwind
C     scheme (QUICK), the Cubic-Upwind scheme (CUS) and Fromm's scheme.
C     These offer generally good resolution, but do not guarantee
C     boundedness and may therefore be prone to spurious oscillations
C     such as "wiggles" around shock waves or negative values of KE or
C     EP.
C
C     NON-LINEAR FLUX-LIMITED SCHEMES: - These schemes also offer
C     improved resolution but are designed to limit the higher-order
C     correction so as to guarantee positive coefficients and therefore
C     bounded behaviour. This limiting may reduce the accuracy of the
C     discretization to some extent. These schemes modify themselves
C     according to the local flow conditions and are therefore also
C     termed "non-linear". The flux limiters that have been implemented
C     are:
C       SMART       (bounded QUICK)
C       Koren       (bounded CUS, TVD)
C       MUSCL       (van Leer: bounded Fromm, TVD)
C       HQUICK      (based on QUICK,   smooth/harmonic)
C       OSPRE       (smooth/continuous)
C       van Leer harmonic (smooth/harmonic, TVD)
C       van Albada  (smooth/continuous)
C       Minmod      (= Zhu/Rodi SOUCUP scheme, TVD)
C       Superbee    (compressive limiter for sharp gradients, TVD)
C       UMIST       (bounded QUICK, TVD )
C       HCUS        (based on CUS, smooth/harmonic)
C       CHARM       (based on QUICK and UDS, polymonial ratio)
C
C     These two classes, as implemented, offer a choice of 17
C     different schemes, though not all of these are in fact recom-
C     mended for general use. Those recommended are:
C
C     (i)  Linear: either CUS or QUICK. LUS is less accurate, but much
C          more stable.
C
C     (ii) Non-linear: SMART and Koren give good accuracy with strong
C          relaxation. MUSCL and H-QUICK are less accurate but more
C          stable. OSPRE is the most stable but again a little less
C          accurate. van Leer harmonic, van Albada, Minmod and Superbee
C          are classical limiters that are included for completeness.
C
C   GXHOCS is called from Group 1, Sections 1 and 2, and also from
C   Group 13 of GREX3:
C
C     ISC=13 Section 13 :    COVAL(HOCS,PHI,FIXFLU,GRND1)
C                            Linear - KAPPA schemes
C
C     ISC=14 Section 14 :    COVAL(HOCS,PHI,FIXFLU,GRND2)
C                            Non-Linear - Flux-Limited Positive schemes
C
C   The integer array INLCS(1) in GXHOCS contains the scheme number
C   for each variable PHI. The LINEAR schemes available are:
C       INLCS(phi) =  1 -> LUS     =  2 ->  FROMM   =  3  -> CUS
C                  =  4 -> QUICK   =  5 ->  CDS
C     The NON-LINEAR schemes available are:
C       INLCS(phi) =  6 -> SMART   =  7 ->  KOREN   =  8 ->  MUSCL
C                  =  9 -> H-QUICK =  10 -> OSPRE
C                  =  11 ->  VAN LEER HARMONIC/HLPA
C                  =  12 ->  VAN ALBADA =  13 ->  MINMOD/SOUCUP
C                  =  14 ->  SUPERBEE   =  15 ->  UMIST
C                  =  16 ->  H-CUS      =  17 ->  CHARM
C
C   Q1-file activation:
C
C   A particular scheme is assigned to one or more dependent variables
C   by the SCHEME command, whose syntax is:
C
C     SCHEME(NAME,variable name 1,variable name 2,...etc.)
C
C   The first argument NAME identifies the required discretisation
C   scheme, as follows:
C   NAME = LUS, FROMM, CUS, QUICK, CDS, SMART, KOREN, MUSCL, HQUICK,
C          OSPRE, VANLH, VANALB, MINMOD, SUPBEE, UMIST, HCUS or CHARM.
C
C   The second argument permits the user to specify those SOLVEd
C   variables which will use the selected scheme. If ALL is entered as
C   the second argument, then the selected scheme is applied to all
C   SOLVEd-for variables. For example,
C
C     SCHEME(QUICK,U1,V1);SCHEME(SMART,H1,C1,C2)
C
C   selects QUICK for U1 and V1, and SMART for H1,C1 and C2, and
C   Upstream differencing (as DIFCUT=0) for any SOLVEd variables which
C   do not appear in a SCHEME command. The INLCS array is determined
C   by the SCHEME command, and then transmitted to GXHOCS via the
C   SPEDAT facility. The schemes may also be activated via the MENU.
C
C   The schemes are described under the Encyclopaedia entry 'SCHEMES
C   FOR CONVECTION DISCRETISATION' ( see also the HELP entry for the PIL
C   command SCHEME ).
C
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
      INCLUDE 'patcmn'
      INCLUDE 'farray'
      INCLUDE 'satear'
      COMMON /IDA10/  INLCS(150)
      COMMON /LB/     LB1(16),LZDZ,LB2(8),LZDZG,LB3(3),DXU2D,DYV2D,
     1                LB4(8),DXG2D,DYG2D,LB5(96),LRHO1,LRHO1H,LB6,
     1                LRHO2,LRHO2H,LB7(3),AEAST,ANORTH,AHIGH,LB8(46),
     1                LCEA,LCNA,LCHA,LCEA2,LCNA2,LCHA2,LB9(58),CON1E,
     1                CON2E,CON1N,CON2N,CON1H,CON2H,LB10(36)
     1       /GENI/   NXNY,IGN1(8),NFM,IGN2(38),NPHI4,IGN3(11)
     1       /IGE/    IXF,IXL,IYF,IYL,IGE1(2),IGR,ISC,IRUN,IZ,
     1                IGE2,ISWEEP,ISTEP,INDVAR,VAL,IGE4(10)
     1       /CGE/    NPATCH
     1       /BFCIN1/ L0B(106),LULAST
     1       /BFCINT/ APROJE,APROJN,APROJH,IBFC1(7),DHXPE,IBFC2,DHYPN,
     1                IBFC3,DHZPH,IBFC4(12),DZGPH,IBFC5(39),XP,YP,ZP,
     1                IBFC6(36)
     1       /GHOCS/  L0PHI,L0VAL,L0RUA,L0RHO,L0RHOH,L0R,L0FPR,L0PRP,
     1                L0CVL,L0UCT,L0VCT,L0WCT,L0XP,L0YP,L0ZP,L0DXG,
     1                L0DYG,L0DZG,KAPP(150),KAPM(150),U1ORU2,V1ORV2,
     1                W1ORW2,UVW,STACON,ISCHEM,DEBHOC,debhcs
     1       /IZLIM/  IZF,IZL
     1       /RSG/    RGFIL1(19),RSG20,RGFIL2(180)
     1       /ISG/    ISGF1(57),ISCHM1,ISGF2(42)
      COMMON/LMAIN2/  DBFIL1,DBGSWP,DBGTZ,DBFIL2(8)
      COMMON/BFICRT/  ICRT(6),IUCMP(2),IVCMP(2),IWCMP(2)
      COMMON /NAMFN/  NAMFUN,NAMSUB
      COMMON/MDFCNV/LIMCNV ! convection was modified by in-form
      LOGICAL LIMCNV
C
      CHARACTER NPATCH*8,NAMFUN*6,NAMSUB*6,CTEMP*8
      logical debhcs,debhoc,dbgswp,dbgtz,dbfil1,dbfil2
      LOGICAL INIRUA,U1ORU2,V1ORV2,W1ORW2,UVW,FLUID1,STACON,SOLVE
      INTEGER VAL,AEAST,ANORTH,AHIGH,APROJN,
     1        APROJH,DXU2D,DYV2D,DXG2D,DYG2D,DHXPE,DHYPN,DHZPH,DZGPH,
     1        CON1E,CON1N,CON1H,CON2E,CON2N,CON2H,APROJE,XP,YP,ZP
      REAL KAPP,KAPM,KAPPA
      CHARACTER*6 SCNAM(17)
      DATA SCNAM /'LUS', 'FROMM', 'CUS', 'QUICK',
     1            'CDS', 'SMART', 'KOREN', 'VANL1',
     1            'HQUICK', 'OSPRE', 'VANL2', 'VANALB',
     1            'MINMOD', 'SUPBEE', 'UMIST', 'HCUS',
     1            'CHARM'/
      SAVE INIRUA
 
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      NAMSUB='GXHOCS'
C
      IF(IGR==1 .AND. ISC==1) THEN
C
C>>>>>>>>>>>>>>>>>>>>> Group  1 ------- Section  1 <<<<<<<<<<<<<<<<<<<<<
C
C.... Construct INLCS array from SPEDAT facility & then set parameter
C     values for KAPPA schemes.
        if(.not.nullpr) WRITE(LUPR1,101)
 101    FORMAT(2X,'NO',3X,'VARIABLE',4X,'SCHEME')
        ONED3=1./3.
        DO MPHI=1,NPHI
          INLCS(MPHI)=0
          IF(MPHI>2 .AND. SOLVE(MPHI)) THEN
            WRITE(CTEMP,'(I3.3)') MPHI
            CALL GETSDI('SCHEME','INLCS'//CTEMP,INLCS(MPHI))
            IF(INLCS(MPHI)/=0) THEN
              IF(.NOT.NULLPR) WRITE(LUPR1,102)
     1              MPHI,NAME(MPHI),SCNAM(INLCS(MPHI))
            ELSE
              CYCLE
            ENDIF
            IF(INLCS(MPHI)<6.OR.INLCS(MPHI)==10.OR.
     1        INLCS(MPHI)==12) THEN
              IF(.NOT.NULLPR) THEN
                CALL WRIT40('Scheme '//SCNAM(INLCS(MPHI)))
                CALL WRIT40('Use of this scheme is not recommended')
                CALL WRIT40('See PHOENICS Encyclopaedia entry on  ')
                CALL WRIT40('"Schemes of discretization"          ')
              ENDIF
            ENDIF
            IF(INLCS(MPHI)==1) THEN
              KAPPA=-1.
            ELSEIF(INLCS(MPHI)==2) THEN
              KAPPA=0.
            ELSEIF(INLCS(MPHI)==3) THEN
              KAPPA=ONED3
            ELSEIF(INLCS(MPHI)==4) THEN
              KAPPA=0.5
            ELSEIF(INLCS(MPHI)==5) THEN
              KAPPA=1.
            ENDIF
            IF(INLCS(MPHI)>=1.AND.INLCS(MPHI)<=5) THEN
              KAPP(MPHI) = 0.5*(1.0 + KAPPA)
              KAPM(MPHI) = 0.5*(1.0 - KAPPA)
            ELSE
              KAPP(MPHI) = 0.
              KAPM(MPHI) = 0.
            ENDIF
          ENDIF
        ENDDO
        IF(ISCHM1>1) THEN
          WRITE(LUPR1,'('' Schemes become active after sweep '',I6)')
     1                                                           ISCHM1
        ELSE
          WRITE(LUPR1,'('' Schemes are active for all sweeps '')')
        ENDIF
        CALL WRITST
 102    FORMAT(1X,I3,5X,A4,5X,A6)
C
C.... STACON=T signifies that in-line convections for velocity
C     will be computed by averaging scalar convection fluxes.
        STACON=UCONV.AND.NAMGRD=='CONV'
        CALL GETSDL('SCHEME','DEBHOC',DEBHOC)
      ELSEIF(IGR==13) THEN
C
C>>>>>>>>>>>>>>>>>>>>> Group  13 ------- Section 13-14 <<<<<<<<<<<<<<<<
C
        IF(NPATCH(1:4)/='HOCS') RETURN
C
        ISCHEM=INLCS(INDVAR)
        STACON = STACON.OR.LIMCNV
        debhcs=debhoc.and.dbgswp.and.dbgtz.and.dbgphi(indvar)
C.... Initialize values of source term to zero on slab
        CALL FN1(VAL,0.0)
        IF(STEADY.AND.FIINIT(1)/=READFI) THEN ! always active for transient & restart
          IF(ISWEEP1) L0RHOH = L0F(LRHO1H)
          IF(.NOT.ONEPHS) L0R = L0F(9)
        ELSE
          L0RHO = L0F(DEN2)
          IF(NZ>1) L0RHOH = L0F(LRHO2H)
          IF(.NOT.ONEPHS) L0R = L0F(10)
        ENDIF
C.... Set cell-geometry addresses
        IF(NX>1) THEN
          L0DXU=L0F(DXU2D); L0DXG=L0F(DXG2D); L0AE=L0F(AEAST)
          IF(BFC) L0DXU = L0B(DHXPE)+IZNXNY
        ENDIF
        IF(NY>1) THEN
          L0DYV=L0F(DYV2D); L0DYG=L0F(DYG2D); L0AN=L0F(ANORTH)
          IF(BFC) L0DYV = L0B(DHYPN)+IZNXNY
        ENDIF
        IF(NZ>1) THEN
          IF(.NOT.BFC) THEN
            L0AH = L0F(AHIGH)
            L0DZW=L0F(LZDZ); L0DZG=L0F(LZDZG)
          ELSE
            L0AH = L0B(APROJH)+IZNXNY
            L0DZG = L0B(DZGPH)+IZNXNY
            L0DZW = L0B(DHZPH)+IZNXNY
          ENDIF
        ENDIF
C
C....   N.B. Distance of cell centre from cell face in covariant
C       components is not available: DHXPE etc. are NORMAL to cell faces.
C       (DHXPE,DHYPN,DHZPH are accessed here but not used.)
C       Hence in BFC case it is assumed faces are midway between centres.
        IF(BFC) THEN
          L0DXU=L0DXG; L0DYV=L0DYG; L0DZW=L0DZG
        ENDIF
C
C.... Set offset indices for 1st/2nd-phase velocities and mass fluxes
        IF(FLUID1(INDVAR)) THEN
          IADD = 0; IAD2 = 0
        ELSE
          IADD = 1; IAD2 = 6
        ENDIF
C.... Store original PATCH limits
C     (Limits IZF,IZL are not defined and hence GETPTC required)
        CALL GETPTC(NPATCH,PATYP,IXF,IXL,IYF,IYL,IZF,IZL,ITF,ITL)
        IXFP=IXF; IXLP=IXL; IYFP=IYF; IYLP=IYL; IZFP=IZF; IZLP=IZL
C
C.... Check for momentum variables and reset patch limits appropriately
        U1ORU2=.FALSE.
        V1ORV2=.FALSE.
        W1ORW2=.FALSE.
        IF(INDVAR==3.OR.INDVAR==4) THEN
          U1ORU2=.TRUE.
          IXL = MIN(IXL,(NX-1))
        ELSEIF(INDVAR==5.OR.INDVAR==6) THEN
          V1ORV2=.TRUE.
          IYL = MIN(IYL,(NY-1))
        ELSEIF(INDVAR==7.OR.INDVAR==8) THEN
          W1ORW2=.TRUE.
          IZL = MIN(IZL,(NZ-1))
        ENDIF
        UVW = U1ORU2.OR.V1ORV2.OR.W1ORW2
        IF(BFC.AND.UVW) THEN
          IF(FLUID1(INDVAR)) THEN
            IF(NX>1) JUCMP = IUCMP(1)
            IF(NY>1) JVCMP = IVCMP(1)
            IF(NZ>1) JWCMP = IWCMP(1)
            IUCRT=ICRT(1); IVCRT=ICRT(3); IWCRT=ICRT(5)
          ELSE
            IF(NX>1) JUCMP = IUCMP(2)
            IF(NY>1) JVCMP = IVCMP(2)
            IF(NZ>1) JWCMP = IWCMP(2)
            IUCRT=ICRT(2); IVCRT=ICRT(4); IWCRT=ICRT(6)
          ENDIF
          L0UCT=L0F(IUCRT); L0VCT=L0F(IVCRT); L0WCT=L0F(IWCRT)
          L0XP=L0B(XP); L0YP=L0B(YP); L0ZP=L0B(ZP)
          L0XP = L0XP + IZNXNY
          L0YP = L0YP + IZNXNY
          L0ZP = L0ZP + IZNXNY
        ENDIF
C
        INIRUA=(STEADY.AND.ISWEEP==FSWEEP).OR.(.NOT.STEADY.AND.
     1          ISTEP==FSTEP.AND.ISWEEP==FSWEEP)
C... NORTH cell-face contribution to HOCS source term
        IF(NY>1) THEN
          IF(.NOT.(BFC.AND.UVW)) THEN
            L0CVL = L0F(5+IADD)
          ELSE
            L0CVL = L0F(JVCMP)
          ENDIF
          IF(FLUID1(INDVAR)) THEN
            LCONN = ITWO(LCNA,CON1N,NZ==1)
            L0RUA = L0F(LCONN) + IZNXNY
          ELSE
            LCONN = ITWO(LCNA2,CON2N,NZ==1)
            L0RUA = L0F(LCONN) + IZNXNY
          ENDIF
          IF(INIRUA.AND.NZ>1.AND.W1ORW2) THEN
             DO J=1,NXNY
               F(L0RUA+NXNY+J)=0.0
             ENDDO
          ENDIF
          IYF=IYF+1; IYL=IYL-2
          CALL GXHOEN(ISC,L0DYG,L0DYV,L0AN,L0DXG,V1ORV2,U1ORU2,1,NY,1)
          IYF=IYFP; IYL=IYLP
        ENDIF
C
C... HIGH/LOW cell-face contribution to HOCS source term
        IF(NZ>1) THEN
          IF(BFC.AND.W1ORW2) THEN
            L0CVL = L0F(JWCMP)
          ELSE
            L0CVL = L0F(7+IADD)
          ENDIF
C.... HIGH face
C     ---------
          IF(FLUID1(INDVAR)) THEN
            L0RUA = L0F(CON1H) + IZNXNY
          ELSE
            L0RUA = L0F(CON2H) + IZNXNY
          ENDIF
          IF(INIRUA.AND.NZ>1.AND.W1ORW2) THEN
            DO J=1,NXNY
              F(L0RUA+NXNY+J)=0.0
            ENDDO
          ENDIF
          IF(IZ>=IZF+1 .AND. IZ<=IZL-2) THEN
            CALL GXHOHL(5,ISC,L0DZW,L0DZG,L0AH)
          ENDIF
C.... LOW face
C     --------
          IF(FLUID1(INDVAR)) THEN
            L0RUA = L0F(CON1H) + (IZ-2)*NXNY
          ELSE
            L0RUA = L0F(CON2H) + (IZ-2)*NXNY
          ENDIF
c          L0CVL = L0CVL-NFM
c          L0PHI = L0PHI-NFM
c          IF(BFC)
c     1      CALL SUB3(L0UCT,L0UCT-NFM,L0VCT,L0VCT-NFM,L0WCT,L0WCT-NFM)
          IF(IZ>=IZF+2 .AND. IZ<=IZL-1) THEN
            L0CVL = L0CVL-NFM
            L0PHI = L0PHI-NFM
            IF(BFC) THEN
              L0UCT=L0UCT-NFM; L0VCT=L0VCT-NFM; L0WCT=L0WCT-NFM
            ENDIF
            CALL GXHOHL(6,ISC,L0DZW,L0DZG,L0AH)
            L0CVL = L0CVL+NFM
            L0PHI = L0PHI+NFM
            IF(BFC) THEN
              L0UCT=L0UCT+NFM; L0VCT=L0VCT+NFM; L0WCT=L0WCT+NFM
            ENDIF
          ENDIF
          IZF=IZFP; IZL=IZLP
        ENDIF
C
C... EAST cell-face contribution to HOCS source term
        IF(NX>1) THEN
          IF(.NOT.(BFC.AND.UVW)) THEN
            L0CVL = L0F(3+IADD)
          ELSE
            L0CVL = L0F(JUCMP)
          ENDIF
          IF(FLUID1(INDVAR)) THEN
            LCONE = ITWO(LCEA,CON1E,NZ==1)
            L0RUA = L0F(LCONE) + IZNXNY
          ELSE
            LCONE = ITWO(LCEA2,CON2E,NZ==1)
            L0RUA = L0F(LCONE) + IZNXNY
          ENDIF
          IF(INIRUA.AND.NZ>1.AND.W1ORW2) THEN
             DO J=1,NXNY
               F(L0RUA+NXNY+J)=0.0
             ENDDO
          ENDIF
          IXF=IXF+1; IXL=IXL-2
          CALL GXHOEN(ISC,L0DXG,L0DXU,L0AE,L0DYG,U1ORU2,V1ORV2,NY,1,3)
          IXF=IXFP; IXL=IXLP
        ENDIF
      ENDIF
      NAMSUB='gxhocs'
      END
C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
      SUBROUTINE GXHOEN(ISC,L0DG,L0DV,L0AF,L0DGN,UVI,UVN,INDI,INDN,IDIR)
C *********************************************************************
C     Purpose: to calculate higher-order correction for EAST or NORTH
C              faces on a single slab of cells
C                       east face     north face
C     Input:  L0DG  =     L0DXG         L0DYG
C             L0DV  =     L0DXU         L0DYV
C             L0AF  =     L0AE          L0AN
C             L0DGN =     L0DYG         L0DXG
C             UVI   =     U1ORU2        V1ORV2
C             UVN   =     V1ORV2        U1ORU2
C             INDI  =     NY            1
C             INDN  =     1             NY
C             IDIR  =     3             1
C     Called by: GXHOCS
C *********************************************************************
      INCLUDE 'farray'
      COMMON
     1       /IDATA/  NX,NY,NZ,LUPR1,IDAT1(116)
     1       /LDATA/  CARTES,LD1(2),ONEPHS,LD2(3),XCYCLE,LD3(10),STEADY,
     1                BFC,LD4(44),NONORT,LD5(19)
     1       /IGE/    IXF,IXL,IYF,IYL,IGE1(5),IZ,
     1                IGE3,ISWEEP,IGE4,INDVAR,IGE5(11)
     1       /GENI/   NXNY,IGN1(8),NFM,IGN2(38),NPHI4,IGN3(11)
     1       /GHOCS/  L0PHI,L0VAL,L0RUA,L0RHO,L0RHOH,L0R,L0FPR,L0PRP,
     1                L0CVL,L0UCT,L0VCT,L0WCT,L0XP,L0YP,L0ZP,LBDXG,
     1                LBDYG,LBDZG,KAPP(150),KAPM(150),U1ORU2,V1ORV2,
     1                W1ORW2,UVW,STACON,ISCHEM,DEBHOC,debhcs
     1       /NAMFN/  NAMFUN,NAMSUB
     1       /RDATA/  TINY,RDAT(84)
 
C
      logical debhoc,debhcs
      LOGICAL BLKSLD,DOCOR,U1ORU2,V1ORV2,W1ORW2,CARTES,UVI,UVN,UVW,LD1,
     1        XCYCLE,LD2,STEADY,BFC,LD3,LD4,LD5,NONORT,QGE,ONEPHS,
     1        STACON
      REAL KAPP,KAPM
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB = 'GXHOEN'
C
C.... Cycle over the patch.
      debhcs=.false.
      if(debhcs) call banner(1,namsub,110817)
      IZNXNY=(IZ-1)*NXNY
      DO 10 IX=IXF,IXL
        IADD = (IX-1)*NY
        DO 10 IY=IYF,IYL
          if(debhcs) call writ2i('ix      ',ix,'iy      ',iy)
          IXY = IY+IADD
C.... Check for blocked faces or solid cells within stencil
          DOCOR=.NOT.BLKSLD(UVI,UVN,W1ORW2,INDI,INDN,NXNY,IXY+IZNXNY,
     1        IDIR)
          if(debhcs) call writ1l('docor   ',docor)
C.... If no blocked/solid cells/faces within stencil then continue
          IF(DOCOR) THEN
C
C.... Set mass-flux coefficient for face
            RUA = F(L0RUA+IXY)
            if(debhcs) call writ3l('uvi   ',uvi,'uvn   ',uvn,
     1                             'stacon  ',stacon)
            IF(UVI) THEN
              IF(STACON) THEN
                RUA = 0.5*(RUA + F(L0RUA+IXY+INDI))
              ELSE
                IF(F(L0CVL+IXY)>0.0) THEN
                  GRO = F(L0RHO+IXY)
                  IF(.NOT.ONEPHS) GVF = F(L0R+IXY)
                ELSE
                  GRO = F(L0RHO+IXY+INDI)
                  IF(.NOT.ONEPHS) GVF = F(L0R+IXY+INDI)
                ENDIF
                GUF = 0.5*(F(L0CVL+IXY) + F(L0CVL+IXY+INDI))
                RUA = GRO*GUF*F(L0AF+IXY)
                IF(.NOT.ONEPHS) RUA=RUA*GVF
              ENDIF
            ELSEIF(UVN) THEN
              RUA = 0.5*(RUA + F(L0RUA+IXY+INDN))
            ELSEIF(W1ORW2) THEN
C N109 fails because f(l0rua+ixy+nxny) has 1e19 on 1st sweep
c see mrm 06.05.05 initialise con1n & con1e for 1st sweep
              RUA = 0.5*(RUA + F(L0RUA+IXY+NXNY))
            ENDIF
C
C.... Values of convected variable on each side of face
            FP  = F(L0PHI+IXY)
            FEN = F(L0PHI+IXY+INDI)
C
C.... Calculate required cell-geometry data
            IF(UVI) THEN
              DXC = F(L0DV+IXY+INDI)
              if(dxc<=0.0) then
                call writ3i('indvar  ',indvar,'ix      ',ix,
     1                      'iy      ',iy)
                CALL SET_ERR(444, 'Error. See reslt file.',1)
C                stop
              endif
              IF(QGE(RUA,0.0)) THEN
                DXU = F(L0DV+IXY)
                GDX = 0.5*F(L0DV+IXY+INDI)
              ELSE
                DXU = F(L0DV+IXY+2*INDI)
                GDX = -0.5*F(L0DV+IXY+INDI)
              ENDIF
            ELSE
              DXC = F(L0DG+IXY)
              IF(QGE(RUA,0.0)) THEN
                DXU = F(L0DG+IXY-INDI)
                GDX = 0.5*F(L0DV+IXY)
              ELSE
                DXU = F(L0DG+IXY+INDI)
                GDX = -0.5*F(L0DV+IXY+INDI)
              ENDIF
            ENDIF
C
C.... Calculate gradients of convected variable
            IF(.NOT.(BFC.AND.NONORT.AND.UVW)) THEN
C
              GRPHC = (FEN - FP)/(DXC+tiny)
              IF(QGE(RUA,0.0)) THEN
                fpnei = F(L0PHI+IXY-INDI)
                GRPHU = (FP - F(L0PHI+IXY-INDI))/(DXU+tiny)
              ELSE
                fpnei = F(L0PHI+IXY+2*INDI)
                GRPHU = (F(L0PHI+IXY+2*INDI) - FEN)/(DXU+tiny)
              ENDIF
C
            ELSE
C
              IF(UVI) THEN
                INDB = INDI
                DP = F(L0DG+IXY)
                DE = F(L0DG+IXY+INDI)
              if(debhcs) call writ2r('dp     1',dp,'de      ',de)
              ELSEIF(UVN) THEN
                INDB = INDN
                DP = F(L0DGN+IXY)
                DE = F(L0DGN+IXY+INDI)
                if(debhcs) then
                  call writ2r('dp     2',dp,'de      ',de)
                  call writ3i('l0dgn   ',l0dgn,'ixy    ',ixy,
     1                        'indi    ',indi)
                endif
              ELSE
                INDB = NXNY
                DP = F(LBDZG+IXY)
                DE = F(LBDZG+IXY+INDI)
              if(debhcs) call writ2r('dp     3',dp,'de      ',de)
              ENDIF
C
              if(de<=0.0.or.dp<=0) then
                call writ3i('indvar  ',indvar,'ix      ',ix,
     1                      'iy      ',iy)
                call writ2r('dp      ',dp,'de      ',de)
                CALL SET_ERR(511, 'Error. See result file.',1)
C                call wayout(1)
              endif
              DP=DP+TINY
              DE=DE+TINY
              DPX = (F(L0XP+IXY+INDB)-F(L0XP+IXY))/DP
              DPY = (F(L0YP+IXY+INDB)-F(L0YP+IXY))/DP
              DPZ = (F(L0ZP+IXY+INDB)-F(L0ZP+IXY))/DP
C
              DEX = (F(L0XP+IXY+INDI+INDB)-F(L0XP+IXY+INDI))/DE
              DEY = (F(L0YP+IXY+INDI+INDB)-F(L0YP+IXY+INDI))/DE
              DEZ = (F(L0ZP+IXY+INDI+INDB)-F(L0ZP+IXY+INDI))/DE
C
              FENP = F(L0UCT+IXY+INDI)*DPX+F(L0VCT+IXY+INDI)*DPY+
     1               F(L0WCT+IXY+INDI)*DPZ
C
              GRPHC = (FENP - FP)/(DXC+tiny)
              IF(QGE(RUA,0.0)) THEN
                FUPW = F(L0UCT+IXY-INDI)*DPX+F(L0VCT+IXY-INDI)*DPY+
     1                 F(L0WCT+IXY-INDI)*DPZ
                GRPHU = (FP - FUPW)/(DXU+tiny)
              ELSE
                FUPW = F(L0UCT+IXY+2*INDI)*DPX+F(L0VCT+IXY+2*INDI)*DPY+
     1                 F(L0WCT+IXY+2*INDI)*DPZ
                GRPHU = (FUPW - FENP)/(DXU+tiny)
              ENDIF
C
              FPE = F(L0UCT+IXY)*DEX+F(L0VCT+IXY)*DEY+F(L0WCT+IXY)*DEZ
C
              GRPHCI = (FEN - FPE)/(DXC+tiny)
              IF(QGE(RUA,0.0)) THEN
                FUPW = F(L0UCT+IXY-INDI)*DEX+F(L0VCT+IXY-INDI)*DEY+
     1                 F(L0WCT+IXY-INDI)*DEZ
                GRPHUI = (FPE - FUPW)/(DXU+tiny)
              ELSE
                FUPW = F(L0UCT+IXY+2*INDI)*DEX+F(L0VCT+IXY+2*INDI)*DEY+
     1                 F(L0WCT+IXY+2*INDI)*DEZ
                GRPHUI = (FUPW - FEN)/(DXU+tiny)
              ENDIF
C
            ENDIF
C
C.... Use gradients to calculate higher-order correction
            IF(ISC==13) THEN
C.... Kappa schemes
              FLUXO = RUA*(KAPP(INDVAR)*GRPHC+KAPM(INDVAR)*GRPHU)*GDX
              FLUXI = FLUXO
              IF(BFC.AND.NONORT.AND.UVW) FLUXI=RUA*
     1          (KAPP(INDVAR)*GRPHCI+KAPM(INDVAR)*GRPHUI)*GDX
            ELSE
C.... Flux-limited positive schemes
              CALL GXFLPS(GRPHC,GRPHU,ISCHEM,PSI)
              FLUXO = RUA*PSI*GRPHU*GDX
              FLUXI = FLUXO
              IF(BFC.AND.NONORT.AND.UVW) THEN
                CALL GXFLPS(GRPHCI,GRPHUI,ISCHEM,PSI)
                FLUXI = RUA*PSI*GRPHUI*GDX
              ENDIF
            ENDIF
C
            F(L0VAL+IXY)      = F(L0VAL+IXY)      - FLUXO
            F(L0VAL+IXY+INDI) = F(L0VAL+IXY+INDI) + FLUXI
          ENDIF
  10  CONTINUE
      NAMSUB = 'gxhoen'
      if(debhcs) call banner(2,namsub,110817)
      END
C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
      SUBROUTINE GXHOHL(NDIR,ISC,L0DZW,L0DZG,L0AH)
C *********************************************************************
C     Purpose: to calculate higher-order correction for HIGH/LOW faces
C              on a single slab of cells
C                       low face      high face
C     Input:  NDIR  =       6            5
C     Called by: GXHOCS
C *********************************************************************
      INCLUDE 'farray'
      COMMON
     1       /IDATA/  NX,NY,NZ,LUPR1,IDAT1(116)
     1       /LDATA/  CARTES,LD1(2),ONEPHS,LD2(3),XCYCLE,LD3(10),STEADY,
     1                BFC,LD4(44),NONORT,LD5(19)
     1       /IGE/    IXF,IXL,IYF,IYL,IGE1(2),IGR,ISC1,IRUN,IZ,
     1                IGE2,ISWEEP,IGE3,INDVAR,IGE4(11)
     1       /IZLIM/  IZF,IZL
     1       /GENI/   NXNY,IGN1(8),NFM,IGN2(38),NPHI4,IGN3(11)
     1       /GHOCS/  L0PHI,L0VAL,L0RUA,L0RHO,L0RHOH,L0R,L0FPR,L0PRP,
     1                L0CVL,L0UCT,L0VCT,L0WCT,L0XP,L0YP,L0ZP,LBDXG,
     1                LBDYG,LBDZG,KAPP(150),KAPM(150),U1ORU2,V1ORV2,
     1                W1ORW2,UVW,STACON,ISCHEM,debhoc,debhcs
     1       /NAMFN/  NAMFUN,NAMSUB
     1       /RDATA/  TINY,RDAT(84)
C
      LOGICAL BLKSLD,DOCOR,U1ORU2,V1ORV2,W1ORW2,CARTES,UVW,LD1,XCYCLE,
     1        LD2,STEADY,BFC,LD3,NONORT,LD4,QGE,LD5,ONEPHS,STACON
      logical debhoc,debhcs
      REAL KAPP,KAPM
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB = 'GXHOHL'
C
C.... Set indices for cell geometry
      IF(.NOT.BFC) THEN
        LDVXY = L0DZW+IZ
        LDGXY = L0DZG+IZ
      ELSE
        LDZG = L0DZG
        LDZ  = L0DZW
      ENDIF
C
C.... Adjust indices by one slab if LOW face
      ISLAB=(IZ-1)*NXNY
      IF(NDIR==6) THEN
        ISLAB=(IZ-2)*NXNY
        IF(.NOT.BFC) THEN
          LDVXY = LDVXY-1
          LDGXY = LDGXY-1
        ELSE
          LDZG = LDZG-NXNY
          LDZ  = LDZ-NXNY
          L0AH = L0AH-NXNY
        ENDIF
      ENDIF
C
      INDI = NFM
      INDZ = 1
      IF(BFC) INDZ = NXNY
C
      IF(U1ORU2) THEN
        IXLP= IXL
        IXL = MIN(IXL,(NX-1))
      ELSEIF(V1ORV2) THEN
        IYLP= IYL
        IYL = MIN(IYL,(NY-1))
      ENDIF
C.... Cycle over the patch.
      DO 10 IX=IXF,IXL
        IADD = (IX-1)*NY
        DO 10 IY=IYF,IYL
          IXY = IY+IADD
C
C.... Check for blocked faces or solid cells within stencil
          DOCOR=.NOT.BLKSLD(W1ORW2,U1ORU2,V1ORV2,NXNY,NY,1,IXY+ISLAB,5)
C.... If no blocked/solid cells/faces within stencil then continue
          IF(DOCOR) THEN
C
C.... Set mass-flux coefficient for face
            RUA = F(L0RUA+IXY)
            IF(W1ORW2) THEN
              IF(STACON) THEN
                RUA = 0.5*(RUA + F(L0RUA+IXY+NXNY))
              ELSE
                IF(F(L0CVL+IXY)>0.0) THEN
                  GRO = F(L0RHO+IXY)
                  IF(.NOT.ONEPHS) GVF = F(L0R+IXY)
                ELSE
                  GRO = F(L0RHOH+IXY)
                  IF(.NOT.ONEPHS) GVF = F(L0R+IXY+INDI)
                ENDIF
                GUF = 0.5*(F(L0CVL+IXY) + F(L0CVL+IXY+INDI))
                GAF = F(L0AH+IXY)
                RUA = GRO*GUF*GAF
                IF(.NOT.ONEPHS) RUA=RUA*GVF
              ENDIF
            ELSEIF(U1ORU2) THEN
              RUA = 0.5*(RUA + F(L0RUA+IXY+NY))
            ELSEIF(V1ORV2) THEN
              RUA = 0.5*(RUA + F(L0RUA+IXY+1))
            ENDIF
C
C.... Values of convected variable on each side of face
            FP  = F(L0PHI+IXY)
            FH  = F(L0PHI+IXY+INDI)
C
C.... Calculate required cell-geometry data
            IF(BFC) THEN
              LDGXY = LDZG+IXY
              LDVXY = LDZ +IXY
            ENDIF
            IF(W1ORW2) THEN
              DZC = F(LDVXY+INDZ)
              IF(QGE(RUA,0.0)) THEN
                DZU = F(LDVXY)
                GDZ = 0.5*F(LDVXY+INDZ)
              ELSE
                DZU = F(LDVXY+2*INDZ)
                GDZ = -0.5*F(LDVXY+INDZ)
              ENDIF
            ELSE
              DZC = F(LDGXY)
              IF(QGE(RUA,0.0)) THEN
                DZU = F(LDGXY-INDZ)
                GDZ = 0.5*F(LDVXY)
              ELSE
                DZU = F(LDGXY+INDZ)
                GDZ = -0.5*F(LDVXY+INDZ)
              ENDIF
            ENDIF
C
C.... Calculate gradients of convected variable
            IF(.NOT.(BFC.AND.NONORT.AND.UVW)) THEN
C
              dzc=dzc+tiny
              dzu=dzu+tiny
              GRPHC = (FH-FP)/DZC
              IF(QGE(RUA,0.0)) THEN
                fpnei=F(L0PHI+IXY-INDI)
                GRPHU = (FP-F(L0PHI+IXY-INDI))/DZU
              ELSE
                fpnei= F(L0PHI+IXY+2*INDI)
                GRPHU = (F(L0PHI+IXY+2*INDI)-FH)/DZU
              ENDIF
C
            ELSE
C
              IF(W1ORW2) THEN
                INDB = NXNY
                DP = F(LBDZG+IXY)
              ELSEIF(U1ORU2) THEN
                INDB = NY
                DP = F(LBDXG+IXY)
              ELSEIF(V1ORV2) THEN
                INDB = 1
                DP = F(LBDYG+IXY)
              ENDIF
C
              dp=dp+tiny
              DPX = (F(L0XP+IXY+INDB)-F(L0XP+IXY))/DP
              DPY = (F(L0YP+IXY+INDB)-F(L0YP+IXY))/DP
              DPZ = (F(L0ZP+IXY+INDB)-F(L0ZP+IXY))/DP
C
              IF(NDIR==5) THEN
                FH = F(L0UCT+IXY+INDI)*DPX+F(L0VCT+IXY+INDI)*DPY+
     1               F(L0WCT+IXY+INDI)*DPZ
              ELSE
                FP = F(L0UCT+IXY)*DPX+F(L0VCT+IXY)*DPY+F(L0WCT+IXY)*DPZ
              ENDIF
C
              GRPHC = (FH - FP)/DZC
              IF(QGE(RUA,0.0)) THEN
                FUPW = F(L0UCT+IXY-INDI)*DPX+F(L0VCT+IXY-INDI)*DPY+
     1                 F(L0WCT+IXY-INDI)*DPZ
                GRPHU = (FP - FUPW)/DZU
              ELSE
                FUPW = F(L0UCT+IXY+2*INDI)*DPX+F(L0VCT+IXY+2*INDI)*DPY+
     1                 F(L0WCT+IXY+2*INDI)*DPZ
                GRPHU = (FUPW - FH)/DZU
              ENDIF
C
            ENDIF
C
C.... Use gradients to calculate higher-order correction
            IF(ISC==13) THEN
C.... Kappa schemes
              FLUX  = RUA*(KAPP(INDVAR)*GRPHC+KAPM(INDVAR)*GRPHU)*GDZ
            ELSE
C.... Flux-limited positive schemes
              CALL GXFLPS(GRPHC,GRPHU,ISCHEM,PSI)
              FLUX = RUA*PSI*GRPHU*GDZ
            ENDIF
C
            IF(NDIR==5) THEN
              F(L0VAL+IXY) = F(L0VAL+IXY) - FLUX
            ELSE
              F(L0VAL+IXY) = F(L0VAL+IXY) + FLUX
            ENDIF
          ENDIF
C
  10    CONTINUE
C
      IF(U1ORU2) THEN
        IXL = IXLP
      ELSEIF(V1ORV2) THEN
        IYL = IYLP
      ENDIF
C
      NAMSUB = 'gxhohl'
C
      END
C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
      SUBROUTINE GXFLPS(GRC,GRU,ISC,PSI)
C *********************************************************************
C     Purpose: to calculate flux-limiter functions for Flux-Limited
C              Positive Schemes. Ten limiters have been implemented:
C              MUSCL, van Leer harmonic, van Albada, Koren, SMART,
C              Minmod, OSPRE, H-Quick, Superbee and UMIST. Others may
C              be implemented by imitation.
C
C     Input:   GRC - gradient of convected variable across face
C              GRU - gradient of convected variable upwind of face
C              ISC - scheme number which indicates scheme selection
C
C     Output:   PSI - Flux-limiter function
C     Called by: GXHOEN, GXHOHL
C *********************************************************************
C
c.... initialisation in case of inappropriate entry
      PSI=0.0
      RP=(GRC+1.E-10)/(GRU+1.E-10)
C
      IF(ISC==6) THEN
C
C..SMART: PLNS limiter based on QUICK (Gaskell & Lau [1988])
C
        R1=AMIN1(0.25+0.75*RP,4.0)
        R1=AMIN1(R1,2.0*RP)
        PSI=AMAX1(0.0,R1)
C
      ELSEIF(ISC==7) THEN
C
C..Koren: PLNS limiter based on CUS (Koren [1993])
C
        R1=AMIN1(0.33333+0.66666*RP,2.0)
        R1=AMIN1(R1,2.0*RP)
        PSI=AMAX1(1.E-20,R1)
C
      ELSEIF(ISC==8) THEN
C
C..van Leer1/MUSCL : PLS limiter based on Fromm & identical to Noll'S
C                    MLU scheme (van Leer[1979])
C
        R1=AMIN1(0.5+0.5*RP,2.0)
        R1=AMIN1(R1,2.0*RP)
        PSI=AMAX1(1.E-20,R1)
C
      ELSEIF(ISC==9) THEN
C
C..H-QUICK: HNS limiter based on QUICK (Deconinck & Waterson [1995])
C
        IF(RP>0.0) THEN
          PSI=4.0*RP/(RP+3.0)
        ELSE
          PSI=1.0E-20
        ENDIF
C
      ELSEIF(ISC==10) THEN
C
C..OSPRE: PRS limiter (Deconinck & Waterson [1995])
C
        RP2=RP*RP
        PSI=1.5*(RP2+RP)/(RP2+RP+1.0)
C
      ELSEIF(ISC==11) THEN
C
C..van Leer2: PRS harmonic limiter based on Fromm, which is
C             identical to Zhu HLPA scheme  (see Hirsch[1990])
C
        IF(RP>0.0) THEN
          PSI=2.0*RP/(RP+1.0)
        ELSE
          PSI=1.0E-20
        ENDIF
C
      ELSEIF(ISC==12) THEN
C
C..van Albada: 'classical' limiter (van Albada et al [1982])
C
        PSI=RP*(RP+1.0)/(RP*RP+1.0)
C
      ELSEIF(ISC==13) THEN
C
C..Minmod: 'classical' limiter, which is identical to Zhu/Rodi SOUCUP
C          scheme  (see Roe [1986])
C
        RP=AMIN1(RP,1.0)
        PSI=AMAX1(0.0,RP)
C
      ELSEIF(ISC==14) THEN
C
C..Superbee: 'classical' limiter  (see Roe [1986])
C
        R1=AMIN1(2*RP,1.0)
        R2=AMIN1(RP,2.0)
        R2=AMAX1(R1,R2)
        PSI=AMAX1(0.0,R2)
C
      ELSEIF(ISC==15) THEN
C
C..UMIST: PLS limiter based mainly on QUICK Lien & Leschziner [1994]
C
        R1=AMIN1(0.75+0.25*RP,2.0)
        R2=AMIN1(0.25+0.75*RP,2.0*RP)
        R1=AMIN1(R1,R2)
        PSI=AMAX1(0.0,R1)
C
      ELSEIF(ISC==16) THEN
C
C..H-CUS: HNS limiter based on CUS (Deconinck & Waterson [1995])
C
        IF(RP>0.0) THEN
          PSI=3.0*RP/(RP+2.0)
        ELSE
          PSI=1.0E-20
        ENDIF
C
      ELSEIF(ISC==17) THEN
C
C..CHARM: PRNS limiter based on QUICK (Zhou [1995])
C
        IF(RP>0.0) THEN
          PSI=RP*(3.*RP+1.)/(RP*RP+2.0*RP+1.0)
        ELSE
          PSI=1.0E-20
        ENDIF
      ENDIF
C
      END
C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      FUNCTION BLKSLD(VELI,VELN1,VELN2,INDI,INDN1,INDN2,IJK,IDIR)
C *********************************************************************
C     Purpose: to detect presence of blocked faces or solid cells
C              within stencil used to calculate higher-order correction
C              for a single cell face. If BLKSLD is .TRUE. then no
C              higher-order correction will be calculated.
C
C     Input:   VELI  - .TRUE. if variable is "in-line" velocity
C              VELN1 - .TRUE. if variable is NOT "in-line" velocity
C              VELN2 - .TRUE. if variable is NOT "in-line" velocity
C              INDI  - Offset index in dirn. normal to face
C              INDN1 - Offset index corresponding to VELN1
C              INDN2 - Offset index corresponding to VELN2
C              IJK   - IY+(IX-1)*NY + (IZ-1)*NXNY
C              IDIR  -
C
C     Output:  BLKSLD = .TRUE. if blockage/solid detected
C     Called by: GXHOEN, GXHOHL
C *********************************************************************
C
      INCLUDE 'farray'
C
      LOGICAL BLKSLD,VELI,VELN1,VELN2,LSLDIR,LSOLID
C
      BLKSLD = LSOLID(IJK)
      IF(BLKSLD) RETURN
C.... Check for blocked faces
      BLKSLD = LSLDIR(IJK,IDIR)
      BLKSLD = BLKSLD.OR.LSLDIR(IJK+INDI,IDIR)
      BLKSLD = BLKSLD.OR.LSLDIR(IJK-INDI,IDIR)
      IF(BLKSLD) RETURN
C
      IF(VELI) THEN
        BLKSLD = LSLDIR(IJK+2*INDI,IDIR)
      ELSEIF(VELN1) THEN
        BLKSLD = LSLDIR(IJK+INDN1,IDIR)
        BLKSLD = BLKSLD.OR.LSLDIR(IJK+INDI+INDN1,IDIR)
        BLKSLD = BLKSLD.OR.LSLDIR(IJK-INDI+INDN1,IDIR)
      ELSEIF(VELN2) THEN
        BLKSLD = LSLDIR(IJK+INDN2,IDIR)
        BLKSLD = BLKSLD.OR.LSLDIR(IJK+INDI+INDN2,IDIR)
        BLKSLD = BLKSLD.OR.LSLDIR(IJK-INDI+INDN2,IDIR)
      ENDIF
C
      END
c