c

C
C***********************************************************************
C     The following is a summary of the variables used by the
C     scalar-equation method,  together with the PIL commands
C     necessary to activate the model.
C
C     * SURN
C     This variable is the fluid 'marker', of which the value
C     is one in a cell which is filled with the heavier fluid
C     only. A zero is assigned to a cell whose content is the
C     lighter fluid.  Fractional values are in cells in which
C     both fluids exist.
C
C     * VFOL
C     This variable stores the values of SURN at the old time
C     step, and provides information on the inflow volumetric
C     fluxes.
C
C     * IPRPSA and IPRPSB
C     These parameters are given PRPS indices of the heavier
C     and the lighter fluid respectively in the Q1 file.
C
C     * GALA
C     Logical variable GALA must be set .true. in the Q1 file
C
C     * TERMS command
C     This command is used to suppress all the built-in terms
C     in the SURN-equation. It takes the form:
C
C         TERMS(SURN,N,N,N,N,P,P)
C
C     * Time PATCH/COVAL
C     This is always needed,  and should be in the Q1 file in
C     the form of
C
C         PATCH(PNAME,CELL,1,NX,1,NY,1,NZ,1,LSTEP)
C         COVAL(PNAME,SURN,GRND,GRND)
C
C     * Inflow condition for variable SURN must always be set
C     in the form of a net flux, ie the product of SURN value
C     and inflow velocity. In addition, a COVAL statement for
C     VFOL must be specified in which the VALue is 1.00 / rho
C     where rho is equal to the density of the incoming fluid.
C
C     For more detailed description of SEM, users are referred
C     to the PHOENICS Encyclopaedia.
C
C     Definitions of FNxxs are given where they appear in the
C     code.  The FUNCTion entry of the Encyclopaedia provides
C     full explanations.
C
C***********************************************************************
c
      SUBROUTINE GXSURF
      INCLUDE 'farray'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE 'grdbfc'
C
      COMMON /INDAUX/L0ISL,L0IST,L0SL,L0ST,NSTO,NSOL,L0SLRS,L0TTRS,
     1        L0TTFL,L0TTLS,L0MAXC,L0MAXV,L0MINV,L0RATE,L0MAXI,L0NETT,
     1        L0TTFLM,L0TTFLX,L0TTFLY,L0TTFLZ
      COMMON/GENI/NXNY,IGFIL0(1),NXNYST,IGFIL1(6),NFM,IFGIL2(45),IPRPS,
     1                                                IGFIL3(4)
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6  NAMFUN,NAMSUB
      LOGICAL NPRP,QEQ,TCOVAL,NEZ
      SAVE TCOVAL,L0E,L0N,L0H,L0A,IVFOL,ISURN,L0X,L0Y,NPRP,L0VOLO
C
      NAMSUB='GXSURF'
C>>>>>>>>>>>>>>>>>>>>> Group  1 ------- Sections  3 <<<<<<<<<<<<<<<<
C     Logical initializations and storage allocations
C.......................................................................
      IF(IGR.EQ.1.AND.ISC.EQ.3) THEN
        USOURC = .TRUE.
        TCOVAL = .FALSE.
        IF(NX.GT.1) THEN
          CALL GXMAKE(L0E,NXNY,'CONE')
        ENDIF
        IF(NY.GT.1) THEN
          CALL GXMAKE(L0N,NXNY,'CONN')
        ENDIF
        IF(NZ.GT.1) THEN
          CALL GXMAKE(L0H,NXNY,'CONH')
          CALL GXMAKE(L0A,NXNY,'AREA')
        ENDIF
        IF(NEZ(W1AD)) CALL GXMAKE(L0VOLO,NXNY*NZ,'VOLO')
C>>>>>>>>>>>>>>>>>>>>> Group  1 ------- Section  2 <<<<<<<<<<<<<<<<<<<<<
C     F-array zero locations
C.......................................................................
      ELSEIF(IGR.EQ.1.AND.ISC.EQ.2) THEN
        IVFOL = LBNAME('VFOL')
        ISURN = LBNAME('SURN')
        IF(NX.GT.1) L0X = L0F(DXU2D)
        IF(NY.GT.1) L0Y = L0F(DYV2D)
C>>>>>>>>>>>>>>>>>>>>> Group  8 ------- Section  7 <<<<<<<<<<<<<<<<<<<<<
C     Volumetric source for GALA
C.......................................................................
      ELSEIF(IGR.EQ.8.AND.ISC.EQ.7) THEN
