c
c
c     file name         GXPOTFLO.HTM     070405
C**** SUBROUTINE GXPOTC is called from GREX3, group 19, section 5,
C     and is entered only when POTVEL and POTCMP are set to T, and
C     only for  ISWEEP > 1.
C
C.... The library case 117 makes use of it.
C
      SUBROUTINE GXPOTC(EPOR,NPOR,HPOR,FACTOR,EXPNNT,NZ)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      COMMON/GENI/NXNY,IGFIL(59)
      COMMON/RDA1/DTFALS(150)
      INTEGER EPOR,NPOR,HPOR
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
C
      NAMSUB = 'GXPOTC'
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... Only the high porosity is modified here, and by reference
C     to the W velocity alone. If help is needed in generalising the
C     coding to take account of other directions, please contact CHAM
C     user support. See core library case 117 for exemplification.
C
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      IF(STORE(HPOR)) THEN
        IF(STORE(W1).AND.IZ.NE.NZ) THEN
          SETPOR=.TRUE.
          L0HPOR=L0F(HPOR)
          L0W1=L0F(W1)
c....                              under-relaxation factor
          IF(DTFALS(HPOR).LT.0.0) THEN
            FAC1=-DTFALS(HPOR)
          ELSE
            FAC1=1.0
          ENDIF
c....                              adjust HPOR
          DO 100 I=1,NXNY
            IF(F(L0HPOR+I).GT.1.E-6) THEN
              TERM=1.0 - FACTOR*F(L0W1+I)**2
              IF(TERM.GT.1.E-5) THEN
                TERM=TERM  ** EXPNNT
              ELSE
                CALL WRIT40('!! Velocity supersonic. Method invalid!!')
         CALL SET_ERR(230,'Error. Velocity supersonic. Method invalid.'
     *              ,1)
C                CALL EAROUT(1)
              ENDIF
              F(L0HPOR+I)=F(L0HPOR+I) + FAC1*(TERM - F(L0HPOR+I))
            ENDIF
  100     CONTINUE
        ENDIF
      ELSE
        CALL WRIT40('Store HPOR for POTVEL = T and NZ > 1    ')
        CALL SET_ERR(231,'Error. Store HPOR for POTVEL = T and NZ > 1.'
     *              ,1)
C        CALL EAROUT(1)
      ENDIF
      NAMSUB = 'gxpotc'
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C**** SUBROUTINE GXPOTV is called from GREX3, group 19, section 5,
C     and is entered only when POTVEL is set to T, and called only
C     for ISWEEP.GT.1
C
C     The logical functions XF and XF0 are used to detect whether the
C     x (ie N, E, or H) face is a solid fluid boundary or has a zero
C     porosity value.
C
C.... Core-library case 116 exemplifies the use of this subroutine.
C
      SUBROUTINE GXPOTV(EPOR,NPOR,HPOR,NZ,DZGG,BFC)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      COMMON/IDATA/NX,NY,IDFL(118)
      COMMON/GENI/NXNY,IGFIL1(59)
      COMMON/LDATA/LFIL1(7),XCYCLE,LFIL2(76)
      COMMON /NAMFN/NAMFUN,NAMSUB
      INTEGER EPOR,NPOR,HPOR
      CHARACTER*6 NAMFUN,NAMSUB
      LOGICAL BFC,XCYCLE,EF,NF,HF,EF0,NF0,HF0,EQZ,SLD
C
      NAMSUB = 'GXPOTV'
      LBPOT =  LBNAME('POT')
      L0POT=  L0F(LBPOT)
