c
C**** SUBROUTINE GXPORA is called from GREX3, group 19, section 3 and
C when IPORIA is set to greater than zero in SATELLITE.
C
C.... SETPOR must be set to T when the porosities are modified and
C reset, in order that EARTH knows that geometrical quantities must
C be re-calculated
C
C The examples provided are useful for simulating free-surface flow
C phenomena, with gravity operating in the negative z-direction.
C Surface waves are represented by allowing the porosities of the
C upper cells to increase with pressure.
C
C.... The library cases 457-458 exemplify the use of it.
C
C The following auxiliary subroutines are used:
C fn0(y,x) y = x
C fn10(y,x1,x2,a,b1,b2) y = a+b1*x1+b2*x2
C fn47(y,a,b) y = amax1(a,amin1(b,y))
C
SUBROUTINE GXPORA(SETPOR)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'satear'
COMMON /IGE/IFIL1(9),IZSTEP,IFIL2,ISWEEP,IFIL3(13)
LOGICAL SETPOR
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
NAMSUB = 'GXPORA'
C
IF(VPOR.NE.0) THEN
SETPOR=.FALSE.
IF(STORE(P2)) THEN
C.... Set ?POR = ?POR*PORIA + PORIB*(P1-P2)
CALL FN10(VPOR,P1,P2,PORIA,PORIB,-PORIB)
SETPOR = .TRUE.
ELSE
C.... Set ?POR = PORIA + PORIB * P1 for IPORIA < or = IZSTEP
IF(IZSTEP.LE.IPORIA) THEN
L0VPOR=L0F(VPOR)
L0P1=L0F(P1)
RELX=MAX(0.0,MIN(1.,ABS(DTFALS(VPOR))))
IF(ISWEEP==1) RELX=1.0
RELXM1=1.-RELX
DO 10 I=1,NX*NY
IF(F(L0VPOR+I).GT.1.E-6) THEN
VPR=F(L0VPOR+I)
F(L0VPOR+I)=RELXM1*VPR+RELX*(PORIA+PORIB*F(L0P1+I))
ENDIF
10 CONTINUE
SETPOR=.TRUE.
ENDIF
ENDIF
IF(SETPOR) THEN
CALL FN47(VPOR,VARMIN(VPOR),VARMAX(VPOR))
IF(EPOR.NE.0.AND.NX.GT.1) CALL FN0(EPOR,VPOR)
IF(NPOR.NE.0.AND.NY.GT.1) CALL FN0(NPOR,VPOR)
call setblk(izstep)
ENDIF
ENDIF
NAMSUB = 'gxpora'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c file name gxmovpor.htm 120704
c
C.... GXMOBS is called from GREX3 (Gr.19 Sc.1) if IPORIB.NE.0 from Group
c 19, section 1, i.e. at the start of each time step, in
c to compute values of VPOR, EPOR, NPOR or HPOR.
C
C The following examples are provided:
C - IPORIB= 1 creates porosities appropriate to a spherical
C coordinate system with the y-direction radial;
C - IPORIB= 2 resets VPOR to model moving piston.
C
C.... Variables PRBLOK and PRUNBL store indices of materials, which will
C be used by EARTH to fill PRPS-array in just-blocked cells (PRBLOK)
C and/or unblocked cells (PRUNBL). Just blocked cell is the cell for
C which VPOR is set to zero at the current time step.
C By default, PRBLOK= PORPRP; PRUNBL=1.0. The user MUST modify them,
C if other solids or fluids are used in his/her case.
C
SUBROUTINE GXMOBS
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
COMMON /MOVCMN/LMVCMN(5),PRBLOK,PRUNBL
C 1 /MOVPST/PTYP,IXPS,IYPS,IZPS,LXMOVE,LYMOVE,LZMOVE
COMMON /MOVPSTI/IXPS,IYPS,IZPS
1 /MOVPSTL/LXMOVE,LYMOVE,LZMOVE
1 /MOVPSTR/PTYP
1 /SLDPRP/SOLPRP,PORPRP,VACPRP
1 /GENI/ IGFIL2(55),IPRPS,IFIL1(4)
1 /NAMFN/NAMFUN,NAMSUB
LOGICAL LMVCMN,QEQ,LXMOVE,LYMOVE,LZMOVE
CHARACTER*6 NAMFUN,NAMSUB
SAVE IPF,IPL,JPF,JPL,KPF,KPL,ITPF,ITPL
C
NAMSUB = 'GXMOBS'
C.... Set SETPOR to true to activate appropriate resetting of arrays
C used by EARTH test functions
SETPOR= .TRUE.
C.... Set indices for materials, which will be used by EARTH to fill just
C blocked (PRBLOK) or/and unblocked (PRUNBL) cells.
IF(ISTEP.EQ.FSTEP) CALL SUB2R( PRBLOK,PORPRP, PRUNBL,-1.0 )
C.... PISTON: Get the direction of piston movement and the sizes of area
C which could be covered by piston:
IF(IPORIB.EQ.2) THEN ! if the piston-patch is set via spedat
IF(ISTEP.EQ.FSTEP) THEN ! determine initial position of piston
! here i, j and k signify ix, iy, iz
CALL GETPTC('PISTON',PTYP,IPF,IPL,JPF,JPL,KPF,KPL,ITPF,ITPL)
LXMOVE= QEQ(PTYP,2.) .OR. QEQ(PTYP,3.) ! deduce direction of motion
LYMOVE= QEQ(PTYP,4.) .OR. QEQ(PTYP,5.) ! from patch type, e.g. x
LZMOVE= QEQ(PTYP,6.) .OR. QEQ(PTYP,7.) ! if east or west, y if ..
ENDIF
C.... This example blocks one cell every other time step. The user is
C free to introduce other rule of blocking by changing the coding:
ISTPIS= ISTEP/2
IF(ISTPIS*2.EQ.ISTEP ) THEN
PSTPOR= 0.0
ISTPIS= ISTPIS-1
ELSE
PSTPOR= 0.5
ENDIF
CALL SUB3(IXPS,IPF-1, IYPS,JPF-1, IZPS,KPF-1)
IF(LXMOVE) THEN
IXPS= ITWO(IPF+ISTPIS,IPL-ISTPIS,QEQ(PTYP,2.))
ELSEIF(LYMOVE) THEN
IYPS= ITWO(JPF+ISTPIS,JPL-ISTPIS,QEQ(PTYP,4.))
ELSE
IZPS= ITWO(KPF+ISTPIS,KPL-ISTPIS,QEQ(PTYP,6.))
ENDIF
ENDIF
C.... Loop over whole domain to introduce new blockages:
CALL SUB2( L0VPOR,0, L0PRPS,0 )
CALL SUB3( L0EPOR,0, L0NPOR,0, L0HPOR,0 )
DO 70 IZZ= 1,NZ
C.... Get addresses of blockages for IZ-slab:
IF(EPOR.NE.0) L0EPOR= L0F(EPOR)
IF(NPOR.NE.0) L0NPOR= L0F(NPOR)
IF(HPOR.NE.0) L0HPOR= L0F(HPOR)
IF(VPOR.NE.0) L0VPOR= L0F(VPOR)
IF(IPRPS.NE.0) L0PRPS= L0F(IPRPS)
C.... IPORIB switches between various options to set porosities
IF(IPORIB.EQ.1) THEN
C.... Create porosities appropriate to a spherical coordinate system
C with the y-direction radial.
L0XG= L0F(LXXG)
L0YG= L0F(LYYG)
L0YV= L0F(LYYV)
DO 50 IX= 1,NX
SINXG = SIN(F(L0XG+IX))/(RINNER+YVLAST)
SINXU = SIN(F(L0XG+IX))/(RINNER+YVLAST)
L0ECEL= L0EPOR+NY*(IX-1)
L0NCEL= L0NPOR+NY*(IX-1)
L0HCEL= L0HPOR+NY*(IX-1)
L0VCEL= L0VPOR+NY*(IX-1)
IF(EPOR.NE.0) THEN
DO 10 IY= 1,NY
F(L0ECEL+IY)=(RINNER+F(L0YG+IY))*SINXU
10 CONTINUE
ENDIF
IF(NPOR.NE.0) THEN
DO 20 IY= 1,NY
F(L0NCEL+IY)=(RINNER+F(L0YV+IY))*SINXG
20 CONTINUE
ENDIF
IF(HPOR.NE.0) THEN
DO 30 IY= 1,NY
F(L0HCEL+IY)=(RINNER+F(L0YG+IY))*SINXG
30 CONTINUE
ENDIF
IF(VPOR.NE.0) THEN
DO 40 IY= 1,NY
F(L0VCEL+IY)=(RINNER+F(L0YG+IY))*SINXG
40 CONTINUE
ENDIF
50 CONTINUE
ELSEIF(IPORIB.EQ.2) THEN ! the moving piston
IF(ITPF.LE.ISTEP .AND. ISTEP.LE.ITPL) THEN
IF(KPF.LE.IZZ .AND. IZZ.LE.KPL) THEN
DO 60 IX= IPF,IPL
DO 60 IY= JPF,JPL
IJ= IY+(IX-1)*NY
IF(IX.LT.IXPS .OR. IY.LT.IYPS .OR. IZZ.LT.IZPS) THEN
F(L0VPOR+IJ)= 0.0
ELSEIF(IX.EQ.IXPS.OR.IY.EQ.IYPS.OR.IZZ.EQ.IZPS) THEN
F(L0VPOR+IJ)= PSTPOR
ENDIF
60 CONTINUE
ENDIF
ENDIF
ELSEIF(IPORIB.EQ.3) THEN
c.... north and volume porosities are proportional to IY**2
FNY=FLOAT(NY)
RFNY=1.0/FNY
TERM2=0.0
DO IY=1,NY
TERM1=(FLOAT(IY)*RFNY)**2
F(L0NPOR+IY)=TERM1
F(L0VPOR+IY)=0.5*(TERM1+TERM2)
TERM2=TERM1
ENDDO
ELSE
C.... The user can introduce his/her own option to set blockages here
ENDIF
C.... Reset addresses for use in Z-loop:
IF(NZ.GT.1) CALL RSTADR(IZZ)
70 CONTINUE
NAMSUB = 'gxmobs'
END
c