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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/satear'
      COMMON /IGE/IFIL1(9),IZSTEP,IFIL2(15)
      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)
            DO 10 I=1,NX*NY
              IF(F(L0VPOR+I).GT.1.E-6) F(L0VPOR+I)=PORIA+PORIB*F(L0P1+I)
   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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/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