c..................................................... U1
      IF(STORE(U1) .AND. NX.GT.1) THEN
        L0U1=L0F(U1)
        IF(EPOR.NE.0) L0EPOR=L0F(EPOR)
        IF(BFC) THEN
          L0DXG=L0B(DXGPE) + (IZ-1)*NXNY
        ELSE
          L0DXG=L0F(DXG2D)
        ENDIF
        DO 10 IX=1,NX
          IPLUS=NY
          IF(IX.EQ.NX) IPLUS=IPLUS-NXNY
          DO 10 IY=1,NY
            I=IY+(IX-1)*NY
            IE=I+IPLUS
            F(L0U1+I)=(F(L0POT+I)-F(L0POT+IE))/F(L0DXG+I)
            IF(EF(I).OR.EF0(I)) THEN
              F(L0U1+I) = 0.0
            ELSEIF(EPOR.NE.0) THEN
              IF(EQZ(F(L0EPOR+I))) F(L0U1+I) = 0.0
            ENDIF
   10   CONTINUE
        IF(STORE(P1)) THEN
          L0P1=L0F(P1)
          DO IX=1,NX
            DO  IY=1,NY
              I=IY+(IX-1)*NY
              IW = I - NY
              IF(.NOT.XCYCLE) THEN
                IF(IX.EQ.1) THEN
                  F(L0P1+I)=F(L0P1+I) - 0.5*F(L0U1+I)**2
                ELSEIF(IX.EQ.NX) THEN
                  F(L0P1+I)=F(L0P1+I) - 0.5*F(L0U1+IW)**2
                ELSE
                  F(L0P1+I)=F(L0P1+I) - 0.125*(F(L0U1+I)+F(L0U1+IW))**2
                ENDIF
              ELSE
                IF(IX.EQ.1) IW=IW + NXNY
                F(L0P1+I)=F(L0P1+I) - 0.125*(F(L0U1+I)+F(L0U1+IW))**2
              ENDIF
              IF(SLD(I)) F(L0P1+I)=0.0
            ENDDO
          ENDDO
        ENDIF
      ENDIF
c
      IF(STORE(V1) .AND. NY.GT.1) THEN
        L0V1=L0F(V1)
        IF(NPOR.NE.0) L0NPOR=L0F(NPOR)
        IF(BFC) THEN
          L0DYG=L0B(DYGPN) + (IZ-1)*NXNY
        ELSE
          L0DYG=L0F(DYG2D)
        ENDIF
        DO 20 IX=1,NX
          DO 20 IY=1,NY
            I=IY+(IX-1)*NY
            IN=I+1
            F(L0V1+I)=(F(L0POT+I)-F(L0POT+IN))/F(L0DYG+I)
            IF(IY.EQ.NY.OR.NF(I).OR.NF0(I)) THEN
              F(L0V1+I) = 0.0
            ELSEIF(NPOR.NE.0) THEN
              IF(EQZ(F(L0NPOR+I))) F(L0V1+I) = 0.0
            ENDIF
   20   CONTINUE
        IF(STORE(P1)) THEN
          L0P1=L0F(P1)
          DO IX=1,NX
            DO  IY=1,NY
              I=IY+(IX-1)*NY
              IS=I-1
              IF(IY.EQ.1) THEN
                F(L0P1+I)=F(L0P1+I) - 0.5*F(L0V1+I)**2
              ELSEIF(IY.EQ.NY) THEN
                F(L0P1+I)=F(L0P1+I) - 0.5*F(L0V1+IS)**2
              ELSE
                F(L0P1+I)=F(L0P1+I) - 0.125*(F(L0V1+I) + F(L0V1+IS)) **2
              ENDIF
              IF(SLD(I)) F(L0P1+I)=0.0
            ENDDO
          ENDDO
        ENDIF
      ENDIF