C>>>>>>>>>>>>>>>>>>>>> Group  8 ------- Section  8 <<<<<<<<<<<<<<<<<<<<<
C     Algebraic convective fluxes ndirca=1 for north
C                                 ndirca=3 for  east
C                                 ndirca=5 for  high
C.......................................................................
      ELSEIF(IGR.EQ.8.AND.ISC.EQ.8) THEN
        IF(INDVAR.NE.R1)  RETURN
        IF(NDIRCA.EQ.1) THEN
          CALL FN0(-L0N,LD5)
        ELSEIF(NDIRCA.EQ.3) THEN
          CALL FN0(-L0E,LD5)
        ELSEIF(NDIRCA.EQ.5) THEN
          CALL FN0(-L0H,LD5)
        ENDIF
C>>>>>>>>>>>>>>>>>>>>> Group  8 ------- Section 12 <<<<<<<<<<<<<<<<<<<<<
C     Compute convective and time fluxes of SURN in different directions
C.......................................................................
      ELSEIF(IGR.EQ.8.AND.ISC.EQ.12) THEN
        IF(INDVAR.NE.ISURN) RETURN
        L0S = L0F(LSU); L0C = L0F(IVFOL); III=ISL(ISURN)
C.... Calculate convective flux of SURN in y direction
        IF(NY.GT.1) THEN
          LV=L0F(V1)
          IF(BFC) L0Y=L0B(DHYPN)+NY*NX*(IZ-1)
          DO 800 I=1,NX
            IADD=(I-1)*NY
            DO 800 J=1,NY-1
              IND=J+IADD
              IF(F(L0N+IND).GE.0.0) THEN
                IDB=IND-1
                IDM=IND
                IDF=IND+1
                IF(J.EQ.1) IDB=IDM
              ELSE
                IDB=IND
                IDM=IND+1
                IDF=IDM+1
                IF(J.EQ.NY-1) IDF=IDM
              ENDIF
              DYF=F(L0Y+IDF)
              DYM=F(L0Y+IDM)
              DYB=F(L0Y+IDB)
              DPF=F(L0C+IDF)-F(L0C+IDM)
              DPB=F(L0C+IDM)-F(L0C+IDB)
              BAR=F(L0C+IDM)+PHIBAR(F(LV+IND),DPF,DPB,DYF,DYM,DYB,DT)
              F(L0TTFLY+III)=F(L0TTFLY+III)+ABS(F(L0N+IND)*BAR)
              F(L0S+IND)  =F(L0S+IND  )-F(L0N+IND)*BAR
  800         F(L0S+IND+1)=F(L0S+IND+1)+F(L0N+IND)*BAR
        ENDIF
C.... Calculate convective flux of SURN in x direction
        IF(NX.GT.1) THEN
          LU=L0F(U1)
          IF(BFC) L0X=L0B(DHXPE)+NY*NX*(IZ-1)
          JXM=NX-1
          IF(XCYCLE) JXM=NX
          DO 802 I=1,JXM
            IADD=(I-1)*NY
            DO 802 J=1,NY
              IND=J+IADD
              IF(F(L0E+IND).GE.0.0) THEN
                IDB=IND-NY
                IDM=IND
                IDF=IND+NY
                IF(I.EQ.1) IDB=IDM
                IF(XCYCLE) THEN
                  IF(I.EQ.1) IDB=IND+(NX-1)*NY
                  IF(I.EQ.NX) IDF=IND-IADD
                ENDIF
              ELSE
                IDB=IND
                IDM=IDB+NY
                IDF=IDM+NY
                IF(I.EQ.NX-1) IDF=IDM
                IF(XCYCLE) THEN
                  IF(I.EQ.NX) THEN
                    IDM=IDB-IADD
                    IDF=IDM+NY
                  ELSEIF(I.EQ.NX-1) THEN
                    IDF=IDB-IADD
                  ENDIF
                ENDIF
              ENDIF
              DXF=F(L0X+IDF)
              DXM=F(L0X+IDM)
              DXB=F(L0X+IDB)
              DPF=F(L0C+IDF)-F(L0C+IDM)
              DPB=F(L0C+IDM)-F(L0C+IDB)
              BAR=F(L0C+IDM)+PHIBAR(F(LU+IND),DPF,DPB,DXF,DXM,DXB,DT)
              F(L0S+IND)   =F(L0S+IND)   -F(L0E+IND)*BAR
              F(L0TTFLX+III)=F(L0TTFLX+III)+ABS(F(L0E+IND)*BAR)
              JNYINC=NY
              IF(XCYCLE.AND.I.EQ.NX) JNYINC=-IADD
  802     F(L0S+IND+JNYINC)=F(L0S+IND+JNYINC)+F(L0E+IND)*BAR
        ENDIF
