cc
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 'farray'
INCLUDE 'grdear'
INCLUDE 'grdloc'
INCLUDE '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 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE '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 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE '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