c
      IF(STORE(W1) .AND. NZ.GT.1 .AND. IZ.NE.NZ) THEN
        L0POTH =  L0F(HIGH(LBPOT))
        L0W1=L0F(W1)
        IF(BFC) L0DZG=L0B(DZGPH) + (IZ-1)*NXNY
        IF(HPOR.NE.0) L0HPOR=L0F(HPOR)
        DO 30 I=1,NXNY
          IF(BFC) THEN
            DELZ=F(L0DZG+I)
          ELSE
            DELZ=DZG
          ENDIF
          F(L0W1+I)=(F(L0POT+I)-F(L0POTH+I))/DELZ
          IF(HF(I).OR.HF0(I)) THEN
            F(L0W1+I) = 0.0
          ELSEIF(HPOR.NE.0) THEN
            IF(EQZ(F(L0HPOR+I))) F(L0W1+I) = 0.0
          ENDIF
   30   CONTINUE
        IF(STORE(P1)) THEN
          L0P1=L0F(P1)
          L0LOW1=L0F(LOW(W1))
          DO IX=1,NX
            DO  IY=1,NY
              I=IY+(IX-1)*NY
              IF(IZSTEP.EQ.1) THEN
                F(L0P1+I)=F(L0P1+I) - 0.5*F(L0W1+I)**2
              ELSEIF(IZSTEP.EQ.NZ) THEN
                F(L0P1+I)=F(L0P1+I) - 0.5*F(L0LOW1+I)**2
              ELSE
                F(L0P1+I)=F(L0P1+I) - 0.125*(F(L0W1+I) + F(L0LOW1+I))**2
              ENDIF
              IF(SLD(I)) F(L0P1+I)=0.0
            ENDDO
          ENDDO
        ENDIF
      ENDIF
      NAMSUB = 'gxpotv'
      END
c
C**** SUBROUTINE GXPSIV is called from GREX3, group 19, section 5,
C     and is entered only when LPSIV is set to T, and called only
C     for ISWEEP.GT.1
C
      SUBROUTINE GXPSIV(BFC)
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      COMMON/IDATA/NX,NY,IDFL(118)
      COMMON/GENI/NXNY,IGFIL1(59)
C      COMMON/LDATA/LFIL1(7),XCYCLE,LFIL2(76)
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      LOGICAL BFC
C
      NAMSUB = 'GXPSIV'
      LBPSIV =  LBNAME('PSIV')
      IF(LBPSIV.NE.0) L0PSIV=  L0F(LBPSIV)
c..................................................... U1
      IF(STORE(U1) .AND. NX.GT.1) THEN
        L0U1=L0F(U1)
        IF(BFC) THEN
          L0DYG=L0B(DYGPN) + (IZ-1)*NXNY
        ELSE
          L0DYV=L0F(DYV2D)
        ENDIF
        DO IX=1,NX
          DO IY=1,NY
            I=IY+(IX-1)*NY
            IF(IY.EQ.1) THEN
              F(L0U1+I)=0.0
            ELSE
              F(L0U1+I)=( F(L0PSIV+I) - F(L0PSIV+I-1) ) / F(L0DYV+I)
            ENDIF
          ENDDO
        ENDDO
      ENDIF
c..................................................... V1
      IF(STORE(V1) .AND. NY.GT.1) THEN
        L0V1=L0F(V1)
        IF(BFC) THEN
          L0DXG=L0B(DXGPE) + (IZ-1)*NXNY
        ELSE
          L0DXU=L0F(DXU2D)
        ENDIF
        DO IX=1,NX
          DO IY=1,NY
            I=IY+(IX-1)*NY
            IF(IX.EQ.1)  THEN
              F(L0V1+I) = 0.0
            ELSE
              F(L0V1+I)=( F(L0PSIV+I-NY) - F(L0PSIV+I) ) / F(L0DXU+I)
            ENDIF
          ENDDO
        ENDDO
      ENDIF
c....................................................PSID for display
      LBPSID=LBNAME('PSID')
      IF(LBPSID.NE.0) THEN
        L0PSID=L0F(LBPSID)
        DO IX=2,NX
          DO IY=2,NY
            I=IY+(IX-1)*NY
            F(L0PSID+I)= 0.25*( F(L0PSIV+I) + F(L0PSIV+I-1)
     1                        + F(L0PSIV+I-NY) + F(L0PSIV+I-1-NY) )
          ENDDO
        ENDDO
      ENDIF
      namsub = 'grex3 '
      END
c