C.... Calculate convective flux of SURN in z direction
        IF(NZ.GT.1) THEN
          LW=L0F(W1)
          IF(NEZ(W1AD)) THEN
            L0WGRID=L0F(WGRDNZ)
          ELSE
            L0WGRID=0
          ENDIF
          IVFOLL=L0F( LOW(IVFOL))
          IVFOLH=L0F(HIGH(IVFOL))
          JVFOLL=L0F(ANYZ(IVFOL,MIN(IZ+1,NZ)))
          JVFOLH=L0F(ANYZ(IVFOL,MIN(IZ+2,NZ)))
          IF(BFC) THEN
            LZ=L0B(DHZPH)
          ELSE
            LZ=L0F(DZWNZ)
          ENDIF
          DO 804 I=1,NX
            DO 804 J=1,NY
              IND=J+(I-1)*NY
              IF(F(L0H+IND).GE.0.0) THEN
                LCB=IVFOLL
                LCP=L0C
                LCF=IVFOLH
                IADD=ITWO(-NFM,0,IZ.GT.1)
                IADDZ=ITWO(-1,0,IZ.GT.1)
                IF(IZ.EQ.1) LCB=LCP
                IF(BFC) THEN
                  LZF=LZ+NXNY
                  LZP=LZ
                  LZB=LZ-NXNY
                  IF(IZ.EQ.1) LZB=LZ
                  DZF=F(LZF+IND)
                  DZP=F(LZP+IND)
                  DZB=F(LZB+IND)
                ELSE
                  DZF=F(LZ+MIN(IZ+1,NZ))
                  DZP=F(LZ+IZ)
                  DZB=F(LZ+MAX(IZ-1, 1))
                ENDIF
              ELSE
                LCB=L0C
                LCP=JVFOLL
                LCF=JVFOLH
                IF(IZ.EQ.NZ-1) LCF=LCP
                IADD=ITWO(NFM,0,IZ.LT.NZ)
                IADDZ=ITWO(1,0,IZ.GT.1)
                IF(BFC) THEN
                  LZB=LZ
                  LZP=LZ+NXNY
                  LZF=NXNY+LZP
                  IF(IZ.EQ.NZ-1) LZF=LZP
                  DZB=F(LZB+IND)
                  DZP=F(LZP+IND)
                  DZF=F(LZF+IND)
                ELSE
                  DZF=F(LZ+MIN(IZ+2,NZ))
                  DZP=F(LZ+MIN(IZ+1,NZ))
                  DZB=F(LZ+IZ)
                ENDIF
              ENDIF
              DPF=F(LCF+IND)-F(LCP+IND)
              DPB=F(LCP+IND)-F(LCB+IND)
              IF(L0WGRID.EQ.0) THEN
C                GWVEL=F(LW+IND)
                GWVEL=0.5*(F(LW+IND)+F(LW+IND+IADD))
              ELSE
C                GWVEL=F(LW+IND)-F(L0WGRID+IZ)
                GWVEL=0.5*(F(LW+IND)-F(L0WGRID+IZ)+
     1                     F(LW+IND+IADD)-F(L0WGRID+IZ+IADDZ))
              ENDIF
              BAR=F(LCP+IND)+PHIBAR(GWVEL,DPF,DPB,DZF,DZP,DZB,DT)
              F(L0TTFLZ+III)=F(L0TTFLZ+III)+ABS(F(L0H+IND)*BAR)
              IF(IZ.LT.NZ) F(L0S+IND)=F(L0S+IND)-F(L0H+IND)*BAR
              IF(IZ.GT. 1) F(L0S+IND)=F(L0S+IND)+F(L0A+IND)
              F(L0A+IND) = F(L0H+IND)*BAR
  804     CONTINUE
        ENDIF
Come    here to add the transient term to the equation when
C       this is not already done in Sec. 1 & 12 of Group 13.
        IF(.NOT.TCOVAL) THEN
          CALL FN34(LAP,VOL,1./DT)
          CALL FN53(LSU,VOL,IVFOL, 1./DT)
          CALL FN53(LSU,VOL,ISURN,-1./DT)
        ENDIF
C>>>>>>>>>>>>>>>>>>>>> Group 13 ------- Section  1 <<<<<<<<<<<<<<<<<<<<<
C     Set coefficient in source expression representing the omitted
C     transient term in SURM equation by TERMS(SURN,...) in Q1 file.
C.......................................................................
      ELSEIF(IGR.EQ.13.AND.ISC.EQ. 1) THEN
        TCOVAL = .TRUE.
        CALL FN2(CO,VOL,0.,1./DT)
C>>>>>>>>>>>>>>>>>>>>> Group 13 ------- Section  2 <<<<<<<<<<<<<<<<<<<<<
C     Set coefficient in fixed-pressure type of outflow boundary to
C     the local fluid density.
C.......................................................................
      ELSEIF(IGR.EQ.13.AND.ISC.EQ. 2) THEN
        IF(INDVAR.EQ.P1) CALL FN0(CO,DEN1)
C>>>>>>>>>>>>>>>>>>>>> Group 13 ------- Section 12 <<<<<<<<<<<<<<<<<<<<<
C     Set VAL in source for the omitted transient term in SURN equation
C     to SURN of previous time step, i.e. IVFOL.
C.......................................................................
      ELSEIF(IGR.EQ.13.AND.ISC.EQ.12) THEN
        IF(NEZ(W1AD)) THEN
          CALL FN56(VAL,IVFOL,-(L0VOLO+(IZ-1)*NXNY),VOL,1.0)
        ELSE
          CALL FN0(VAL,IVFOL)
        ENDIF
C>>>>>>>>>>>>>>>>>>>>> Group 13 ------- Section 13 <<<<<<<<<<<<<<<<<<<<<
C     Set outflow boundary condition for SURN.
C.......................................................................
      ELSEIF(IGR.EQ.13.AND.ISC.EQ.13) THEN
        IF(INDVAR.EQ.ISURN) THEN
          LITER(DEN1)=1
          CALL FN56(VAL,LMOUT1,ISURN,DEN1,-1.0)
          CALL GETPTC(NPATCH,GTYP,IIF,IIL,JJF,JJL,KKF,KKL,ITF,ITL)
          IF(QEQ(GTYP,2.0) .OR. QEQ(GTYP,3.0)) THEN
            CALL FN27(VAL,AEAST)
          ELSEIF(QEQ(GTYP,4.0) .OR. QEQ(GTYP,5.0)) THEN
            CALL FN27(VAL,ANORTH)
          ELSEIF(QEQ(GTYP,6.0) .OR. QEQ(GTYP,7.0)) THEN
            CALL FN27(VAL,AHIGH)
          ELSE
            CALL WRYT40('* Please Use Area Type For Exit Boundary')
            CALL WRYT40('   Where Outflow of Liquid is Expected *')
            CALL SET_ERR(246,'Error. See result file.',1)
C            CALL EAROUT(2)
          ENDIF
        ENDIF
C>>>>>>>>>>>>>>>>>>>>> Group 19 ------- Section  1 <<<<<<<<<<<<<<<<<<<<<
C.......................................................................
      ELSEIF(IGR.EQ.19.AND.ISC.EQ.1) THEN
        NPRP=.TRUE.
        IF(NEZ(W1AD)) THEN
          DO IZZ=1,NZ
            GDZ=F(L0F(DZWNZ)+IZZ)
            CALL FN2(-(L0VOLO+(IZZ-1)*NXNY), AHIGH, 0.0, GDZ)
          ENDDO
        ENDIF
 
C>>>>>>>>>>>>>>>>>>>>> Group 19 ------- Section  3 <<<<<<<<<<<<<<<<<<<<<
C     Compute material properties according to PRPS.
C.......................................................................
      ELSEIF(IGR.EQ.19.AND.ISC.EQ.3) THEN
        IF(NPRP) THEN
          CALL FN0(IVFOL,ISURN)
          CALL FNPRPS(IPRPS,IVFOL,RLOLIM,RUPLIM,IPRPSB,IPRPSA)
          IF(IZ.EQ.NZ) NPRP=.FALSE.
        ENDIF
      ENDIF
      NAMSUB='gxsurf'
      END
C***********************************************************************
      FUNCTION PHIBAR(VEL,DPF,DPB,DISF,DISM,DISB,DT)
C-----------------------------------------------------------------------
C.... This function returns the value of a scalar at a
C     cell face according to VAN LEER advection method
C-----------------------------------------------------------------------
      IF(VEL.GE.0.0) THEN
        DISS=+DISM
      ELSE
        DISS=-DISM
      ENDIF
      DPF=2.0*DPF/(DISM+DISF)
      DPB=2.0*DPB/(DISB+DISM)
C.... Compute phibar at cell face using VAN LEER scheme
C     Monotonicity conditions employed to limit dphi/dx
      S=0.0
      SGNPHI=1.0
      IF(DPF.LT.0.0) SGNPHI= -1.0
      IF(DPF*DPB.GT.0.0) S=SGNPHI
      DPHI=S*MIN(2.0*ABS(DPB),0.5*(ABS(DPB)+ABS(DPF)),2.0*ABS(DPF))
      PHIBAR=0.5*DPHI*(DISS-VEL*DT)
      END
c